123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960 |
- #!/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})"
- }
- }
|