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