Алгоритм Тонелли — Шенкса

import java.math.BigInteger;

import static java.math.BigInteger.ONE;
import static java.math.BigInteger.ZERO;
import static java.math.BigInteger.valueOf;

public class TonelliShanksAlgorithm {

    static final BigInteger TWO = valueOf(2);
    static final BigInteger THREE = valueOf(3);
    static final BigInteger FOUR = valueOf(4);

    public static void main(String[] args) {
        printRoots(103, 1789);
        printRoots("24296403966", "32416188697");
        printRoots("84235199442121369408907697", "241573142393627673576957439049");
        printRoots("6706901388311723110433121899414477003899", "45994811347886846310221728895223034301839");
    }

    static void printRoots(String n, String p) {
        printRoots(new BigInteger(n), new BigInteger(p));
    }

    static void printRoots(int n, int p) {
        printRoots(valueOf(n), valueOf(p));
    }

    static void printRoots(BigInteger n, BigInteger p) {
        BigInteger a = getRoot(n, p);
        a = a.min(p.subtract(a));
        System.out.println(a + ", " + p.subtract(a));
    }

    static BigInteger getRoot(BigInteger n, BigInteger p) {
        BigInteger s = ZERO;
        BigInteger q = p.subtract(ONE);

        while (true) {
            BigInteger[] quotientAndRem = q.divideAndRemainder(TWO);

            if (quotientAndRem[1].equals(ZERO)) {
                q = quotientAndRem[0];
                s = s.add(ONE);
            } else
                break;
        }

        if (s.equals(ONE))
            return n.modPow(p.add(ONE).divide(FOUR), p);

        BigInteger z = TWO;
        for (; ; z = z.add(ONE))
            if (jacobiSymbol(z, p) == -1)
                break;

        BigInteger c = z.modPow(q, p);
        BigInteger r = n.modPow(q.add(ONE).divide(TWO), p);
        BigInteger t = n.modPow(q, p);
        BigInteger m = s;

        while (true) {
            if (t.equals(ONE))
                return r;

            int i = 1;
            BigInteger d = t;
            for (; i < m.intValue(); i++) {
                d = d.multiply(d).mod(p);
                if (d.equals(ONE))
                    break;
            }

            BigInteger degree = TWO.modPow(m.subtract(valueOf(i + 1)), p.subtract(ONE));
            BigInteger b = c.modPow(degree, p);
            r = r.multiply(b).mod(p);
            c = b.multiply(b).mod(p);
            t = t.multiply(c).mod(p);
            m = valueOf(i);
        }
    }

    static int jacobiSymbol(BigInteger a, BigInteger n) {
        if (a.equals(ZERO))
            return (n.equals(ONE)) ? 1 : 0;

        if (a.equals(TWO))
            switch (n.mod(valueOf(8)).intValue()) {
                case 1:
                case 7:
                    return 1;
                case 3:
                case 5:
                    return -1;
            }

        if (a.compareTo(n) >= 0)
            return jacobiSymbol(a.mod(n), n);

        if (a.mod(TWO).equals(ZERO))
            return jacobiSymbol(TWO, n) * jacobiSymbol(a.divide(TWO), n);

        return (a.mod(FOUR).equals(THREE) && n.mod(FOUR).equals(THREE))
                ? -jacobiSymbol(n, a) : jacobiSymbol(n, a);
    }
}