61 lines
1.7 KiB
Racket
61 lines
1.7 KiB
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))))
|