Tonelli-Shanks algorithm: Difference between revisions
Thundergnat (talk | contribs) (→{{header|Perl 6}}: Add Perl 6) |
m (→{{header|Sidef}}: code simplifications) |
||
| Line 542: | Line 542: | ||
legendre(n, p) == 1 || die "not a square (mod p)" |
legendre(n, p) == 1 || die "not a square (mod p)" |
||
var q = p-1 |
var q = p-1 |
||
var s = |
var s = valuation(q, 2) |
||
| ⚫ | |||
while (q.is_even) { |
|||
| ⚫ | |||
q >>= 1 |
|||
s += 1 |
|||
} |
|||
if (s == 1) { |
|||
| ⚫ | |||
} |
|||
var c |
|||
for z in (2 ..^ p) { |
|||
| ⚫ | |||
c = powmod(z, q, p) |
|||
break |
|||
} |
|||
} |
|||
var r = powmod(n, (q + 1) >> 1, p) |
var r = powmod(n, (q + 1) >> 1, p) |
||
var t = powmod(n, q, p) |
var t = powmod(n, q, p) |
||
Revision as of 13:58, 25 October 2016
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
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)
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>
Perl 6
Translation of the Wikipedia pseudocode, heavily influenced by Sidef and Python.
<lang perl6># Legendre operator (𝑛│𝑝) sub infix:<│> (Int \𝑛, Int \𝑝 where (𝑝.is-prime and ?(𝑝 != 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 $𝑐;
for 2 .. 𝑝 {
$𝑐 = .expmod($𝑄, 𝑝) and last if ($_│𝑝) < 0;
}
my $𝑅 = 𝑛.expmod( ($𝑄+1) +> 1, 𝑝 );
my $𝑡 = 𝑛.expmod( $𝑄, 𝑝 );
my $𝑀 = $𝑆;
my $b;
while (($𝑡-1) % 𝑝) {
my $𝑡2 = $𝑡² % 𝑝;
for 1 .. $𝑀 {
if ($𝑡2-1) %% 𝑝 {
$b = $𝑐.expmod( 1 +< ( $𝑀 - 1 - $_ ), 𝑝 );
$𝑀 = $_;
last;
}
$𝑡2 = $𝑡2² % 𝑝;
}
$𝑅 = ($𝑅 * $b) % 𝑝;
$𝑐 = $b² % 𝑝;
$𝑡 = ($𝑡 * $𝑐) % 𝑝;
}
$𝑅;
}
my @tests = (
(10, 13),
(56, 101),
(1030, 10009),
(44402, 100049),
(665820697, 1000000009),
(881398088036, 1000000000039),
(41660815127637347468140745042827704103445750172002,
100000000000000000000000000000000000000000000000577)
);
for @tests -> ($n, $p) {
my $t = tonelli-shanks($n, $p);
die if ($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 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
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
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 ...