Z<-QCF.D4F |Quadratic expansion as continued fraction 0 r #do "Z<-EDIT ""QCF.HLP""" #ifndef "\$QCF" \$QCF<-"" QM<-5 |Default Z<-ADD_ALIAS QCF.D4S "QCF.D4S" edqcf|WW.ED "qcf.d4f ",ARGS ldqcf|#copy"qcf.d4f" & QCF.D4F pell|RP_PELL 1 t #fi ARGS pelltab|LIST.PELL 2000 "QCF.HLP" >>EOF This set of functions evaluates the continued fraction expansion for the square roots of integers. The algorithm can be extended to work with multiple precision integers. A continued fraction looks like this sqrt(3) = 1+ 1 1 1 1 1 1 -- -- -- -- -- -- etc. 1 +2 +1 +2 +1 +2 Successive terms are 1+1/1=2, 1+1/(1+1/2)=5/3, 1+1/(1+1/(2+1/1))=1+1/(1+1/3)=1+3/4=7/4. Notice that 7*7-3*4*4=1. Continued fractions are important in working out integer solutions to the equation X*X-M*Y*Y = 1 for a given number M. This is called PELL's EQUATION. In modern times continued fractions have been used to try and factorise large numbers, and thereby crack RSA encryption. Functions QCF M continued fraction of square root of M LIST.CF N tabulation of square roots up to N LIST.CF A B tabulation of square roots of x: A<=x<=B LIST.PELL A,B tabulate Pell's Equation solutions for A<=M(0=Y<-(c/Y) #mod Z)/0 Z<-EA32 Y,Z Z<-X MULT Y |Multiply two values (a+b*sqrt(q))/c X<-X,(2=rX)/1 & Y<-Y,(2=rY)/1 Z<-,X[0 1] #outer * Y[0 1] Z<-(Z[0]+QM*Z[3]),Z[1]+Z[2] Z<-QREDUCE Z, X[2]*Y[2] Z<-QREDUCE X;H;D Z<-X % EA32 X[2], EA32 2tX Z<-INV X | X=(a,b,c) so want c/a+b*sqrt(q) = (c*a-b*sqrt(q))/(a*a-q*b*b) Z<-(X[2]*Z*1 _1), +/(1,-QM)*Z*Z<-2tX Z<-QREDUCE Z**Z[2] Z<-INT X;T |Integer part (a+b*sqrt(q))/c Z<-X[0]+ISQRT QM*X[1]*X[1] Z<-Z%X[2] Z<-FRAC X;S |Fractional part S<-INT X Z<-(X[0]-X[2]*S), X[1] Z<-QREDUCE Z,X[2] Z<-QCF M;A;J;QD;X;X0 |Continued fraction J<-0 & QD<-ISQRT QM<-M Z<-A<-INT 0 1 1 ->(M=QD*QD)/0 X0<-X<-FRAC 0 1 1 AA: X<-INV X Z<-Z,A<-INT X X<-FRAC X ->(~f/X=X0)/AA XX: Z<-LIST.CF N;I;J;K;T |List continued fractions for sqrt(k) 1(2=rN)/A2 |One arguement J<-i1+K<-ISQRT N Z<-Nr1 & Z[J*J]<-0 K<-0 & N<-rZ<-Z/irZ ->AA A2: |Range Z<-N[0]+i1+N[1]-N[0] J<-J*J<-i1+ISQRT N[1] K<-0 & N<-rZ<-(~Z e J)/Z Z AA: ->(K=N)/XX T<-QCF I<-Z[K] D.W "sqrt(",(zI),") = ", zT K<-K+1 & ->AA XX:Z<-"" Z<-EVAL_CF A;K;P;Q |Compute convergents of a[0]+1/a[1]+1/a[2]+1/a[3]... Z<-0 & ->(0=rA)/0 Z<-A[0],1 & ->(1=rA)/0 |Apply formula explained Hardy & Wright P<-A[0],1+A[0]*A[1] Q<-1,A[1] A<-2dA AA: ->(0=rA)/XX K<-1tA & A<-1dA P<-P,+/(1,K)*_2tP Q<-Q,+/(1,K)*_2tQ ->AA XX: Z<-P,:Q Z<-PELL_EQN M;A;D;P;Q;T;J;N |Solve Pell's equation. y*y-M*x*x=1 Z<-1 0 ->(1=N<-rA<-QCF M)/0 | M is a square J<-1 & P<-A[0],1+A[0]*A[1] & Q<-1,A[1] AA: Z<-(_1tP), _1tQ ->(1=+/(Z*Z)*1,-M)/XX T<-A[1+J #mod N-1] & J++ P<-(1dP),+/P*1,T Q<-(1dQ),+/Q*1,T ->AA XX: Z<-LIST.PELL N;I;V |Tabulate interval N<-((1=rN)/1),N Z<-1 #sed":d" & I<-1tN & N<-1 d N AA: ->(I=N)/XX Z<-RP_PELL I #print V<-("##### " z I), Z & V & I++ & ->AA XX: Z<-RP_PELL M;A;D;P0;P1;Q0;Q1;S;T;J;N |Solve Pell's equation. y*y-M*x*x=1 |Multiple precision Z<-"1 0" ->(1=N<-rA<-QCF M)/0 | M is a square J<-1 P0<-z A[0] & P1<-z 1+A[0]*A[1] Q0<-"1" & Q1<-z A[1] AA: Z<-P1," ",Q1 & T<-(P1*P1)-Q1*Q1*zM |(_48t"Z= ",Z), " d=",12tT |->(27=k)/XX ->("-1" #equiv T)/BB ->("1" #equiv T)/XX T<-z A[1+J #mod N-1] & J++ S<-P1 & P1<-P0+P1*T & P0<-S S<-Q1 & Q1<-Q0+Q1*T & Q0<-S ->AA BB: | Skip a few steps S<-"2"*P1*Q1 & P1<-(P1*P1)+Q1*Q1*zM Z<-P1," ",S XX: