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<B

EOF

Z<-D.W X
|Print in window, then show
#print X
Z<- #sed":s"

Z<-ISQRT N
|Integer part of #sqrt N
Z<-"N" #av #sqrt N

Z<-EA32 Y
|Euclidian Algorithm on Y[0 1]
Z<-f/Y<-#abs Y
->(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<k<=n
|List for A <= k <= B
->(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: