cipolla_algorithm.sf

#!/usr/bin/ruby

# Cipolla's algorithm, for solving a congruence of the form:
#   x^2 == n (mod p)
# where p is an odd prime.

# See also:
#   https://en.wikipedia.org/wiki/Cipolla's_algorithm
#   https://rosettacode.org/wiki/Cipolla%27s_algorithm

func cipolla(n, p) {

    legendre(n, p) == 1 || return nil

    var (a = 0, ω2 = 0)
    loop {
        ω2 = ((a*a - n) % p)
        if (kronecker(ω2, p) == -1) {
            break
        }
        ++a
    }

    struct point { x, y }

    func mul(a, b) {
        point((a.x*b.x + a.y*b.y*ω2) % p, (a.x*b.y + b.x*a.y) % p)
    }

    var r = point(1, 0)
    var s = point(a, 1)

    for (var n = ((p+1) >> 1); n > 0; n >>= 1) {
        r = mul(r, s) if n.is_odd
        s = mul(s, s)
    }

    r.y == 0 ? r.x : nil
}

var tests = [
    [10, 13],
    [56, 101],
    [8218, 10007],
    [8219, 10007],
    [331575, 1000003],
    [665165880, 1000000007],
    [881398088036 1000000000039],
    [34035243914635549601583369544560650254325084643201, 10**50 + 151],
]

for n,p in tests {
    var r = cipolla(n, p)
    if (defined(r)) {
        say "Roots of #{n} are (#{r} #{p-r}) mod #{p}"
    } else {
        say "No solution for (#{n}, #{p})"
    }
}