

 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


// version 1.1.3
 
import java.math.BigInteger
 
data class Solution(val root1: BigInteger, val root2: BigInteger, val exists: Boolean)
 
val bigZero = BigInteger.ZERO
val bigOne  = BigInteger.ONE
val bigTwo  = BigInteger.valueOf(2L)
val bigFour = BigInteger.valueOf(4L)
val bigTen  = BigInteger.TEN
 
fun ts(n: Long, p: Long) = ts(BigInteger.valueOf(n), BigInteger.valueOf(p))
 
fun ts(n: BigInteger, p: BigInteger): Solution {
 
    fun powModP(a: BigInteger, e: BigInteger) = a.modPow(e, p)
 
    fun ls(a: BigInteger) = powModP(a, (p - bigOne) / bigTwo)
 
    if (ls(n) != bigOne) return Solution(bigZero, bigZero, false)
    var q = p - bigOne
    var ss = bigZero
    while (q.and(bigOne) == bigZero) {
        ss = ss + bigOne
        q = q.shiftRight(1)
    }
 
    if (ss == bigOne) {
        val r1 = powModP(n, (p + bigOne) / bigFour)
        return Solution(r1, p - r1, true)
    }
 
    var z = bigTwo
    while (ls(z) != p - bigOne) z = z + bigOne
    var c = powModP(z, q)
    var r = powModP(n, (q + bigOne) / bigTwo)
    var t = powModP(n, q)
    var m = ss
 
    while (true) {
        if (t == bigOne) return Solution(r, p - r, true)
        var i = bigZero
        var zz = t
        while (zz != bigOne && i < m - bigOne) {
            zz  = zz * zz % p
            i = i + bigOne
        }
        var b = c
        var e = m - i - bigOne
        while (e > bigZero) {
            b = b * b % p
            e = e - bigOne
        }
        r = r * b % p
        c = b * b % p
        t = t * c % p
        m = i
    }
}
 
fun main(args: Array<String>) {
    val pairs = listOf<Pair<Long, Long>>(
        10L to 13L, 
        56L to 101L, 
        1030L to 10009L,
        1032L to 10009L,
        44402L to 100049L,
        665820697L to 1000000009L,
        881398088036L to 1000000000039L
    )
 
    for (pair in pairs) {
        val (n, p) = pair
        val (root1, root2, exists) = ts(n, p)
        println("n = $n")
        println("p = $p")
        if (exists) {
            println("root1 = $root1")
            println("root2 = $root2")
        }
        else println("No solution exists")
        println()
    }
 
    val bn = BigInteger("41660815127637347468140745042827704103445750172002")
    val bp = bigTen.pow(50) + BigInteger.valueOf(577L)
    val (broot1, broot2, bexists) = ts(bn, bp)
    println("n = $bn")
    println("p = $bp")
    if (bexists) {
        println("root1 = $broot1")
        println("root2 = $broot2")
    }
    else println("No solution exists")    
}