Tonelli-Shanks algorithm: Difference between revisions
Added solution for D |
Added Java |
||
Line 910: | Line 910: | ||
41660815127637347468140745042827704103445750172002x tosh (10^50x)+577 |
41660815127637347468140745042827704103445750172002x tosh (10^50x)+577 |
||
32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069</lang> |
32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069</lang> |
||
=={{header|Java}}== |
|||
{{trans|Kotlin}} |
|||
{{works with|Java|9}} |
|||
<lang Java>import java.math.BigInteger; |
|||
import java.util.List; |
|||
import java.util.Map; |
|||
import java.util.function.BiFunction; |
|||
import java.util.function.Function; |
|||
public class TonelliShanks { |
|||
private static final BigInteger ZERO = BigInteger.ZERO; |
|||
private static final BigInteger ONE = BigInteger.ONE; |
|||
private static final BigInteger TEN = BigInteger.TEN; |
|||
private static final BigInteger TWO = BigInteger.valueOf(2); |
|||
private static final BigInteger FOUR = BigInteger.valueOf(4); |
|||
private static class Solution { |
|||
private BigInteger root1; |
|||
private BigInteger root2; |
|||
private boolean exists; |
|||
Solution(BigInteger root1, BigInteger root2, boolean exists) { |
|||
this.root1 = root1; |
|||
this.root2 = root2; |
|||
this.exists = exists; |
|||
} |
|||
} |
|||
private static Solution ts(Long n, Long p) { |
|||
return ts(BigInteger.valueOf(n), BigInteger.valueOf(p)); |
|||
} |
|||
private static Solution ts(BigInteger n, BigInteger p) { |
|||
BiFunction<BigInteger, BigInteger, BigInteger> powModP = (BigInteger a, BigInteger e) -> a.modPow(e, p); |
|||
Function<BigInteger, BigInteger> ls = (BigInteger a) -> powModP.apply(a, p.subtract(ONE).divide(TWO)); |
|||
if (!ls.apply(n).equals(ONE)) return new Solution(ZERO, ZERO, false); |
|||
BigInteger q = p.subtract(ONE); |
|||
BigInteger ss = ZERO; |
|||
while (q.and(ONE).equals(ZERO)) { |
|||
ss = ss.add(ONE); |
|||
q = q.shiftRight(1); |
|||
} |
|||
if (ss.equals(ONE)) { |
|||
BigInteger r1 = powModP.apply(n, p.add(ONE).divide(FOUR)); |
|||
return new Solution(r1, p.subtract(r1), true); |
|||
} |
|||
BigInteger z = TWO; |
|||
while (!ls.apply(z).equals(p.subtract(ONE))) z = z.add(ONE); |
|||
BigInteger c = powModP.apply(z, q); |
|||
BigInteger r = powModP.apply(n, q.add(ONE).divide(TWO)); |
|||
BigInteger t = powModP.apply(n, q); |
|||
BigInteger m = ss; |
|||
while (true) { |
|||
if (t.equals(ONE)) return new Solution(r, p.subtract(r), true); |
|||
BigInteger i = ZERO; |
|||
BigInteger zz = t; |
|||
while (!zz.equals(BigInteger.ONE) && i.compareTo(m.subtract(ONE)) < 0) { |
|||
zz = zz.multiply(zz).mod(p); |
|||
i = i.add(ONE); |
|||
} |
|||
BigInteger b = c; |
|||
BigInteger e = m.subtract(i).subtract(ONE); |
|||
while (e.compareTo(ZERO) > 0) { |
|||
b = b.multiply(b).mod(p); |
|||
e = e.subtract(ONE); |
|||
} |
|||
r = r.multiply(b).mod(p); |
|||
c = b.multiply(b).mod(p); |
|||
t = t.multiply(c).mod(p); |
|||
m = i; |
|||
} |
|||
} |
|||
public static void main(String[] args) { |
|||
List<Map.Entry<Long, Long>> pairs = List.of( |
|||
Map.entry(10L, 13L), |
|||
Map.entry(56L, 101L), |
|||
Map.entry(1030L, 10009L), |
|||
Map.entry(1032L, 10009L), |
|||
Map.entry(44402L, 100049L), |
|||
Map.entry(665820697L, 1000000009L), |
|||
Map.entry(881398088036L, 1000000000039L) |
|||
); |
|||
for (Map.Entry<Long, Long> pair : pairs) { |
|||
Solution sol = ts(pair.getKey(), pair.getValue()); |
|||
System.out.printf("n = %s\n", pair.getKey()); |
|||
System.out.printf("p = %s\n", pair.getValue()); |
|||
if (sol.exists) { |
|||
System.out.printf("root1 = %s\n", sol.root1); |
|||
System.out.printf("root2 = %s\n", sol.root2); |
|||
} else { |
|||
System.out.println("No solution exists"); |
|||
} |
|||
System.out.println(); |
|||
} |
|||
BigInteger bn = new BigInteger("41660815127637347468140745042827704103445750172002"); |
|||
BigInteger bp = TEN.pow(50).add(BigInteger.valueOf(577)); |
|||
Solution sol = ts(bn, bp); |
|||
System.out.printf("n = %s\n", bn); |
|||
System.out.printf("p = %s\n", bp); |
|||
if (sol.exists) { |
|||
System.out.printf("root1 = %s\n", sol.root1); |
|||
System.out.printf("root2 = %s\n", sol.root2); |
|||
} else { |
|||
System.out.println("No solution exists"); |
|||
} |
|||
} |
|||
}</lang> |
|||
{{out}} |
|||
<pre>n = 10 |
|||
p = 13 |
|||
root1 = 7 |
|||
root2 = 6 |
|||
n = 56 |
|||
p = 101 |
|||
root1 = 37 |
|||
root2 = 64 |
|||
n = 1030 |
|||
p = 10009 |
|||
root1 = 1632 |
|||
root2 = 8377 |
|||
n = 1032 |
|||
p = 10009 |
|||
No solution exists |
|||
n = 44402 |
|||
p = 100049 |
|||
root1 = 30468 |
|||
root2 = 69581 |
|||
n = 665820697 |
|||
p = 1000000009 |
|||
root1 = 378633312 |
|||
root2 = 621366697 |
|||
n = 881398088036 |
|||
p = 1000000000039 |
|||
root1 = 791399408049 |
|||
root2 = 208600591990 |
|||
n = 41660815127637347468140745042827704103445750172002 |
|||
p = 100000000000000000000000000000000000000000000000577 |
|||
root1 = 32102985369940620849741983987300038903725266634508 |
|||
root2 = 67897014630059379150258016012699961096274733366069</pre> |
|||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
Revision as of 01:54, 17 February 2018
Tonelli–Shanks algorithm
In computational number theory, the Tonelli–Shanks algorithm is a technique for solving an equation of the form:
- x2 ≡ n (mod p)
─── where n is an integer which is a quadratic residue (mod p), p is an odd prime, and
- x,n Є Fp = {0, 1, ... p-1}
It is used in cryptography techniques.
To apply the algorithm, we need the Legendre symbol.
Legendre symbol
- The Legendre symbol ( a | p) denotes the value of a ^ ((p-1)/2) (mod p)
- (a | p) ≡ 1 if a is a square (mod p)
- (a | p) ≡ -1 if a is not a square (mod p)
- (a | p) ≡ 0 if a ≡ 0
- Algorithm pseudo-code
- (copied from Wikipedia):
All ≡ are taken to mean (mod p) unless stated otherwise.
- Input : p an odd prime, and an integer n .
- Step 0. Check that n is indeed a square : (n | p) must be ≡ 1
- Step 1. [Factors out powers of 2 from p-1] Define q -odd- and s such as p-1 = q * 2^s
- if s = 1 , i.e p ≡ 3 (mod 4) , output the two solutions r ≡ +/- n^((p+1)/4) .
- Step 2. Select a non-square z such as (z | p) = -1 , and set c ≡ z^q .
- Step 3. Set r ≡ n ^((q+1)/2) , t ≡ n^q, m = s .
- Step 4. Loop.
- if t ≡ 1 output r, p-r .
- Otherwise find, by repeated squaring, the lowest i , 0 < i< m , such as t^(2^i) ≡ 1
- Let b ≡ c^(2^(m-i-1)), and set r ≡ r*b, t ≡ t*b^2 , c ≡ b^2 and m = i.
- Numerical Example
- n=10, p= 13. See Wikipedia
- Task
Implement the above.
Find solutions (if any) for
- n = 10 p = 13
- n = 56 p = 101
- n = 1030 p = 10009
- n = 1032 p = 10009
- n = 44402 p = 100049
- Extra credit
- n = 665820697 p = 1000000009
- n = 881398088036 p = 1000000000039
- n = 41660815127637347468140745042827704103445750172002 p = 10^50 + 577
- See also
D
<lang D>import std.bigint; import std.stdio; import std.typecons;
alias Pair = Tuple!(long, "n", long, "p");
enum BIGZERO = BigInt("0"); enum BIGONE = BigInt("1"); enum BIGTWO = BigInt("2"); enum BIGTEN = BigInt("10");
struct Solution {
BigInt root1, root2; bool exists;
}
/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method BigInt modPow(BigInt b, BigInt e, BigInt n) {
if (n == 1) return BIGZERO; BigInt result = 1; b = b % n; while (e > 0) { if (e % 2 == 1) { result = (result * b) % n; } e >>= 1; b = (b*b) % n; } return result;
}
Solution ts(long n, long p) {
return ts(BigInt(n), BigInt(p));
}
Solution ts(BigInt n, BigInt p) {
auto powMod(BigInt a, BigInt e) { return a.modPow(e, p); }
auto ls(BigInt a) { return powMod(a, (p-1)/2); }
if (ls(n) != 1) return Solution(BIGZERO, BIGZERO, false); auto q = p - 1; auto ss = BIGZERO; while ((q & 1) == 0) { ss = ss + 1; q = q >> 1; }
if (ss == BIGONE) { auto r1 = powMod(n, (p + 1) / 4); return Solution(r1, p - r1, true); }
auto z = BIGTWO; while (ls(z) != p - 1) z = z + 1; auto c = powMod(z, q); auto r = powMod(n, (q + 1) / 2); auto t = powMod(n, q); auto m = ss;
while (true) { if (t == 1) return Solution(r, p - r, true); auto i = BIGZERO; auto zz = t; while (zz != 1 && i < m - 1) { zz = zz * zz % p; i = i + 1; } auto b = c; auto e = m - i - 1; while (e > 0) { b = b * b % p; e = e - 1; } r = r * b % p; c = b * b % p; t = t * c % p; m = i; }
}
void main() {
auto pairs = [ Pair( 10L, 13L), Pair( 56L, 101L), Pair( 1_030L, 10_009L), Pair( 1_032L, 10_009L), Pair( 44_402L, 100_049L), Pair( 665_820_697L, 1_000_000_009L), Pair(881_398_088_036L, 1_000_000_000_039L), ];
foreach (pair; pairs) { auto sol = ts(pair.n, pair.p);
writeln("n = ", pair.n); writeln("p = ", pair.p); if (sol.exists) { writeln("root1 = ", sol.root1); writeln("root2 = ", sol.root2); } else writeln("No solution exists"); writeln(); }
auto bn = BigInt("41660815127637347468140745042827704103445750172002"); auto bp = BIGTEN ^^ 50 + 577L; auto sol = ts(bn, bp); writeln("n = ", bn); writeln("p = ", bp); if (sol.exists) { writeln("root1 = ", sol.root1); writeln("root2 = ", sol.root2); } else writeln("No solution exists");
}</lang>
- Output:
n = 10 p = 13 root1 = 7 root2 = 6 n = 56 p = 101 root1 = 37 root2 = 64 n = 1030 p = 10009 root1 = 1632 root2 = 8377 n = 1032 p = 10009 No solution exists n = 44402 p = 100049 root1 = 30468 root2 = 69581 n = 665820697 p = 1000000009 root1 = 378633312 root2 = 621366697 n = 881398088036 p = 1000000000039 root1 = 791399408049 root2 = 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 root1 = 32102985369940620849741983987300038903725266634508 root2 = 67897014630059379150258016012699961096274733366069
EchoLisp
<lang scheme> (require 'bigint)
- test equality mod p
(define-syntax-rule (mod= a b p) (zero? (% (- a b) p)))
- assign mod p
(define-syntax-rule (mod:≡ s v p) (set! s (% v p)))
(define (Legendre a p) (powmod a (/ (1- p) 2) p))
(define (Tonelli n p)
(unless (= 1 (Legendre n p)) (error "not a square (mod p)" (list n p))) (define q (1- p)) (define s 0)
(while (even? q) (/= q 2) (++ s)) (if (= s 1) (powmod n (/ (1+ p) 4) p) (begin (define z (for ((z (in-range 2 p))) #:break (= (1- p) (Legendre z p)) => z ))
(define c (powmod z q p)) (define r (powmod n (/ (1+ q) 2) p)) (define t (powmod n q p)) (define m s) (define t2 0) (while #t #:break (mod= 1 t p) => r (mod:≡ t2 (* t t) p) (define i (for ((i (in-range 1 m))) #:break (mod= t2 1 p) => i (mod:≡ t2 (* t2 t2) p))) (define b (powmod c (expt 2 (- m i 1)) p)) (mod:≡ r (* r b) p) (mod:≡ c (* b b) p) (mod:≡ t (* t c) p) (set! m i))))) </lang>
- Output:
(define ttest `((10 13) (56 101) (1030 10009) (44402 100049) (665820697 1000000009) (881398088036 1000000000039) (41660815127637347468140745042827704103445750172002 ,(+ 1e50 577)))) (define (task ttest) (for ((test ttest)) (define n (first test)) (define p (second test)) (define r (Tonelli n p)) (assert (mod= (* r r) n p)) (printf "n = %d p = %d" n p) (printf "\t roots : %d %d" r (- p r)))) (task ttest) n = 10 p = 13 roots : 7 6 n = 56 p = 101 roots : 37 64 n = 1030 p = 10009 roots : 1632 8377 n = 44402 p = 100049 roots : 30468 69581 n = 665820697 p = 1000000009 roots : 378633312 621366697 n = 881398088036 p = 1000000000039 roots : 791399408049 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 roots : 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069 (Tonelli 1032 10009) ❌ error: not a square (mod p) (1032 10009)
FreeBASIC
LongInt version
<lang FreeBASIC>' version 11-04-2017 ' compile with: fbc -s console ' maximum for p is 17 digits to be on the save side
' TRUE/FALSE are built-in constants since FreeBASIC 1.04 ' But we have to define them for older versions.
- Ifndef TRUE
#Define FALSE 0 #Define TRUE Not FALSE
- EndIf
Function mul_mod(a As ULongInt, b As ULongInt, modulus As ULongInt) As ULongInt
' returns a * b mod modulus Dim As ULongInt x, y = a Mod modulus
While b > 0 If (b And 1) = 1 Then x = (x + y) Mod modulus End If y = (y Shl 1) Mod modulus b = b Shr 1 Wend
Return x
End Function
Function pow_mod(b As ULongInt, power As ULongInt, modulus As ULongInt) As ULongInt
' returns b ^ power mod modulus Dim As ULongInt x = 1
While power > 0 If (power And 1) = 1 Then ' x = (x * b) Mod modulus x = mul_mod(x, b, modulus) End If ' b = (b * b) Mod modulus b = mul_mod(b, b, modulus) power = power Shr 1 Wend
Return x
End Function
Function Isprime(n As ULongInt, k As Long) As Long
' miller-rabin prime test If n > 9223372036854775808ull Then ' limit 2^63, pow_mod/mul_mod can't handle bigger numbers Print "number is to big, program will end" Sleep End End If
' 2 is a prime, if n is smaller then 2 or n is even then n = composite If n = 2 Then Return TRUE If (n < 2) OrElse ((n And 1) = 0) Then Return FALSE
Dim As ULongInt a, x, n_one = n - 1, d = n_one Dim As UInteger s
While (d And 1) = 0 d = d Shr 1 s = s + 1 Wend
While k > 0 k = k - 1 a = Int(Rnd * (n -2)) +2 ' 2 <= a < n x = pow_mod(a, d, n) If (x = 1) Or (x = n_one) Then Continue While For r As Integer = 1 To s -1 x = pow_mod(x, 2, n) If x = 1 Then Return FALSE If x = n_one Then Continue While Next If x <> n_one Then Return FALSE Wend Return TRUE
End Function
Function legendre_symbol (a As LongInt, p As LongInt) As LongInt
Dim As LongInt x = pow_mod(a, ((p -1) \ 2), p) If p -1 = x Then Return x - p Else Return x End If
End Function
' ------=< MAIN >=------
Dim As LongInt b, c, i, k, m, n, p, q, r, s, t, z
For k = 1 To 7
Read n, p Print "Find solution for n ="; n; " and p =";p
If legendre_symbol(n, p) <> 1 Then Print n;" is not a quadratic residue" Print Continue For End If
If p = 2 OrElse Isprime(p, 15) = FALSE Then Print p;" is not a odd prime" Print Continue For End If
s = 0 : q = p -1 Do s += 1 q \= 2 Loop Until (q And 1) = 1
If s = 1 And (p Mod 4) = 3 Then r = pow_mod(n, ((p +1) \ 4), p) Print "Solution found:"; r; " and"; p - r Print Continue For End If
z = 1 Do z += 1 Loop Until legendre_symbol(z, p) = -1 c = pow_mod(z, q, p) r = pow_mod(n, (q +1) \ 2, p) t = pow_mod(n, q, p) m = s
Do i = 0 If (t Mod p) = 1 Then Print "Solution found:"; r; " and"; p - r Print Continue For End If
Do i += 1 If i >= m Then Continue For Loop Until pow_mod(t, 2 ^ i, p) = 1 b = pow_mod(c, (2 ^ (m - i -1)), p) r = mul_mod(r, b, p) c = mul_mod(b, b, p) t = mul_mod(t, c, p)' t = t * b ^ 2 m = i Loop
Next
Data 10, 13, 56, 101, 1030, 10009, 1032, 10009, 44402, 100049 Data 665820697, 1000000009, 881398088036, 1000000000039
' empty keyboard buffer While InKey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
Find solution for n = 10 and p = 13 Solution found: 7 and 6 Find solution for n = 56 and p = 101 Solution found: 37 and 64 Find solution for n = 1030 and p = 10009 Solution found: 1632 and 8377 Find solution for n = 1032 and p = 10009 1032 is not a quadratic residue Find solution for n = 44402 and p = 100049 Solution found: 30468 and 69581 Find solution for n = 665820697 and p = 1000000009 Solution found: 378633312 and 621366697 Find solution for n = 881398088036 and p = 1000000000039 Solution found: 791399408049 and 208600591990
GMP version
<lang freebasic>' version 12-04-2017 ' compile with: fbc -s console
- Include Once "gmp.bi"
Data "10", "13", "56", "101", "1030", "10009", "1032", "10009" Data "44402", "100049", "665820697", "1000000009" Data "881398088036", "1000000000039" Data "41660815127637347468140745042827704103445750172002" ' p = 10^50 + 577
' ------=< MAIN >=------
Dim As uLong k Dim As ZString Ptr zstr Dim As String n_str, p_str
Dim As Mpz_ptr b, c, i, m, n, p, q, r, s, t, z, tmp b = Allocate(Len(__Mpz_struct)) : Mpz_init(b) c = Allocate(Len(__Mpz_struct)) : Mpz_init(c) i = Allocate(Len(__Mpz_struct)) : Mpz_init(i) m = Allocate(Len(__Mpz_struct)) : Mpz_init(m) n = Allocate(Len(__Mpz_struct)) : Mpz_init(n) p = Allocate(Len(__Mpz_struct)) : Mpz_init(p) q = Allocate(Len(__Mpz_struct)) : Mpz_init(q) r = Allocate(Len(__Mpz_struct)) : Mpz_init(r) s = Allocate(Len(__Mpz_struct)) : Mpz_init(s) t = Allocate(Len(__Mpz_struct)) : Mpz_init(t) z = Allocate(Len(__Mpz_struct)) : Mpz_init(z) tmp = Allocate(Len(__Mpz_struct)) : Mpz_init(tmp)
For k = 1 To 8
Read n_str Mpz_set_str(n, n_str, 10) If k < 8 Then Read p_str Mpz_set_str(p, p_str, 10) Else p_str = "10^50 + 577" Mpz_set_str(p, "1" + String(50, "0"), 10) Mpz_add_ui(p, p, 577) End If
Print "Find solution for n = "; n_str; " and p = "; p_str
If Mpz_legendre(n, p) <> 1 Then Print n_str; " is not a quadratic residue" Print Continue For End If
If Mpz_tstbit(p, 0) = 0 OrElse Mpz_probab_prime_p(p, 20) = 0 Then Print p_str; "is not a odd prime" Print Continue For End If
Mpz_set_ui(s, 0) : Mpz_set(q, p) : Mpz_sub_ui(q, q, 1) ' q = p -1 Do Mpz_add_ui(s, s, 1) Mpz_fdiv_q_2exp(q, q, 1) Loop Until Mpz_tstbit(q, 0) = 1
If Mpz_cmp_ui(s, 1) = 0 Then If Mpz_tstbit(p, 1) = 1 Then Mpz_add_ui(tmp, p, 1) Mpz_fdiv_q_2exp(tmp, tmp, 2) ' tmp = p +1 \ 4 Mpz_powm(r, n, tmp, p) zstr = Mpz_get_str(0, 10, r) Print "Solution found: "; *zstr; Mpz_sub(r, p, r) zstr = Mpz_get_str(0, 10, r) Print " and "; *zstr Print Continue For End If End If
Mpz_set_ui(z, 1) Do Mpz_add_ui(z, z, 1) Loop Until Mpz_legendre(z, p) = -1 Mpz_powm(c, z, q, p) Mpz_add_ui(tmp, q, 1) Mpz_fdiv_q_2exp(tmp, tmp, 1) Mpz_powm(r, n, tmp, p) Mpz_powm(t, n, q, p) Mpz_set(m, s)
Do Mpz_set_ui(i, 0) Mpz_mod(tmp, t, p) If Mpz_cmp_ui(tmp, 1) = 0 Then zstr = Mpz_get_str(0, 10, r) Print "Solution found: "; *zstr; Mpz_sub(r, p, r) zstr = Mpz_get_str(0, 10, r) Print " and "; *zstr Print Continue For End If
Mpz_set_ui(q, 1) Do Mpz_add_ui(i, i, 1) If Mpz_cmp(i, m) >= 0 Then Continue For end if Mpz_mul_ui(q, q, 2) ' q = 2^i Mpz_powm(tmp, t, q, p) Loop Until Mpz_cmp_ui(tmp, 1) = 0
Mpz_set_ui(q, 2) Mpz_sub(tmp, m, i) : Mpz_sub_ui(tmp, tmp, 1) : Mpz_powm(tmp, q, tmp, p) Mpz_powm(b, c, tmp, p) Mpz_mul(r, r, b) : Mpz_mod(r, r, p) Mpz_mul(tmp, b, b) : Mpz_mod(c, tmp, p) Mpz_mul(tmp, t, c) : Mpz_mod(t, tmp, p) Mpz_set(m, i) Loop
Next
Mpz_clear(b) : Mpz_clear(c) : Mpz_clear(i) : Mpz_clear(m) Mpz_clear(n) : Mpz_clear(p) : Mpz_clear(q) : Mpz_clear(r) Mpz_clear(s) : Mpz_clear(t) : Mpz_clear(z) : Mpz_clear(tmp)
' empty keyboard buffer While InKey <> "" : Wend Print : Print "hit any key to end program" Sleep End</lang>
- Output:
Find solution for n = 10 and p = 13 Solution found: 7 and 6 Find solution for n = 56 and p = 101 Solution found: 37 and 64 Find solution for n = 1030 and p = 10009 Solution found: 1632 and 8377 Find solution for n = 1032 and p = 10009 1032 is not a quadratic residue Find solution for n = 44402 and p = 100049 Solution found: 30468 and 69581 Find solution for n = 665820697 and p = 1000000009 Solution found: 378633312 and 621366697 Find solution for n = 881398088036 and p = 1000000000039 Solution found: 791399408049 and 208600591990 Find solution for n = 41660815127637347468140745042827704103445750172002 and p = 10^50 + 577 Solution found: 32102985369940620849741983987300038903725266634508 and 67897014630059379150258016012699961096274733366069
Go
int
Implementation following Wikipedia, using similar variable names, and using the int type for simplicity. <lang go>package main
import "fmt"
// Arguments n, p as described in WP // If Legendre symbol != 1, ok return is false. Otherwise ok return is true, // R1 is WP return value R and for convenience R2 is p-R1. func ts(n, p int) (R1, R2 int, ok bool) {
// a^e mod p powModP := func(a, e int) int { s := 1 for ; e > 0; e-- { s = s * a % p } return s } // Legendre symbol, returns 1, 0, or -1 mod p -- that's 1, 0, or p-1. ls := func(a int) int { return powModP(a, (p-1)/2) } // argument validation if ls(n) != 1 { return 0, 0, false } // WP step 1, factor out powers two. // variables Q, S named as at WP. Q := p - 1 S := 0 for Q&1 == 0 { S++ Q >>= 1 } // WP step 1, direct solution if S == 1 { R1 = powModP(n, (p+1)/4) return R1, p - R1, true } // WP step 2, select z, assign c z := 2 for ; ls(z) != p-1; z++ { } c := powModP(z, Q) // WP step 3, assign R, t, M R := powModP(n, (Q+1)/2) t := powModP(n, Q) M := S // WP step 4, loop for { // WP step 4.1, termination condition if t == 1 { return R, p - R, true } // WP step 4.2, find lowest i... i := 0 for z := t; z != 1 && i < M-1; { z = z * z % p i++ } // WP step 4.3, using a variable b, assign new values of R, t, c, M b := c for e := M - i - 1; e > 0; e-- { b = b * b % p } R = R * b % p c = b * b % p // more convenient to compute c before t t = t * c % p M = i }
}
func main() {
fmt.Println(ts(10, 13)) fmt.Println(ts(56, 101)) fmt.Println(ts(1030, 10009)) fmt.Println(ts(1032, 10009)) fmt.Println(ts(44402, 100049))
}</lang>
- Output:
7 6 true 37 64 true 1632 8377 true 0 0 false 30468 69581 true
big.Int
For the extra credit, we use big.Int from the math/big package of the Go standard library. While the method call syntax is not as easy on the eyes as operator syntax, the package provides modular exponentiation and even the Legendre symbol as the Jacobi function. <lang go>package main
import (
"fmt" "math/big"
)
func ts(n, p big.Int) (R1, R2 big.Int, ok bool) {
if big.Jacobi(&n, &p) != 1 { return } var one, Q big.Int one.SetInt64(1) Q.Sub(&p, &one) S := 0 for Q.Bit(0) == 0 { S++ Q.Rsh(&Q, 1) } if S == 1 { R1.Exp(&n, R1.Rsh(R1.Add(&p, &one), 2), &p) R2.Sub(&p, &R1) return R1, R2, true } var z, c big.Int for z.SetInt64(2); big.Jacobi(&z, &p) != -1; z.Add(&z, &one) { } c.Exp(&z, &Q, &p) var R, t big.Int R.Exp(&n, R.Rsh(R.Add(&Q, &one), 1), &p) t.Exp(&n, &Q, &p) M := S for { if t.Cmp(&one) == 0 { R2.Sub(&p, &R) return R, R2, true } i := 0 // reuse z as a scratch variable for z.Set(&t); z.Cmp(&one) != 0 && i < M-1; { z.Mod(z.Mul(&z, &z), &p) i++ } // and instead of a new scratch variable b, continue using z z.Set(&c) for e := M - i - 1; e > 0; e-- { z.Mod(z.Mul(&z, &z), &p) } R.Mod(R.Mul(&R, &z), &p) c.Mod(c.Mul(&z, &z), &p) t.Mod(t.Mul(&t, &c), &p) M = i }
}
func main() {
var n, p big.Int n.SetInt64(665820697) p.SetInt64(1000000009) R1, R2, ok := ts(n, p) fmt.Println(&R1, &R2, ok)
n.SetInt64(881398088036) p.SetInt64(1000000000039) R1, R2, ok = ts(n, p) fmt.Println(&R1, &R2, ok) n.SetString("41660815127637347468140745042827704103445750172002", 10) p.SetString("100000000000000000000000000000000000000000000000577", 10) R1, R2, ok = ts(n, p) fmt.Println(&R1) fmt.Println(&R2)
}</lang>
- Output:
378633312 621366697 true 791399408049 208600591990 true 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069
Library
It gets better; the library has a ModSqrt function that uses Tonelli-Shanks internally. Output is same as above. <lang go>package main
import (
"fmt" "math/big"
)
func main() {
var n, p, R1, R2 big.Int n.SetInt64(665820697) p.SetInt64(1000000009) R1.ModSqrt(&n, &p) R2.Sub(&p, &R1) fmt.Println(&R1, &R2)
n.SetInt64(881398088036) p.SetInt64(1000000000039) R1.ModSqrt(&n, &p) R2.Sub(&p, &R1) fmt.Println(&R1, &R2)
n.SetString("41660815127637347468140745042827704103445750172002", 10) p.SetString("100000000000000000000000000000000000000000000000577", 10) R1.ModSqrt(&n, &p) R2.Sub(&p, &R1) fmt.Println(&R1) fmt.Println(&R2)
}</lang>
J
Implementation:
<lang J>leg=: dyad define
x (y&|)@^ (y-1)%2
)
tosh=:dyad define
assert. 1=1 p: y [ 'y must be prime' assert. 1=x leg y [ 'x must be square mod y' pow=. y&|@^ if. 1=m=. {.1 q: y-1 do. r=. x pow (y+1)%4 else. z=. 1x while. 1>: z leg y do. z=.z+1 end. c=. z pow q=. (y-1)%2^m r=. x pow (q+1)%2 t=. x pow q while. t~:1 do. n=. t i=. 0 whilst. 1~:n do. n=. n pow 2 i=. i+1 end. r=. y|r*b=. c pow 2^m-i+1 m=. i t=. y|t*c=. b pow 2 end. end. y|(,-)r
)</lang>
Task examples:
<lang J> 10 tosh 13 7 6
56 tosh 101
37 64
1030 tosh 10009
1632 8377
1032 tosh 10009
|assertion failure: tosh | 1=x leg y['x must be square mod y'
44402 tosh 100049
30468 69581
665820697x tosh 1000000009x
378633312 621366697
881398088036 tosh 1000000000039x
791399408049 208600591990
41660815127637347468140745042827704103445750172002x tosh (10^50x)+577
32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069</lang>
Java
<lang Java>import java.math.BigInteger; import java.util.List; import java.util.Map; import java.util.function.BiFunction; import java.util.function.Function;
public class TonelliShanks {
private static final BigInteger ZERO = BigInteger.ZERO; private static final BigInteger ONE = BigInteger.ONE; private static final BigInteger TEN = BigInteger.TEN; private static final BigInteger TWO = BigInteger.valueOf(2); private static final BigInteger FOUR = BigInteger.valueOf(4);
private static class Solution { private BigInteger root1; private BigInteger root2; private boolean exists;
Solution(BigInteger root1, BigInteger root2, boolean exists) { this.root1 = root1; this.root2 = root2; this.exists = exists; } }
private static Solution ts(Long n, Long p) { return ts(BigInteger.valueOf(n), BigInteger.valueOf(p)); }
private static Solution ts(BigInteger n, BigInteger p) { BiFunction<BigInteger, BigInteger, BigInteger> powModP = (BigInteger a, BigInteger e) -> a.modPow(e, p); Function<BigInteger, BigInteger> ls = (BigInteger a) -> powModP.apply(a, p.subtract(ONE).divide(TWO));
if (!ls.apply(n).equals(ONE)) return new Solution(ZERO, ZERO, false);
BigInteger q = p.subtract(ONE); BigInteger ss = ZERO; while (q.and(ONE).equals(ZERO)) { ss = ss.add(ONE); q = q.shiftRight(1); }
if (ss.equals(ONE)) { BigInteger r1 = powModP.apply(n, p.add(ONE).divide(FOUR)); return new Solution(r1, p.subtract(r1), true); }
BigInteger z = TWO; while (!ls.apply(z).equals(p.subtract(ONE))) z = z.add(ONE); BigInteger c = powModP.apply(z, q); BigInteger r = powModP.apply(n, q.add(ONE).divide(TWO)); BigInteger t = powModP.apply(n, q); BigInteger m = ss;
while (true) { if (t.equals(ONE)) return new Solution(r, p.subtract(r), true); BigInteger i = ZERO; BigInteger zz = t; while (!zz.equals(BigInteger.ONE) && i.compareTo(m.subtract(ONE)) < 0) { zz = zz.multiply(zz).mod(p); i = i.add(ONE); } BigInteger b = c; BigInteger e = m.subtract(i).subtract(ONE); while (e.compareTo(ZERO) > 0) { b = b.multiply(b).mod(p); e = e.subtract(ONE); } r = r.multiply(b).mod(p); c = b.multiply(b).mod(p); t = t.multiply(c).mod(p); m = i; } }
public static void main(String[] args) { List<Map.Entry<Long, Long>> pairs = List.of( Map.entry(10L, 13L), Map.entry(56L, 101L), Map.entry(1030L, 10009L), Map.entry(1032L, 10009L), Map.entry(44402L, 100049L), Map.entry(665820697L, 1000000009L), Map.entry(881398088036L, 1000000000039L) );
for (Map.Entry<Long, Long> pair : pairs) { Solution sol = ts(pair.getKey(), pair.getValue()); System.out.printf("n = %s\n", pair.getKey()); System.out.printf("p = %s\n", pair.getValue()); if (sol.exists) { System.out.printf("root1 = %s\n", sol.root1); System.out.printf("root2 = %s\n", sol.root2); } else { System.out.println("No solution exists"); } System.out.println(); }
BigInteger bn = new BigInteger("41660815127637347468140745042827704103445750172002"); BigInteger bp = TEN.pow(50).add(BigInteger.valueOf(577)); Solution sol = ts(bn, bp); System.out.printf("n = %s\n", bn); System.out.printf("p = %s\n", bp); if (sol.exists) { System.out.printf("root1 = %s\n", sol.root1); System.out.printf("root2 = %s\n", sol.root2); } else { System.out.println("No solution exists"); } }
}</lang>
- Output:
n = 10 p = 13 root1 = 7 root2 = 6 n = 56 p = 101 root1 = 37 root2 = 64 n = 1030 p = 10009 root1 = 1632 root2 = 8377 n = 1032 p = 10009 No solution exists n = 44402 p = 100049 root1 = 30468 root2 = 69581 n = 665820697 p = 1000000009 root1 = 378633312 root2 = 621366697 n = 881398088036 p = 1000000000039 root1 = 791399408049 root2 = 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 root1 = 32102985369940620849741983987300038903725266634508 root2 = 67897014630059379150258016012699961096274733366069
Kotlin
<lang scala>// version 1.1.3
import java.math.BigInteger
data class Solution(val root1: BigInteger, val root2: BigInteger, val exists: Boolean)
val bigZero = BigInteger.ZERO val bigOne = BigInteger.ONE val bigTwo = BigInteger.valueOf(2L) val bigFour = BigInteger.valueOf(4L) val bigTen = BigInteger.TEN
fun ts(n: Long, p: Long) = ts(BigInteger.valueOf(n), BigInteger.valueOf(p))
fun ts(n: BigInteger, p: BigInteger): Solution {
fun powModP(a: BigInteger, e: BigInteger) = a.modPow(e, p)
fun ls(a: BigInteger) = powModP(a, (p - bigOne) / bigTwo)
if (ls(n) != bigOne) return Solution(bigZero, bigZero, false) var q = p - bigOne var ss = bigZero while (q.and(bigOne) == bigZero) { ss = ss + bigOne q = q.shiftRight(1) }
if (ss == bigOne) { val r1 = powModP(n, (p + bigOne) / bigFour) return Solution(r1, p - r1, true) }
var z = bigTwo while (ls(z) != p - bigOne) z = z + bigOne var c = powModP(z, q) var r = powModP(n, (q + bigOne) / bigTwo) var t = powModP(n, q) var m = ss
while (true) { if (t == bigOne) return Solution(r, p - r, true) var i = bigZero var zz = t while (zz != bigOne && i < m - bigOne) { zz = zz * zz % p i = i + bigOne } var b = c var e = m - i - bigOne while (e > bigZero) { b = b * b % p e = e - bigOne } r = r * b % p c = b * b % p t = t * c % p m = i }
}
fun main(args: Array<String>) {
val pairs = listOf<Pair<Long, Long>>( 10L to 13L, 56L to 101L, 1030L to 10009L, 1032L to 10009L, 44402L to 100049L, 665820697L to 1000000009L, 881398088036L to 1000000000039L )
for (pair in pairs) { val (n, p) = pair val (root1, root2, exists) = ts(n, p) println("n = $n") println("p = $p") if (exists) { println("root1 = $root1") println("root2 = $root2") } else println("No solution exists") println() }
val bn = BigInteger("41660815127637347468140745042827704103445750172002") val bp = bigTen.pow(50) + BigInteger.valueOf(577L) val (broot1, broot2, bexists) = ts(bn, bp) println("n = $bn") println("p = $bp") if (bexists) { println("root1 = $broot1") println("root2 = $broot2") } else println("No solution exists")
}</lang>
- Output:
n = 10 p = 13 root1 = 7 root2 = 6 n = 56 p = 101 root1 = 37 root2 = 64 n = 1030 p = 10009 root1 = 1632 root2 = 8377 n = 1032 p = 10009 No solution exists n = 44402 p = 100049 root1 = 30468 root2 = 69581 n = 665820697 p = 1000000009 root1 = 378633312 root2 = 621366697 n = 881398088036 p = 1000000000039 root1 = 791399408049 root2 = 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 root1 = 32102985369940620849741983987300038903725266634508 root2 = 67897014630059379150258016012699961096274733366069
Perl 6
Translation of the Wikipedia pseudocode, heavily influenced by Sidef and Python.
<lang perl6># Legendre operator (𝑛│𝑝) sub infix:<│> (Int \𝑛, Int \𝑝 where 𝑝.is-prime && (𝑝 != 2)) {
given 𝑛.expmod( (𝑝-1) div 2, 𝑝 ) { when 0 { 0 } when 1 { 1 } default { -1 } }
}
sub tonelli-shanks ( \𝑛, \𝑝 where (𝑛│𝑝) > 0 ) {
my $𝑄 = 𝑝 - 1; my $𝑆 = 0; $𝑄 +>= 1 and $𝑆++ while $𝑄 %% 2; return 𝑛.expmod((𝑝+1) div 4, 𝑝) if $𝑆 == 1; my $𝑐 = ((2..𝑝).first: (*│𝑝) < 0).expmod($𝑄, 𝑝); my $𝑅 = 𝑛.expmod( ($𝑄+1) +> 1, 𝑝 ); my $𝑡 = 𝑛.expmod( $𝑄, 𝑝 ); while ($𝑡-1) % 𝑝 { my $b; my $𝑡2 = $𝑡² % 𝑝; for 1 .. $𝑆 { if ($𝑡2-1) %% 𝑝 { $b = $𝑐.expmod(1 +< ($𝑆-1-$_), 𝑝); $𝑆 = $_; last; } $𝑡2 = $𝑡2² % 𝑝; } $𝑅 = ($𝑅 * $b) % 𝑝; $𝑐 = $b² % 𝑝; $𝑡 = ($𝑡 * $𝑐) % 𝑝; } $𝑅;
}
my @tests = (
(10, 13), (56, 101), (1030, 10009), (1032, 10009), (44402, 100049), (665820697, 1000000009), (881398088036, 1000000000039), (41660815127637347468140745042827704103445750172002, 100000000000000000000000000000000000000000000000577)
);
for @tests -> ($n, $p) { try my $t = tonelli-shanks($n, $p); say "No solution for ({$n}, {$p})." and next if !$t or ($t² - $n) % $p; say "Roots of $n are ($t, {$p-$t}) mod $p";
}</lang>
- Output:
Roots of 10 are (7, 6) mod 13 Roots of 56 are (37, 64) mod 101 Roots of 1030 are (1632, 8377) mod 10009 No solution for (1032, 10009). Roots of 44402 are (30468, 69581) mod 100049 Roots of 665820697 are (378633312, 621366697) mod 1000000009 Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039 Roots of 41660815127637347468140745042827704103445750172002 are (32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069) mod 100000000000000000000000000000000000000000000000577
PicoLisp
<lang PicoLisp># from @lib/rsa.l (de **Mod (X Y N)
(let M 1 (loop (when (bit? 1 Y) (setq M (% (* M X) N)) ) (T (=0 (setq Y (>> 1 Y))) M ) (setq X (% (* X X) N)) ) ) )
(de legendre (N P)
(**Mod N (/ (dec P) 2) P) )
(de ts (N P)
(and (=1 (legendre N P)) (let (Q (dec P) S 0 Z 0 C 0 R 0 D 0 M 0 B 0 I 0 ) (until (bit? 1 Q) (setq Q (>> 1 Q)) (inc 'S) ) (if (=1 S) (list (setq @@ (**Mod N (/ (inc P) 4) P)) (- P @@) ) (setq Z 2) (until (= (legendre Z P) (dec P)) (inc 'Z) ) (setq C (**Mod Z Q P) R (**Mod N (/ (inc Q) 2) P) D (**Mod N Q P) M S ) (until (=1 D) (zero I) (for (Z D (and (<> Z 1) (< I (dec M))) (setq Z (% (* Z Z) P)) ) (inc 'I) ) (setq B C) (for (Z (- M I 1) (> Z 0) (dec Z) ) (setq B (% (* B B) P)) ) (setq R (% (* R B) P) C (% (* B B) P) D (% (* D C) P) M I ) ) (list R (- P R)) ) ) ) )
(println (ts 10 13)) (println (ts 56 101)) (println (ts 1030 10009)) (println (ts 1032 10009)) (println (ts 44402 100049)) (println (ts 665820697 1000000009)) (println (ts 881398088036 1000000000039)) (println (ts 41660815127637347468140745042827704103445750172002 (+ (** 10 50) 577)))</lang>
- Output:
(7 6) (37 64) (1632 8377) NIL (30468 69581) (378633312 621366697) (791399408049 208600591990) (32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069)
Python
<lang python>def legendre(a, p):
return pow(a, (p - 1) // 2, p)
def tonelli(n, p):
assert legendre(n, p) == 1, "not a square (mod p)" q = p - 1 s = 0 while q % 2 == 0: q //= 2 s += 1 if s == 1: return pow(n, (p + 1) // 4, p) for z in range(2, p): if p - 1 == legendre(z, p): break c = pow(z, q, p) r = pow(n, (q + 1) // 2, p) t = pow(n, q, p) m = s t2 = 0 while (t - 1) % p != 0: t2 = (t * t) % p for i in range(1, m): if (t2 - 1) % p == 0: break t2 = (t2 * t2) % p b = pow(c, 1 << (m - i - 1), p) r = (r * b) % p c = (b * b) % p t = (t * c) % p m = i return r
if __name__ == '__main__':
ttest = [(10, 13), (56, 101), (1030, 10009), (44402, 100049),
(665820697, 1000000009), (881398088036, 1000000000039),
(41660815127637347468140745042827704103445750172002, 10**50 + 577)] for n, p in ttest: r = tonelli(n, p) assert (r * r - n) % p == 0 print("n = %d p = %d" % (n, p)) print("\t roots : %d %d" % (r, p - r))</lang>
- Output:
n = 10 p = 13 roots : 7 6 n = 56 p = 101 roots : 37 64 n = 1030 p = 10009 roots : 1632 8377 n = 44402 p = 100049 roots : 30468 69581 n = 665820697 p = 1000000009 roots : 378633312 621366697 n = 881398088036 p = 1000000000039 roots : 791399408049 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 roots : 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069
Racket
<lang racket>#lang racket
(require math/number-theory)
(define (Legendre a p)
(modexpt a (quotient (sub1 p) 2)))
(define (Tonelli n p (err (λ (n p) (error "not a square (mod p)" (list n p)))))
(with-modulus p (unless (= 1 (Legendre n p)) (err n p))
(define-values (q s) (let even?-q-loop ((q (sub1 p)) (s 0)) (if (even? q) (even?-q-loop (quotient q 2) (add1 s)) (values q s)))) (cond [(= s 1) (modexpt n (/ (add1 p) 4))] [else (define z (for/first ((z (in-range 2 p)) #:when (= (sub1 p) (Legendre z p))) z)) (let loop ((c (modexpt z q)) (r (modexpt n (quotient (add1 q) 2))) (t (modexpt n q)) (m s)) (cond [(mod= 1 t) r] [else (define-values (t2 m′) (for/fold ((t2 (modsqr t)) (i 1)) ((j (in-range 1 m)) #:final (mod= t2 1)) (values (modsqr t2) j))) (define b (modexpt c (expt 2 (- m m′ 1)))) (define c′ (modsqr b)) (loop c′ (mod* r b) (mod* t c′) m′)]))])))
(module+ test
(require rackunit)
(define ttest `((10 13) (56 101) (1030 10009) (44402 100049) (665820697 1000000009) (881398088036 1000000000039) (41660815127637347468140745042827704103445750172002 ,(+ #e1e50 577))))
(define (task ttest) (for ((test ttest)) (define n (first test)) (define p (second test)) (define r (Tonelli n p)) (printf "n = ~a p = ~a~% roots : ~a ~a~%" n p r (- p r))))
(task ttest)
(check-exn exn:fail? (λ () (Tonelli 1032 1009))))</lang>
- Output:
n = 10 p = 13 roots : 7 6 n = 56 p = 101 roots : 37 64 n = 1030 p = 10009 roots : 1632 8377 n = 44402 p = 100049 roots : 30468 69581 n = 665820697 p = 1000000009 roots : 378633312 621366697 n = 881398088036 p = 1000000000039 roots : 791399408049 208600591990 n = 41660815127637347468140745042827704103445750172002 p = 100000000000000000000000000000000000000000000000577 roots : 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069
Sidef
<lang ruby>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}"
}</lang>
- Output:
Roots of 10 are (7, 6) mod 13 Roots of 56 are (37, 64) mod 101 Roots of 1030 are (1632, 8377) mod 10009 Roots of 44402 are (30468, 69581) mod 100049 Roots of 665820697 are (378633312, 621366697) mod 1000000009 Roots of 881398088036 are (791399408049, 208600591990) mod 1000000000039 Roots of 41660815127637347468140745042827704103445750172002 are (32102985369940620849741983987300038903725266634508, 67897014630059379150258016012699961096274733366069) mod 100000000000000000000000000000000000000000000000577
zkl
<lang zkl>var BN=Import("zklBigNum"); fcn modEq(a,b,p) { (a-b)%p==0 } fcn legendre(a,p){ a.powm((p - 1)/2,p) }
fcn tonelli(n,p){ //(BigInt,Int|BigInt)
_assert_(legendre(n,p)==1, "not a square (mod p)"+vm.arglist); q,s:=p-1,0; while(q.isEven){ q/=2; s+=1; } if(s==1) return(n.powm((p+1)/4,p)); z:=[BN(2)..p].filter1('wrap(z){ legendre(z,p)==(p-1) }); c,r,t,m,t2:=z.powm(q,p), n.powm((q+1)/2,p), n.powm(q,p), s, 0; while(not modEq(t,1,p)){ t2=(t*t)%p; i:=1; while(not modEq(t2,1,p)){ i+=1; t2=(t2*t2)%p; } // assert(i<m) b:=c.powm(BN(1).shiftLeft(m-i-1), p); r,c,t,m = (r*b)%p, (b*b)%p, (t*c)%p, i; } r
}</lang> <lang zkl>ttest:=T(T(10,13), T(56,101), T(1030,10009), T(44402,100049),
T(665820697,1000000009), T(881398088036,1000000000039), T("41660815127637347468140745042827704103445750172002", BN(10).pow(50) + 577), T(1032,10009) );
foreach n,p in (ttest){ n=BN(n);
r:=tonelli(n,p); assert((r*r-n)%p == 0,"(r*r-n)%p == 0 : %s,%s,%s-->%s".fmt(r,n,p,(r*r-n)%p)); println("n=%d p=%d".fmt(n,p)); println(" roots: %d %d".fmt(r, p-r));
}</lang>
- Output:
n=10 p=13 roots: 7 6 n=56 p=101 roots: 37 64 n=1030 p=10009 roots: 1632 8377 n=44402 p=100049 roots: 30468 69581 n=665820697 p=1000000009 roots: 378633312 621366697 n=881398088036 p=1000000000039 roots: 791399408049 208600591990 n=41660815127637347468140745042827704103445750172002 p=100000000000000000000000000000000000000000000000577 roots: 32102985369940620849741983987300038903725266634508 67897014630059379150258016012699961096274733366069 VM#1 caught this unhandled exception: AssertionError : not a square (mod p)L(1032,10009) Stack trace for VM#1 (): bbb.assert addr:13 args(2) reg(0) bbb.tonelli addr:29 args(2) reg(10) R ...