

 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


import java.math.BigInteger;
import java.util.List;
import java.util.Map;
import java.util.function.BiFunction;
import java.util.function.Function;
 
public class TonelliShanks {
    private static final BigInteger ZERO = BigInteger.ZERO;
    private static final BigInteger ONE = BigInteger.ONE;
    private static final BigInteger TEN = BigInteger.TEN;
    private static final BigInteger TWO = BigInteger.valueOf(2);
    private static final BigInteger FOUR = BigInteger.valueOf(4);
 
    private static class Solution {
        private BigInteger root1;
        private BigInteger root2;
        private boolean exists;
 
        Solution(BigInteger root1, BigInteger root2, boolean exists) {
            this.root1 = root1;
            this.root2 = root2;
            this.exists = exists;
        }
    }
 
    private static Solution ts(Long n, Long p) {
        return ts(BigInteger.valueOf(n), BigInteger.valueOf(p));
    }
 
    private static Solution ts(BigInteger n, BigInteger p) {
        BiFunction<BigInteger, BigInteger, BigInteger> powModP = (BigInteger a, BigInteger e) -> a.modPow(e, p);
        Function<BigInteger, BigInteger> ls = (BigInteger a) -> powModP.apply(a, p.subtract(ONE).divide(TWO));
 
        if (!ls.apply(n).equals(ONE)) return new Solution(ZERO, ZERO, false);
 
        BigInteger q = p.subtract(ONE);
        BigInteger ss = ZERO;
        while (q.and(ONE).equals(ZERO)) {
            ss = ss.add(ONE);
            q = q.shiftRight(1);
        }
 
        if (ss.equals(ONE)) {
            BigInteger r1 = powModP.apply(n, p.add(ONE).divide(FOUR));
            return new Solution(r1, p.subtract(r1), true);
        }
 
        BigInteger z = TWO;
        while (!ls.apply(z).equals(p.subtract(ONE))) z = z.add(ONE);
        BigInteger c = powModP.apply(z, q);
        BigInteger r = powModP.apply(n, q.add(ONE).divide(TWO));
        BigInteger t = powModP.apply(n, q);
        BigInteger m = ss;
 
        while (true) {
            if (t.equals(ONE)) return new Solution(r, p.subtract(r), true);
            BigInteger i = ZERO;
            BigInteger zz = t;
            while (!zz.equals(BigInteger.ONE) && i.compareTo(m.subtract(ONE)) < 0) {
                zz = zz.multiply(zz).mod(p);
                i = i.add(ONE);
            }
            BigInteger b = c;
            BigInteger e = m.subtract(i).subtract(ONE);
            while (e.compareTo(ZERO) > 0) {
                b = b.multiply(b).mod(p);
                e = e.subtract(ONE);
            }
            r = r.multiply(b).mod(p);
            c = b.multiply(b).mod(p);
            t = t.multiply(c).mod(p);
            m = i;
        }
    }
 
    public static void main(String[] args) {
        List<Map.Entry<Long, Long>> pairs = List.of(
            Map.entry(10L, 13L),
            Map.entry(56L, 101L),
            Map.entry(1030L, 10009L),
            Map.entry(1032L, 10009L),
            Map.entry(44402L, 100049L),
            Map.entry(665820697L, 1000000009L),
            Map.entry(881398088036L, 1000000000039L)
        );
 
        for (Map.Entry<Long, Long> pair : pairs) {
            Solution sol = ts(pair.getKey(), pair.getValue());
            System.out.printf("n = %s\n", pair.getKey());
            System.out.printf("p = %s\n", pair.getValue());
            if (sol.exists) {
                System.out.printf("root1 = %s\n", sol.root1);
                System.out.printf("root2 = %s\n", sol.root2);
            } else {
                System.out.println("No solution exists");
            }
            System.out.println();
        }
 
        BigInteger bn = new BigInteger("41660815127637347468140745042827704103445750172002");
        BigInteger bp = TEN.pow(50).add(BigInteger.valueOf(577));
        Solution sol = ts(bn, bp);
        System.out.printf("n = %s\n", bn);
        System.out.printf("p = %s\n", bp);
        if (sol.exists) {
            System.out.printf("root1 = %s\n", sol.root1);
            System.out.printf("root2 = %s\n", sol.root2);
        } else {
            System.out.println("No solution exists");
        }
    }
}