%Hey. Contact redheffer at optoutofdatabases dot 33mail dot com if you have questions. This is very poorly documented. I know. Forgive me. Also, the coding style is terrible and I'm wanted in New Jersey for fraud.

%Really, every function here should be in a separate file, or they should at least have an end call. I didn't bother. I'm just trying to package everything together because I'm a lazy piece of shit. Sorry.

function [g]=RedhefferQuick(y,s)
%Does an implicit multiplication by a Redheffer matrix.
sc=size(y,1);
for i=1:s
su=0;
j=0;
while i*j if(i*j>s)
su=su+y(i*j-s);
end;
j=j+1;
end;
g(i,1)=su;
end;


function [d]=redcalcQuickSimple(n)
%Needed because I needed a more computationally efficient version of the algorithm.
y=[0;ones(n-1,1)];
N=floor(log(n)/log(2));
s=n; %s is the size of the non-identity block of Redheffer, sc is its complement.
for i=1:N
sc=s;
sc=s-floor(s/2);
s=floor(s/2);
g=RedhefferQuick(y(s+1:s+sc),s);
v=y(1:s)-g;
y=v;
end;
d=y+1;

function [h]=harmTransMultiply(z,i,n)
%Start at the ith column of the Harmonic matrix of total size n.
m=size(z,1);
h=zeros(n-i+1,1);
for j=1:m
for k=1:n-i+1
%Column i+j of the transposed Harmonic, row k+i
if(k>n)
else
h(k,1)=h(k,1)+z(j,1)*floor((i-1+k)/(i-1+j));
end;
end;
end;

function [h]=harmTransMultiplyContinuous(z,i,n)
%Start at the ith column of the Harmonic matrix of total size n.
m=size(z,1);
h=zeros(n-i+1,1);
for j=1:m
for k=1:n-i+1
%Column i+j of the transposed Harmonic, row k+i
if(k>n)
else
h(k,1)=h(k,1)+z(j,1)*((i-1+k)/(i-1+j));
end;
end;
end;



function [b,a]=harmTrans(n)
% n is the Merten's whatever you want to calculate, a is the calculation per the algorithm, and b is the calculation per the modified algorithm split between a positive and negative vector--this is production code meant to doublecheck/have a sanity check that I'm doing this right.
% Lines with "%%" are for debugging and exploration. These lines and the variables they give may be deleted.
% At the moment, I've been playing with the 'continuous' version. This is theoretically incorrect, I'm just doing this to mess around atm.
[a]=redcalcQuickSimple(n);%%
zp=ones(n,1);
zn=zeros(n,1);
zptemp=zeros(n,1);
zntemp=zeros(n,1);
N=log2(n);
c=0;
i=1;
sum=0;
while c d=2^(i-1);
d=min(n,c+d)-c;
% %[min(zp(c+1:c+d)),max(zp(c+1:c+d))]
% %[log2(min(zp(c+1:c+d))),log2(max(zp(c+1:c+d)))]
m1=harmTransMultiply(zn(c+1:c+d),c+1,n);
m1=harmTransMultiplyContinuous(zn(c+1:c+d),c+1,n);
m1=m1(1+d:end,:);
% m2=harmTransMultiply(zp(c+1:c+d),c+2,n);
% m2=m2(d:end,:);
m2=harmTransMultiply(zn(c+1:c+d),c+2,n);
m2=harmTransMultiplyContinuous(zn(c+1:c+d),c+2,n);
m2=m2(d:end,:);
rx1n=c+d+1;
rx2n=c+d+5;
ry1n=zn(rx1n);
ry2n=zn(rx2n);
rmn=(ry2n-ry1n)/(rx2n-rx1n);
rbn=ry2n-rmn*rx2n;
%if(i>1)
%sum=sum+norm(m1-m2,'inf')/(n*rmn)
%end;

%l11=1;
%l12=norm(m1-m2,'inf')/norm(zn,'inf');
zptemp(c+d+1:end)=zp(c+d+1:end)+m1-m2;

