/*bc program davison */ define positivity(a,b,c,d){ if(a>=0 && b>=0 && c>=0 && d>=0){ return(1) }else{ return(0) } } /* nprod(l,m) finds the least n such that A_n and D=A_0...A_n are non-negative. * The matrix of globabl variables D=[globala,globalb,globalc,globald] is * returned along with n. See Proposition 4.1, J.L. Davison, 'An algorithm for * the continued fraction of e^{l/m}', Proceedings of the Eighth Manitoba * Conference on Numerical Mathematics and Computing (Univ. Manitoba, Winnipeg, * 1978), 169--179, Congress. Numer., XXII, Utilitas Math. */ define nprod(l,m){ auto k,a1,b1,c1,d1,temp1,temp2,s,t,t1,t2 globala=m+l;globalb=m;globalc=m;globald=m-l for(k=1;1;k++){ s=positivity(globala,globalb,globalc,globald) if(s){ break } t=(2*k+1)*m t1=globala+globalb t2=globalc+globald temp1=t*t1 temp2=t*t2 a1=temp1+globala*l b1=temp1-globalb*l c1=temp2+globalc*l d1=temp2-globald*l globala=a1 globalb=b1 globalc=c1 globald=d1 } return(k-1) } define abs(n){ if(n>=0) return(n) return(-n) } /* gcd(m,n) for any integers m and n */ /* Euclid's division algorithm is used. */ /* We use gcd(m,n)=gcd(|m|,|n|) */ define gcd(m,n){ auto a,b,c a=abs(m) if(n==0) return(a) b=abs(n) c=a%b while(c>0){ a=b b=c c=a%b } return(b) } define reduce(a,b,c,d){ auto g g=gcd(a,b) g=gcd(g,c) g=gcd(g,d) return(g) } /* * Input: a non-singular matrix A=[p,q;r,s], p,q,r,s>=0, A!=I_2, A!=[0,1;1,0]. * With L=[1,0;1,1] and R=[1,1;0,1], we express A uniquely as * a product of non-negative powers of L and R, followed by a row-balanced B. * B=[a,b;c,d] is row-balanced if (ad) or (cb) * and a,b,c>=0. We exclude A=I_2 and A=[0,1;1,0]. * See 'On continued fractions and finite automata', G.N. Raney, * Math. Annalen, 206, 265-283 (1973). */ define raney(p,q,r,s){ auto i,j,k k=0 while(1){ i=0 while(p>=r && q>=s){ p=p-r q=q-s i=i+1 } if(i){ globalr=globalr+i flagr=1 if(flagl*flagr){ flagl=0 print globall,"," globall=0 count=count+1 } k=k+1 } j=0 while(r>=p && s>=q){ r=r-p s=s-q j=j+1 } if(j){ globall=globall+j flagl=1 if(flagr*flagl){ flagr=0 print globalr,"," globalr=0 count=count+1 } k=k+1 } if((ps) || (p>r && s>q)){ break } } globala=p globalb=q globalc=r globald=s return(k) } /* We perform the algorithm of J.L. Davison's paper. * With n>=0, we first find the n* of Davison's Proposition 4.1 * and apply raney's factorisation to A_0...A_k, for n*<=k<=n*+n. * The number (count) of partial quotients of e^{l/m} found is returned. * count becomes positive for all large n. */ define davison(l,m,n){ auto g,i,j,k,t count=0 flagr=0 flagl=0 globalr=0 globall=0 k=nprod(l,m) g=reduce(globala,globalb,globalc,globald) if(g>1){ globala=globala/g globalb=globalb/g globalc=globalc/g globald=globald/g } i=k j=k+n while(i<=j){ if(i>k){ t=(2*i+1)*m t1=globala+globalb t2=globalc+globald temp1=t*t1 temp2=t*t2 a1=temp1+globala*l b1=temp1-globalb*l c1=temp2+globalc*l d1=temp2-globald*l globala=a1 globalb=b1 globalc=c1 globald=d1 } i=i+1 t=raney(globala,globalb,globalc,globald); /* the input matrix will not be I_2 or [0,1;1,0] here. */ g=reduce(globala,globalb,globalc,globald) if(g>1){ globala=globala/g globalb=globalb/g globalc=globalc/g globald=globald/g } } print "...\nThe number of partial quotients found for e^(",l,"/",m,") is " return(count) }