Continued Fractions<Text-field layout="Section" style="Section">Finite continued fractions</Text-field>ourcontinuedfraction:=proc(x)
local xprime,
a;
xprime:=x;
a:=[];
while not type(xprime,integer) do
a:=[op(a),floor(xprime)];
xprime:=1/(xprime-floor(xprime));
end do;
a:=[op(a),floor(xprime)];
end proc:ourcontinuedfraction(125/54);ourfromcontinuedfraction:=proc(a)
local c;
c:=a;
while nops(c)>1 do
c[-2]:=c[-2]+1/c[-1];
c:=c[1..-2];
end do;
c[1];
end proc:ourfromcontinuedfraction([2,3,5,1,2]);ourconvergents:=proc(x)
local p, q, a, c, i;
if type(x,list)
then
a:=x;
else
a:=ourcontinuedfraction(x);
end if;
p[0]:=a[1];
q[0]:= 1;
c[0]:= p[0]/q[0];
if nops(a) > 0 then
p[1]:=a[1]*a[2] + 1;
q[1]:=a[2];
c[1]:=p[1]/q[1];
for i from 2 to nops(a)-1 do
p[i]:=a[i+1]*p[i-1]+p[i-2];
q[i]:=a[i+1]*q[i-1]+q[i-2];
c[i]:=p[i]/q[i];
end do;
end if;
[seq(c[i],i=0..nops(a)-1)];
end proc:
ourconvergents(125/54);ourconvergents([2,3,5,1,2]);<Text-field layout="Section" style="Section">Infinite continued fractions</Text-field>ourcontinuedfraction2:=proc(x,n)
local xprime, a, i;
xprime:=x;
a:=[floor(xprime)];
i:=1;
while i<n and not type(xprime,integer) do
xprime:=1/(xprime-floor(xprime));
a:=[op(a),floor(xprime)];
i:=i+1;
end do;
a;
end proc:ourcontinuedfraction2(sqrt(3),6);ourcontinuedfraction2(125/54,8);ourconvergents2:=(x,n)->ourconvergents(ourfromcontinuedfraction(ourcontinuedfraction2(x,n))):
ourconvergents2(125/54,8);ourconvergents2(sqrt(3),8);<Text-field layout="Section" style="Section">Rational approximation of real numbers</Text-field>Digits:=20;ourconvergents2(evalf(Pi),5);array(map(n->[ourconvergents2(Pi,5)[n],evalf(abs(Pi-ourconvergents2(Pi,5)[n]))],[seq(n,n=1..5)]));for b from 1 to 106 do
a:=round(b*Pi):
if abs(evalf(Pi)-a/b) < abs(evalf(Pi)-22/7) then
print([a/b, evalf(abs(Pi-a/b)), evalf(abs(b*Pi-a))]);
end if
end do:<Text-field layout="Section" style="Section">Periodic continued fractions</Text-field>quadirrcontinuedfraction:=proc(x,y,d)
local n, alpha, a, P, Q, ourlist, position, s, k;
n:=0;
P[n]:=x;
Q[n]:=y;
ourlist:=[];
while not member([P[n],Q[n]],ourlist,'position') do
alpha[n]:=(P[n]+sqrt(d))/Q[n];
ourlist:=[op(ourlist),[P[n],Q[n]]];
a[n]:=floor(alpha[n]);
P[n+1]:=a[n]*Q[n]-P[n];
Q[n+1]:=(d-(P[n+1])^2)/Q[n];
n:=n+1;
end do;
s:=n-position+1;
[seq(a[k],k=0..position-2),[seq(a[k],k=position-1..position+s-2)]];
end proc:quadirrcontinuedfraction(10,7,11);<Text-field layout="Section" style="Section">Continued fraction factorization</Text-field>fermatfactor:=proc(d)
local x;
x:=ceil(sqrt(d));
while not sqrt(x^2-d)=floor(sqrt(x^2-d)) do
x:=x+1;
end do;
[x+sqrt(x^2-d),x-sqrt(x^2-d)];
end proc: fermatfactor(320813);fermatfactor(13);cfcongruences:=proc(d,k)
local n, alpha, a, p, P, Q;
P[0]:=0;
Q[0]:=1;
p[-1]:=1;
p[0]:=floor(sqrt(d));
for n from 0 to k do
P[n+1]:=(floor((P[n]+sqrt(d))/Q[n]))*Q[n]-P[n];
Q[n+1]:=(d-(P[n+1])^2)/Q[n];
p[n+1]:=floor((P[n+1]+sqrt(d))/Q[n+1])*p[n]+p[n-1];
print([n,modp(p[n],d),ifactor((-1)^(n-1)*Q[n+1])]);
end do;
end proc:cfcongruences(320813,15);