Tonelli-Shanks algorithm

Tonelli-Shanks algorithm is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Tonelli–Shanks algorithm

Solve x² ≡ n (mod p)

In computational number theory, the Tonelli–Shanks algorithm is a technique for solving an equation of the form x² ≡ n (mod p), where 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 is 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:

• Modular exponentiation
• Cipolla's algorithm

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)

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

  8218 tosh 10007

9872 135

  8219 tosh 10007

|assertion failure: tosh | 1=x leg y['x must be square mod y'

  331575 tosh 1000003

144161 855842

  665165880x tosh 1000000007x

475131702 524868305

  881398088036x tosh 1000000000039x

791399408049 208600591990

  34035243914635549601583369544560650254325084643201x tosh (10^50x) + 151

82563118828090362261378993957450213573687113690751 17436881171909637738621006042549786426312886309400</lang>