Cipolla's algorithm: Difference between revisions

Content added Content deleted

Inline

@@ Line 1,442: / Line 1,442: @@
  {"82563118828090362261378993957450213573687113690751","17436881171909637738621006042549786426312886309400","true"}}
 </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}}==