Diophantine Equations<Text-field layout="Section" style="Section">Gaussian integers</Text-field>numtheory[GIgcd](10+2*I,8-I);ourgcd:=proc(w,z)
local a,x;
a:=[w,z];
while not a[2]=0 do
a:=[a[2],a[1]-a[2]* (round(Re(a[1]/a[2]))+I*round(Im(a[1]/a[2])))];
end do;
op(select(x->evalb(Re(x)>0 and Im(x)>=0),{a[1],-a[1],I*a[1],-I*a[1]}));
end proc:ourgcd(10+2*I,8-I);p:=541:
g:=2:
modp(g&^((p-1)/2),p);x:=modp(g&^((p-1)/4),p);numtheory[GIgcd](p,x+I);<Text-field layout="Section" style="Section">Pell's equation</Text-field>pellxy:=proc(n)
option remember:
if n=0 then
return [1,0]
else
return [2*pellxy(n-1)[1]+3*pellxy(n-1)[2],pellxy(n-1)[1]+2*pellxy(n-1)[2]]
end if;
end proc:[seq(pellxy(n),n=0..10)];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:<Text-field layout="Heading 1" style="Heading 1">Continued fraction solution of Pell's equation</Text-field>quadirrcontinuedfraction(0,1,13);numtheory[cfrac]([3,1,1,1,1,6,1,1,1,1]);<Text-field layout="Heading 1" style="Heading 1">The <Font italic="true">abc</Font><Font italic="false"> conjecture</Font></Text-field>for b from 1 to 50 do
for a from 1 to b do
if gcd(a,b)=1 then
producttemp:=product(ifactors(a*b*(a+b))[2][i][1],
i=1..nops(ifactors(a*b*(a+b))[2]));
if a+b>producttemp then print(a,b,a+b,[producttemp]);
end if
end if;
end do;
end do;