%h1=0;
%for k=1:n-(c+1)+1
% %Column i+j of the transposed Harmonic, row k+i
% if(k>n)
% else
% h1(k,1)=floor((c+1-1+k)/(c+1-1+1));
% end;
%end;
%h2=0;
%for k=1:n-(c+2)+1
% %Column i+j of the transposed Harmonic, row k+i
% if(k>n)
% else
% h2(k,1)=floor((c+2-1+k)/(c+2-1+1));
% end;
%end;
%[m1-m2,m1,m2]
%aa1=1;
%aa2=norm(m1-m2,'inf')/norm(zn,'inf');
%% if(size(zp(c+d+1:end),1)==0)
%% else
%% AP(i,1:13)=[i,c,d,max(zn(c+1:c+d)),min(zn(c+1:c+d)),max(zp(c+1:c+d)),min(zp(c+1:c+d)),max(zp(c+d+1:end)),min(zp(c+d+1:end)),max(m1),min(m1),max(m2),min(m2)];%%
%% end;
% %[norm(zp(c+d+1:end),1),norm(m1,1),norm(m2,1)]
% %[norm(zp(c+d+1:end),'inf'),norm(m1,'inf'),norm(m2,'inf')]
m1=harmTransMultiply(zp(c+1:c+d),c+1,n);
m1=harmTransMultiplyContinuous(zp(c+1:c+d),c+1,n);
m1=m1(1+d:end,:);
m2=harmTransMultiply(zp(c+1:c+d),c+2,n);
m2=harmTransMultiplyContinuous(zp(c+1:c+d),c+2,n);
m2=m2(d:end,:);

%l22=1;
%l21=norm(m1-m2,'inf')/norm(zp,'inf');

zntemp(c+d+1:end)=zn(c+d+1:end)+m1-m2;
%m1-m2<0
%aa3=norm(m1-m2,'inf')/norm(zp,'inf');
%aa4=1;
%[aa1,aa2;aa3,aa4]*[1;-1]
%% if(size(zp(c+d+1:end),1)==0)
%% else
%% AN(i,1:13)=[i,c,d,max(zp(c+1:c+d)),min(zp(c+1:c+d)),max(zn(c+1:c+d)),min(zn(c+1:c+d)),max(zn(c+d+1:end)),min(zn(c+d+1:end)),max(m1),min(m1),max(m2),min(m2)];%%
%% end;
%size(zptemp(c+d+1:end))
%%plot(zptemp(c+d+1:end))
%%plot(h1)
%%plot(h1(2:end)-h2)

rx1p=c+d+1;
rx2p=c+d+5;
ry1p=zp(rx1p);
ry2p=zp(rx2p);
rmp=(ry2p-ry1p)/(rx2p-rx1p);
rbp=ry2p-rmp*rx2p;
[rmn,rmp]
[rbp,rbn]
%for ii=1:n
% cp(ii)=rmp*ii+rbp;
% cn(ii)=rmn*ii+rbn;
%end;
hold on
%plot(cp/(n*(1+sqrt(2)/2)^(i)),"color",[0,0,0])
%plot(cn/(n*(1+sqrt(2)/2)^(i)),"color",[0,0,0])
plot(zp/(n*(1+sqrt(2)/2)^(i)),"color",[0.1,0.1+i/(2*N),0.5])
plot(zn/(n*(1+sqrt(2)/2)^(i)),"color",[0.5,0.1+i/(2*N),0.1])
%norm(cp(n-10:end)'-zp(n-10:end),2)/norm(zp,'inf')
%norm(cn(n-10:end)'-zn(n-10:end),2)/norm(zn,'inf')
%rbp
%rbn

%[norm(zp,'inf'),norm(zn,'inf')]
%[l11,l12;l21,l22]
%l12-l21

hold off
%max(abs((zptemp(c+d+1:end)-zntemp(c+d+1:end))/(max(max([zptemp(c+d+1:end),zntemp(c+d+1:end)])))))

zp=zptemp;
zn=zntemp;
% %norm(zn(c+1:c+d)-zp(c+1:c+d),1)
% %[zp(end)-zn(end),norm(zp(c+d+1:end)-zn(c+d+1:end),1)]
%% blep=zp(c+1:c+d)-zn(c+1:c+d);
% %[zp(c+1:c+d)-zn(c+1:c+d),zp(c+1:c+d),zn(c+1:c+d)]
%% [max(blep),min(blep)]
c=c+d;
i=i+1;
end;
b=zp(end,1)-zn(end,1);