1234567891011121314151617181920212223242526272829303132333435363738394041424344454647 |
- #!/usr/bin/ruby
- #
- ## https://rosettacode.org/wiki/Tonelli-Shanks_algorithm
- #
- 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}"
- }
|