; your code goes here
(define prime? ; strong pseudoprime to prime bases less than 100
(let* ((ps (list 2 3 5 7 11 13 17 19 23 29 31 37
41 43 47 53 59 61 67 71 73 79 83 89 97))
(p100 (apply * ps)))
(lambda (n)
(define (expm b e m)
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (modulo (* b b) m) (quotient e 2)
(if (odd? e) (modulo (* b x) m) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d) (if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
(< r s)))))))
(if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
(define (euclid x y)
(let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y))
(if (zero? w) (values a b g)
(let ((q (quotient g w)))
(loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w)))))))
(define (inverse x m)
(if (not (= (gcd x m) 1))
(error 'inverse "divisor must be coprime to modulus")
(call-with-values
(lambda () (euclid x m))
(lambda (a b g) (modulo a m)))))
(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 (jacobi a n)
(if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
(error 'jacobi "modulus must be positive odd integer")
(let jacobi ((a a) (n n))
(cond ((= a 0) 0)
((= a 1) 1)
((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
((< n a) (jacobi (modulo a n) n))
((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
(else (jacobi n a))))))
(define (mod-fact n m)
(if (<= m n) 0
(let loop ((k 2) (p 1))
(if (zero? p) 0
(if (< n k) p
(loop (+ k 1) (modulo (* p k) m)))))))
(define (mod-sqrt a p)
(define (both n) (list n (- p n)))
(cond ((not (and (odd? p) (prime? p)))
(error 'mod-sqrt "modulus must be an odd prime"))
((not (= (jacobi a p) 1))
(error 'mod-sqrt "must be a quadratic residual"))
(else (let ((a (modulo a p)))
(case (modulo p 8)
((3 7) (both (expm a (/ (+ p 1) 4) p)))
((5) (let* ((x (expm a (/ (+ p 3) 8) p))
(c (expm x 2 p)))
(if (= a c) (both x)
(both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p)))))
(else (let* ((d (let loop ((d 2))
(if (= (jacobi d p) -1) d
(loop (+ d 1)))))
(s (let loop ((p (- p 1)) (s 0))
(if (odd? p) s
(loop (quotient p 2) (+ s 1)))))
(t (quotient (- p 1) (expt 2 s)))
(big-a (expm a t p))
(big-d (expm d t p))
(m (let loop ((i 0) (m 0))
(cond ((= i s) m)
((= (- p 1)
(expm (* big-a (expm big-d m p))
(expt 2 (- s 1 i)) p))
(loop (+ i 1) (+ m (expt 2 i))))
(else (loop (+ i 1) m))))))
(both (modulo (* (expm a (/ (+ t 1) 2) p)
(expm big-d (/ m 2) p)) p)))))))))
(define-syntax (with-modulus stx)
(syntax-case stx ()
((with-modulus e expr ...)
(with-syntax ((modulus (datum->syntax-object (syntax with-modulus) 'modulus))
(== (datum->syntax-object (syntax with-modulus) '== ))
(+ (datum->syntax-object (syntax with-modulus) '+ ))
(- (datum->syntax-object (syntax with-modulus) '- ))
(* (datum->syntax-object (syntax with-modulus) '* ))
(/ (datum->syntax-object (syntax with-modulus) '/ ))
(^ (datum->syntax-object (syntax with-modulus) '^ ))
(! (datum->syntax-object (syntax with-modulus) '! ))
(sqrt (datum->syntax-object (syntax with-modulus) 'sqrt )))
(syntax (letrec ((fold (lambda (op base xs)
(if (null? xs) base
(fold op (op base (car xs)) (cdr xs))))))
(let* ((modulus e)
(mod (lambda (x)
(if (not (integer? x))
(error 'with-modulus "all arguments must be integers")
(modulo x modulus))))
(== (lambda (x y) (= (mod x) (mod y))))
(+ (lambda xs (fold (lambda (x y) (mod (+ x (mod y)))) 0 xs)))
(- (lambda (x . xs)
(if (null? xs)
(mod (- 0 x))
(fold (lambda (x y) (mod (- x (mod y)))) x xs))))
(* (lambda xs (fold (lambda (x y) (mod (* x (mod y)))) 1 xs)))
(/ (lambda (x . xs)
(if (null? xs)
(inverse x e)
(fold (lambda (x y) (* x (inverse y e))) x xs))))
(^ (lambda (base exp) (expm base exp e)))
(! (lambda (n) (mod-fact n modulus)))
(sqrt (lambda (x) (mod-sqrt x e))))
expr ...)))))))
(define (twin? m)
(with-modulus (* m (+ m 2))
(== (* 4 (+ (! (- m 1)) 1))
(- m))))
(define (range . args)
(case (length args)
((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
((3) (let ((le? (if (negative? (caddr args)) >= <=)))
(let loop ((x(car args)) (xs '()))
(if (le? (cadr args) x)
(reverse xs)
(loop (+ x (caddr args)) (cons x xs))))))
(else (error 'range "unrecognized arguments"))))
(display (filter twin? (range 3 1000 2))) (newline)