func tonelli(n, p) {
    legendre(n, p) == 1 || die "not a square (mod p)"
    var q = p-1
    var s = valuation(q, 2)
    s == 1 ? return(powmod(n, (p + 1) >> 2, p)) : (q >>= s)
    var c = powmod(2 ..^ p -> first {|z| legendre(z, p) == -1}, q, p)
    var r = powmod(n, (q + 1) >> 1, p)
    var t = powmod(n, q, p)
    var m = s
    var t2 = 0
    while (!p.divides(t - 1)) {
        t2 = ((t * t) % p)
        var b
        for i in (1 ..^ m) {
            if (p.divides(t2 - 1)) {
                b = powmod(c, 1 << (m - i - 1), p)
                m = i
                break
            }
            t2 = ((t2 * t2) % p)
        }

        r = ((r * b) % p)
        c = ((b * b) % p)
        t = ((t * c) % p)
    }
    return r
}

var tests = [
    [10, 13], [56, 101], [1030, 10009], [44402, 100049],
    [665820697, 1000000009], [881398088036, 1000000000039],
    [41660815127637347468140745042827704103445750172002, 10**50 + 577],
]

for n,p in tests {
    var r = tonelli(n, p)
    assert((r*r - n) % p == 0)
    say "Roots of #{n} are (#{r}, #{p-r}) mod #{p}"
}