(phixonline)-->
with javascript_semantics
include mpfr.e
function ts(string ns, ps)
mpz n = mpz_init(ns),
p = mpz_init(ps),
t = mpz_init(),
r = mpz_init(),
pm1 = mpz_init(),
pm2 = mpz_init()
mpz_sub_ui(pm1,p,1) -- pm1 = p-1
mpz_fdiv_q_2exp(pm2,pm1,1) -- pm2 = pm1/2
mpz_powm(t,n,pm2,p) -- t = mod(n^pm2,p)
if mpz_cmp_si(t,1)!=0 then
return "No solution exists"
end if
mpz q = mpz_init_set(pm1)
integer ss = 0
while mpz_even(q) do
ss += 1
mpz_fdiv_q_2exp(q,q,1) -- q/=2
end while
if ss=1 then
mpz_add_ui(t,p,1)
mpz_fdiv_q_2exp(t,t,2)
mpz_powm(r,n,t,p) -- r = mod(n^((p+1)/4),p)
else
mpz z = mpz_init(2)
while true do
mpz_powm(t,z,pm2,p) -- t = mod(z^pm2,p)
if mpz_cmp(t,pm1)=0 then exit end if
mpz_add_ui(z,z,1) -- z+= 1
end while
mpz {b,c,zz} = mpz_inits(3)
mpz_powm(c,z,q,p) -- c = mod(z^q,p)
mpz_add_ui(t,q,1)
mpz_fdiv_q_2exp(t,t,1)
mpz_powm(r,n,t,p) -- r = mod(n^((q+1)/2),p)
mpz_powm(t,n,q,p) -- t = mod(n^q,p)
integer m = ss
while mpz_cmp_si(t,1) do -- t!=1
integer i = 0
mpz_set(zz,t)
while mpz_cmp_si(zz,1)!=0 and i<m-1 do
mpz_powm_ui(zz,zz,2,p) -- zz = mod(zz^2,p)
i += 1
end while
mpz_set(b,c)
integer e = m-i-1
while e>0 do
mpz_powm_ui(b,b,2,p) -- b = mod(b^2,p)
e -= 1
end while
mpz_mul(r,r,b)
mpz_mod(r,r,p) -- r = mod(r*b,p)
mpz_powm_ui(c,b,2,p) -- c = mod(b^2,p)
mpz_mul(t,t,c)
mpz_mod(t,t,p) -- t = mod(t*c,p)
m = i
end while
end if
mpz_sub(p,p,r)
return mpz_get_str(r)&" and "&mpz_get_str(p)
end function
constant tests = {{"10","13"},
{"56","101"},
{"1030","10009"},
{"1032","10009"},
{"44402","100049"},
{"665820697","1000000009"},
{"881398088036","1000000000039"},
{"41660815127637347468140745042827704103445750172002",
sprintf("1%s577",repeat('0',47))}} -- 10^50+577
for i=1 to length(tests) do
string {p1,p2} = tests[i]
printf(1,"For n = %s and p = %s, %s\n",{p1,p2,ts(p1,p2)})
end for