42 lines
958 B
Common Lisp
42 lines
958 B
Common Lisp
(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)))))
|