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