Mathematica code for the illustrations of the paper
Oliver Knill and Folkert Tangerman
Selfsimilarity and growth in Birkhoff sums for the golden rotation


This is code added for the revision   (http://www.math.harvard.edu/~kniil/birkhoffsum)

Figure 2   (the left graph  in the paper was generated with M=36, L=2000)

M=22;a=(Sqrt[5]-1)/2;r=Last[Convergents[a,M]];q=Denominator[r];g[x_]:=N[Log[2-2 Cos[2*Pi*x]]];
T[{k_,x_}]:=N[{k+1,x+g[k*a]}];b=NestList[T,{1,0},q];bb=Table[b[[k,2]]/(Log[k]*2),{k,2,Length[b]}];
L=500;ss=Table[0,{L}];Do[tt=Floor[bb[[k]]*L+1];If[0<tt≤L,ss[[tt]]++],{k,Length[bb]}];
SS=ListPlot[ss L/q,PlotStyle->{Thickness[0.001]},Joined->True,Filling->Axis,AspectRatio->1]

Figure 2 (the computation for the paper done with M=36, K=500  as a batch job. It takes a couple of hours and saves the coefficients in file)

M=20;(*Fibonnacci level*)K=300;(*number of Fourier coefficients*)PartialFraction[r_,m_]:=Module[{},U=ContinuedFraction[r,m];s=U[[m]];
Do[N[s=(1/s)+U[[kk]],400],{kk,m-1,1,-1}];s];alpha0=(Sqrt[5]-1)/2;alpha1=N[2 Pi alpha0];
qq=Table[Denominator[PartialFraction[alpha0,l]],{l,M}];pq=PartialFraction[alpha0,M];q=Denominator[pq];
T[k_]:=N[Log[2-2 Cos[k*alpha1]]];
ccoslist=ssinlist=Table[0,{M},{K}];
Do[u=0.0;ssin=0.0;ccos=0.0;Do[u+=T[k];v=N[Pi (u/Log[k]-1)];ccos+=N[Cos[l v]];ssin+=N[Sin[l v]];
If[MemberQ[qq,k],kk=Position[qq,k][[1,1]];ssinlist[[kk,l]]=ssin;ccoslist[[kk,l]]=ccos],{k,2,q}],{l,K}];
Save["fourierseries36.out",ssinlist,ccoslist];

Figure 2 (the Fourier coefficients)

<<fourierseries36.out;
M=Length[ccoslist];K=Length[ccoslist[[1]]];  abslist=Abs[ccoslist+I ssinlist];
PartialFraction[r_,m_]:=Module[{},U=ContinuedFraction[r,m];s=U[[m]];
Do[N[s=(1/s)+U[[kk]],400],{kk,m-1,1,-1}];s];alpha0=(Sqrt[5]-1)/2;
qq=Table[Denominator[PartialFraction[alpha0,l]],{l,M}];pq=PartialFraction[alpha0,M];q=Denominator[pq];
c[s_]:=Table[{k,Log[Abs[s[[M,k]]/q]]},{k,K}];A3=c[abslist];
h=Fit[A3,{1,x},x];T3:=Plot[h,{x,0,K}];
S0={Arrow[{{1,0},{K,0}}],Line[{{2,-10},{2,0}}]};
S0text={FontSize->12,Table[Text[k,{0.02,-k}],{k,1,10}]};
S0abs={FontSize->18,Text["log of fourier series",{K/3,-1}]};
S3=Graphics[{S0,S0abs,S0text,{Red,Line[A3]}}];
Show[{T3,S3},PlotRange->{{1,K},{-7,-2}},Axes->False,AspectRatio->1]

Figure (not shownThe first few Fourier coefficients. We computed up to M=42. )

S0={Arrow[{{0,0},{M+1,0}}],Arrow[{{0,-0.5},{0,1}}]};
S0a={Thickness[0.001],Table[Line[{{k,-0.02},{k,0.02}}],{k,1,M,3}]};
S0b={FontSize->12,Table[Text[k,{k,-0.02}],{k,1,M,3}]};
S0c={Thickness[0.001],Table[Line[{{-0.2,k/5},{0.2,k/5}}],{k,-1,4}]};
S0d={FontSize->12,Table[Text[N[k/10],{1,k/10}],{k,-1,7,2}]};
S0abs={FontSize->18,Text["fourier",{4.2,1}]};KK=10;
S3=Graphics[{S0,S0a,S0b,S0c,S0d,S0abs,Table[{Hue[n/KK],Line[Table[{l,abslist[[l,n]]/qq[[l]]},{l,3,M}]]},{n,KK}]}];
S=Show[S3,AspectRatio->1,PlotRange->{{-0.3,M+1.0},{-0.05,0.85}}]

July 31, 2011, OK

Figure (not shown in paper, the Fourier synthesis of the distribution)

Plot[1+(1/q)(2 Sum[(ccoslist[[M,n]] Cos[n x]+ssinlist[[M,n]] Sin[n x]),{n,K}]),{x,-Pi,Pi}]

The Hecke case again

Plot[g[x],{x,0,4Pi}]

a=(Sqrt[5]-1)/2;M=10;r=Last[Convergents[a,M]];q=Denominator[r];g[x_]:=N[Mod[x,2Pi]-Pi];
qq[m_]:=Denominator[Last[Convergents[a,m]]];
T[{k_,x_}]:={k+1,x+g[2Pi k*a]};S[{k_,x_}]:={k+1,x+g[2Pi k*r]};
b1=NestList[T,{1,0},q];b2=NestList[S,{1,0},q];R[s_]:=Table[s[[k,2]],{k,Length[s]}];
P=Graphics[Table[l=qq[k];{{RGBColor[1,0,0],PointSize[0.012],Point[{l,b1[[l,2]]}]},{RGBColor[0,0,1],PointSize[0.012],Point[{l+1,b1[[l+1,2]]}]}},{k,M}]];
S1=ListPlot[R[b1],Joined->True,PlotStyle->{Blue,Thickness[0.001]}];
S2=ListPlot[R[b2],Joined->True,PlotStyle->{Orange,Thickness[0.001]}];
S3=ListPlot[R[b2-b1],PlotStyle->{Green,Thickness[0.001]}];
Show[{S1,S2,S3,P},AspectRatio->1]

M=10;a=(Sqrt[5]-1)/2;q[m_]:=Denominator[Last[Convergents[a,m]]];g[x_]:=N[Mod[x,2Pi]-Pi];
T[{k_,x_}]:={k+1,x+g[2Pi k*a]};b=NestList[T,{1,0},q[M]];R[s_]:=Table[s[[k,2]],{k,Length[s]}];
P=Table[l=q[k];{{RGBColor[1,0,0],Point[{l,b[[l,2]]}]},{RGBColor[0,0,1],Point[{l+1,b[[l+1,2]]}]}},{k,M}];
Show[{ListPlot[R[b],Joined->True],Graphics[{PointSize[0.01],P}]}]

o. knill, August 5, 2011