cipolla_algorithm_simple.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.

# This implementation uses quadratic integers.

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

func cipolla(n, p) {

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

    var (a, ω) = (
        0..Inf -> lazy.map {|a|
            [a, submod(a*a, n, p)]
        }.first_by {|t|
            kronecker(t[1], p) == -1
        }...
    )

    var r = lift(Mod(Quadratic(a, 1, ω), p)**((p+1)>>1))

    r.b == 0 ? r.a : 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})"
    }
}