[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/CAPseLXK    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Nov 23:
; 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)


Output:
1
2
3
(62 51)
(19424 21537)
(162244492740221711333411667492080568 197783291342858236925606294763745561)


Create a new paste based on this one


Comments: