fork download
copy 
  1. ; your code goes here
  2.  
  3. (define prime? ; strong pseudoprime to prime bases less than 100
  4.   (let* ((ps (list 2 3 5 7 11 13 17 19 23 29 31 37
  5.              41 43 47 53 59 61 67 71 73 79 83 89 97))
  6.          (p100 (apply * ps)))
  7.     (lambda (n)
  8.       (define (expm b e m)
  9.         (let loop ((b b) (e e) (x 1))
  10.           (if (zero? e) x
  11.             (loop (modulo (* b b) m) (quotient e 2)
  12.                   (if (odd? e) (modulo (* b x) m) x)))))
  13.       (define (spsp? n a) ; #t if n is a strong pseudoprime base a
  14.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  15.             ((odd? d) (if (= (expm a d n) 1) #t
  16.               (do ((r 0 (+ r 1)))
  17.                   ((or (= (expm a (* d (expt 2 r)) n) (- n 1)) (= r s))
  18.                     (< r s)))))))
  19.       (if (< n 2) #f (if (< 1 (gcd n p100)) (if (member n ps) #t #f)
  20.         (do ((ps ps (cdr ps)))
  21.             ((or (null? ps) (not (spsp? n (car ps)))) (null? ps))))))))
  22.  
  23. (define (euclid x y)
  24.   (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y))
  25.     (if (zero? w) (values a b g)
  26.       (let ((q (quotient g w)))
  27.         (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w)))))))
  28.  
  29. (define (inverse x m)
  30.   (if (not (= (gcd x m) 1))
  31.       (error 'inverse "divisor must be coprime to modulus")
  32.       (call-with-values
  33.         (lambda () (euclid x m))
  34.         (lambda (a b g) (modulo a m)))))
  35.  
  36. (define (expm b e m)
  37.   (define (m* x y) (modulo (* x y) m))
  38.   (cond ((zero? e) 1)
  39.         ((even? e) (expm (m* b b) (/ e 2) m))
  40.         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
  41.  
  42. (define (jacobi a n)
  43.   (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
  44.       (error 'jacobi "modulus must be positive odd integer")
  45.       (let jacobi ((a a) (n n))
  46.         (cond ((= a 0) 0)
  47.               ((= a 1) 1)
  48.               ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
  49.               ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
  50.               ((< n a) (jacobi (modulo a n) n))
  51.               ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
  52.               (else (jacobi n a))))))
  53.  
  54. (define (mod-fact n m)
  55.   (if (<= m n) 0
  56.     (let loop ((k 2) (p 1))
  57.       (if (zero? p) 0
  58.         (if (< n k) p
  59.           (loop (+ k 1) (modulo (* p k) m)))))))
  60.  
  61. (define (mod-sqrt a p)
  62.   (define (both n) (list n (- p n)))
  63.   (cond ((not (and (odd? p) (prime? p)))
  64.           (error 'mod-sqrt "modulus must be an odd prime"))
  65.         ((not (= (jacobi a p) 1))
  66.           (error 'mod-sqrt "must be a quadratic residual"))
  67.         (else (let ((a (modulo a p)))
  68.                 (case (modulo p 8)
  69.                   ((3 7) (both (expm a (/ (+ p 1) 4) p)))
  70.                   ((5) (let* ((x (expm a (/ (+ p 3) 8) p))
  71.                               (c (expm x 2 p)))
  72.                          (if (= a c) (both x)
  73.                            (both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p)))))
  74.                   (else (let* ((d (let loop ((d 2))
  75.                                     (if (= (jacobi d p) -1) d
  76.                                       (loop (+ d 1)))))
  77.                                (s (let loop ((p (- p 1)) (s 0))
  78.                                     (if (odd? p) s
  79.                                       (loop (quotient p 2) (+ s 1)))))
  80.                                (t (quotient (- p 1) (expt 2 s)))
  81.                                (big-a (expm a t p))
  82.                                (big-d (expm d t p))
  83.                                (m (let loop ((i 0) (m 0))
  84.                                     (cond ((= i s) m)
  85.                                           ((= (- p 1)
  86.                                               (expm (* big-a (expm big-d m p))
  87.                                                     (expt 2 (- s 1 i)) p))
  88.                                             (loop (+ i 1) (+ m (expt 2 i))))
  89.                                           (else (loop (+ i 1) m))))))
  90.                           (both (modulo (* (expm a (/ (+ t 1) 2) p)
  91.                                            (expm big-d (/ m 2) p)) p)))))))))
  92.  
  93. (define-syntax (with-modulus stx)
  94.   (syntax-case stx ()
  95.     ((with-modulus e expr ...)
  96.       (with-syntax ((modulus (datum->syntax-object (syntax with-modulus) 'modulus))
  97.                     (==   (datum->syntax-object (syntax with-modulus) '==   ))
  98.                     (+    (datum->syntax-object (syntax with-modulus) '+    ))
  99.                     (-    (datum->syntax-object (syntax with-modulus) '-    ))
  100.                     (*    (datum->syntax-object (syntax with-modulus) '*    ))
  101.                     (/    (datum->syntax-object (syntax with-modulus) '/    ))
  102.                     (^    (datum->syntax-object (syntax with-modulus) '^    ))
  103.                     (!    (datum->syntax-object (syntax with-modulus) '!    ))
  104.                     (sqrt (datum->syntax-object (syntax with-modulus) 'sqrt )))
  105.         (syntax (letrec ((fold (lambda (op base xs)
  106.                                  (if (null? xs) base
  107.                                    (fold op (op base (car xs)) (cdr xs))))))
  108.                   (let* ((modulus e)
  109.                          (mod (lambda (x)
  110.                                 (if (not (integer? x))
  111.                                     (error 'with-modulus "all arguments must be integers")
  112.                                     (modulo x modulus))))
  113.                          (== (lambda (x y) (= (mod x) (mod y))))
  114.                          (+ (lambda xs (fold (lambda (x y) (mod (+ x (mod y)))) 0 xs)))
  115.                          (- (lambda (x . xs)
  116.                               (if (null? xs)
  117.                                   (mod (- 0 x))
  118.                                   (fold (lambda (x y) (mod (- x (mod y)))) x xs))))
  119.                          (* (lambda xs (fold (lambda (x y) (mod (* x (mod y)))) 1 xs)))
  120.                          (/ (lambda (x . xs)
  121.                               (if (null? xs)
  122.                                   (inverse x e)
  123.                                   (fold (lambda (x y) (* x (inverse y e))) x xs))))
  124.                          (^ (lambda (base exp) (expm base exp e)))
  125.                          (! (lambda (n) (mod-fact n modulus)))
  126.                          (sqrt (lambda (x) (mod-sqrt x e))))
  127.                     expr ...)))))))
  128.  
  129. (define (twin? m)
  130.   (with-modulus (* m (+ m 2))
  131.     (== (* 4 (+ (! (- m 1)) 1))
  132.         (- m))))
  133.  
  134. (define (range . args)
  135.   (case (length args)
  136.     ((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
  137.     ((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
  138.     ((3) (let ((le? (if (negative? (caddr args)) >= <=)))
  139.            (let loop ((x(car args)) (xs '()))
  140.              (if (le? (cadr args) x)
  141.                  (reverse xs)
  142.                  (loop (+ x (caddr args)) (cons x xs))))))
  143.     (else (error 'range "unrecognized arguments"))))
  144.  
  145. (display (filter twin? (range 3 1000 2))) (newline)
Runtime error #stdin #stdout #stderr 0.03s 8716KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
ice-9/psyntax.scm:987:26: In procedure scan:
ice-9/psyntax.scm:987:26: Syntax error:
/home/086Bbo/prog.scm:93:0: source expression failed to match any pattern in form (define-syntax (with-modulus stx) (syntax-case stx () ((with-modulus e expr ...) (with-syntax ((modulus (datum->syntax-object (syntax with-modulus) (quote modulus))) (== (datum->syntax-object (syntax with-modulus) (quote ==))) (+ (datum->syntax-object (syntax with-modulus) (quote +))) (- (datum->syntax-object (syntax with-modulus) (quote -))) (* (datum->syntax-object (syntax with-modulus) (quote *))) (/ (datum->syntax-object (syntax with-modulus) (quote /))) (^ (datum->syntax-object (syntax with-modulus) (quote ^))) (! (datum->syntax-object (syntax with-modulus) (quote !))) (sqrt (datum->syntax-object (syntax with-modulus) (quote 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 (quote 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 ...)))))))
https://ideone.com/SZyqeb
language:
Scheme (guile 2.2.4)
created:
4 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!