; tonelli-shanks algorithm
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
(define prime?
(let ((seed 3141592654))
(lambda (n)
(define (rand)
(set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))
(+ (quotient (* seed (- n 2)) 4294967296) 2))
(define (expm b e m)
(define (times x y) (modulo (* x y) m))
(let loop ((b b) (e e) (r 1))
(if (zero? e) r
(loop (times b b) (quotient e 2)
(if (odd? e) (times b r) r)))))
(define (spsp? n a)
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d)
(let ((t (expm a d n)))
(if (or (= t 1) (= t (- n 1))) #t
(do ((s (- s 1) (- s 1))
(t (expm t 2 n) (expm t 2 n)))
((or (zero? s) (= t (- n 1)))
(positive? s))))))))
(if (not (integer? n))
(error 'prime? "must be integer")
(if (< n 2) #f
(do ((a (rand) (rand)) (k 25 (- k 1)))
((or (zero? k) (not (spsp? n a)))
(zero? k))))))))
(define (tonelli a p)
(define (expp b e) (expm b e p))
(define (modp n) (modulo n p))
(let* ((p-1 (- p 1)) (p-1/2 (/ p-1 2)))
(if (even? p) (error 'tonelli "must be odd"))
(if (not (prime? p)) (error 'tonelli "must be prime"))
(if (< p 3) (error 'tonelli "must be greater than two"))
(if (< 1 (gcd a p)) (error 'tonelli "must be coprime"))
(if (= (expp a p-1/2) -1)
(error 'tonelli "must be quadratic residue"))
(let loop ((s p-1) (e 0))
(if (even? s) (loop (/ s 2) (+ e 1))
(let loop ((n 2))
(if (not (= (expp n p-1/2) p-1)) (loop (+ n 1))
(let loop1 ((x (expp a (/ (+ s 1) 2)))
(b (expp a s)) (g (expp n s)) (r e))
(let loop2 ((m 0) (mm 1))
(if (not (= (expp b mm) 1)) (loop2 (+ m 1) (* mm 2))
(if (zero? m) (list x (- p x)) ; found result
(loop1 (modp (* x (expp g (expt 2 (- r m 1)))))
(modp (* b (expp g (expt 2 (- r m)))))
(expp g (expt 2 (- r m))) m)))))))))))
(display (tonelli 2 113)) (newline)
(display (tonelli 5 40961)) (newline)
(display (tonelli 2 360027784083079948259017962255826129)) (newline)