Cipolla's algorithm: Difference between revisions
Content added Content deleted
(→{{header|Phix}}: added gmp version) |
|||
Line 1,442: | Line 1,442: | ||
{"82563118828090362261378993957450213573687113690751","17436881171909637738621006042549786426312886309400","true"}} |
{"82563118828090362261378993957450213573687113690751","17436881171909637738621006042549786426312886309400","true"}} |
||
</pre> |
</pre> |
||
=== gmp === |
|||
<lang Phix>include mpfr.e |
|||
procedure legendre(mpz r, a, p) |
|||
-- Legendre symbol, returns 1, 0 or p - 1 (in r) |
|||
mpz_sub_ui(r,p,1) |
|||
{} = mpz_fdiv_q_ui(r, r, 2) |
|||
mpz_powm(r,a,r,p) |
|||
end procedure |
|||
procedure mul_point(sequence a, b, mpz omega2, p) |
|||
-- (modifies a) |
|||
mpz {xa,ya} = a, |
|||
{xb,yb} = b, |
|||
xaxb = mpz_init(), |
|||
yayb = mpz_init(), |
|||
xayb = mpz_init(), |
|||
xbya = mpz_init() |
|||
mpz_mul(xaxb,xa,xb) |
|||
mpz_mul(yayb,ya,yb) |
|||
mpz_mul(xayb,xa,yb) |
|||
mpz_mul(xbya,xb,ya) |
|||
mpz_mul(yayb,yayb,omega2) |
|||
mpz_add(xaxb,xaxb,yayb) |
|||
mpz_mod(xa,xaxb,p) -- xa := mod(xaxb+yayb*omega2,p) |
|||
mpz_add(xayb,xayb,xbya) |
|||
mpz_mod(ya,xayb,p) -- ya := mod(xayb+xbya,p) |
|||
{xaxb,yayb,xayb,xbya} = mpz_clear({xaxb,yayb,xayb,xbya}) |
|||
end procedure |
|||
function cipolla(object no, po) |
|||
mpz n = mpz_init(no), |
|||
p = mpz_init(po), |
|||
t = mpz_init() |
|||
-- Step 0, validate arguments |
|||
legendre(t,n,p) |
|||
if mpz_cmp_si(t,1)!=0 then return {"0","0","false"} end if |
|||
-- Step 1, find a, omega2 |
|||
integer a = 0 |
|||
mpz omega2 = mpz_init(), |
|||
pm1 = mpz_init() |
|||
mpz_sub_ui(pm1,p,1) |
|||
while true do |
|||
mpz_sub(t,p,n) |
|||
mpz_add_ui(t,t,a*a) |
|||
mpz_mod(omega2,t,p) |
|||
legendre(t,omega2,p) |
|||
if mpz_cmp(t,pm1)=0 then exit end if |
|||
a += 1 |
|||
end while |
|||
-- Step 2, compute power |
|||
sequence r = {mpz_init(1),mpz_init(0)}, |
|||
s = {mpz_init(a),mpz_init(1)} |
|||
mpz nn = mpz_init() |
|||
mpz_add_ui(nn,p,1) |
|||
{} = mpz_fdiv_q_ui(nn, nn, 2) |
|||
mpz_mod(nn,nn,p) |
|||
while mpz_cmp_si(nn,0)>0 do |
|||
if mpz_fdiv_ui(nn,2)=1 then |
|||
mul_point(r,s,omega2,p) |
|||
end if |
|||
mul_point(s,s,omega2,p) |
|||
{} = mpz_fdiv_q_ui(nn, nn, 2) |
|||
end while |
|||
-- Step 3, check x in Fp |
|||
if mpz_cmp_si(r[2],0)!=0 then return {"0","0","false"} end if |
|||
-- Step 5, check x * x = n |
|||
mpz_powm_ui(t,r[1],2,p) |
|||
if mpz_cmp(t,n)!=0 then return {"0","0","false"} end if |
|||
-- Step 4, solutions |
|||
mpz_sub(p,p,r[1]) |
|||
return {mpz_get_str(r[1]), mpz_get_str(p), "true"} |
|||
end function |
|||
constant tests = {{10,13}, |
|||
{56,101}, |
|||
{8218,10007}, |
|||
{8219,10007}, |
|||
{331575,1000003}, |
|||
{665165880,1000000007}, |
|||
{"881398088036","1000000000039"}, |
|||
{"34035243914635549601583369544560650254325084643201", |
|||
"100000000000000000000000000000000000000000000000151"}} |
|||
for i=1 to length(tests) do |
|||
object {n,p} = tests[i] |
|||
?{n,p,cipolla(n,p)} |
|||
end for</lang> |
|||
same output |
|||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |