

 This page uses content from Wikipedia. The original article was at Tonelli-Shanks algorithm. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



In computational number theory, the Tonelli–Shanks algorithm is a technique for solving for x in a congruence of the form:


 x2 ≡ n (mod p)

where n is an integer which is a quadratic residue (mod p), p is an odd prime, and x,n ∈ Fp where Fp = {0, 1, ..., p - 1}.

It is used in cryptography techniques.



To apply the algorithm, we need the Legendre symbol:

The Legendre symbol (a | p) denotes the value of a(p-1)/2 (mod p).

 (a | p) ≡ 1    if a is a square (mod p)
 (a | p) ≡ -1    if a is not a square (mod p)
 (a | p) ≡ 0    if a ≡ 0 (mod p) 


Algorithm pseudo-code


All ≡ are taken to mean (mod p) unless stated otherwise.

 Input: p an odd prime, and an integer n .
 Step 0: Check that n is indeed a square: (n | p) must be ≡ 1 .
 Step 1: By factoring out powers of 2 from p - 1, find q and s such that p - 1 = q2s with q odd .
 If p ≡ 3 (mod 4) (i.e. s = 1), output the two solutions r ≡ ± n(p+1)/4 .
 Step 2: Select a non-square z such that (z | p) ≡ -1 and set c ≡ zq .
 Step 3: Set r ≡ n(q+1)/2, t ≡ nq, m = s .
 Step 4: Loop the following:
 If t ≡ 1, output r and p - r .
 Otherwise find, by repeated squaring, the lowest i, 0 < i < m , such that t2i ≡ 1 .
 Let b ≡ c2(m - i - 1), and set r ≡ rb, t ≡ tb2, c ≡ b2 and m = i .




Task

Implement the above algorithm.

Find solutions (if any) for

 n = 10 p = 13
 n = 56 p = 101
 n = 1030 p = 10009
 n = 1032 p = 10009
 n = 44402 p = 100049 
Extra credit

 n = 665820697 p = 1000000009 
 n = 881398088036 p = 1000000000039 
 n = 41660815127637347468140745042827704103445750172002 p = 10^50 + 577 


See also

 Modular exponentiation
 Cipolla's algorithm


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
