62 lines
1.7 KiB
Text
62 lines
1.7 KiB
Text
with javascript_semantics
|
|
|
|
include mpfr.e
|
|
|
|
procedure fun(mpz a,b,t, integer c)
|
|
-- {a,b} = {b,c*b+a} (and t gets trashed)
|
|
mpz_set(t,a)
|
|
mpz_set(a,b)
|
|
mpz_mul_si(b,b,c)
|
|
mpz_add(b,b,t)
|
|
end procedure
|
|
|
|
function SolvePell(integer n)
|
|
integer x = floor(sqrt(n)), y = x, z = 1, r = x*2
|
|
mpz e1 = mpz_init(1), e2 = mpz_init(),
|
|
f1 = mpz_init(), f2 = mpz_init(1),
|
|
t = mpz_init(0), u = mpz_init(),
|
|
a = mpz_init(1), b = mpz_init(0)
|
|
if x*x!=n then
|
|
while mpz_cmp_si(t,1)!=0 do
|
|
y = r*z - y
|
|
z = floor((n-y*y)/z)
|
|
r = floor((x+y)/z)
|
|
fun(e1,e2,t,r) -- {e1,e2} = {e2,r*e2+e1}
|
|
fun(f1,f2,t,r) -- {f1,f2} = {f2,r*r2+f1}
|
|
mpz_set(a,f2)
|
|
mpz_set(b,e2)
|
|
fun(b,a,t,x) -- {b,a} = {f2,x*f2+e2}
|
|
mpz_mul(t,a,a)
|
|
mpz_mul_si(u,b,n)
|
|
mpz_mul(u,u,b)
|
|
mpz_sub(t,t,u) -- t = a^2-n*b^2
|
|
end while
|
|
end if
|
|
return {a, b}
|
|
end function
|
|
|
|
function split_into_chunks(string x, integer one, rest)
|
|
sequence res = {x[1..one]}
|
|
x = x[one+1..$]
|
|
integer l = length(x)
|
|
while l do
|
|
integer k = min(l,rest)
|
|
res = append(res,x[1..k])
|
|
x = x[k+1..$]
|
|
l -= k
|
|
end while
|
|
return join(res,"\n"&repeat(' ',29))&"\n"&repeat(' ',17)
|
|
end function
|
|
|
|
sequence ns = {4, 61, 109, 181, 277, 8941}
|
|
for i=1 to length(ns) do
|
|
integer n = ns[i]
|
|
mpz {x, y} = SolvePell(n)
|
|
string xs = mpz_get_str(x,comma_fill:=true),
|
|
ys = mpz_get_str(y,comma_fill:=true)
|
|
if length(xs)>97 then
|
|
xs = split_into_chunks(xs,98,96)
|
|
ys = split_into_chunks(ys,99,96)
|
|
end if
|
|
printf(1,"x^2 - %3d*y^2 = 1 for x = %27s and y = %25s\n", {n, xs, ys})
|
|
end for
|