fork download
copy 
  1. ; perrin pseudoprimes
  2.  
  3. (define (make-matrix rows columns . value)
  4.   (do ((m (make-vector rows)) (i 0 (+ i 1)))
  5.       ((= i rows) m)
  6.     (if (null? value)
  7.         (vector-set! m i (make-vector columns))
  8.         (vector-set! m i (make-vector columns (car value))))))
  9.  
  10. (define (matrix-rows x) (vector-length x))
  11.  
  12. (define (matrix-cols x) (vector-length (vector-ref x 0)))
  13.  
  14. (define (matrix-ref m i j) (vector-ref (vector-ref m i) j))
  15.  
  16. (define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))
  17.  
  18. (define-syntax for
  19.   (syntax-rules ()
  20.     ((for (var first past step) body ...)
  21.       (let ((ge? (if (< first past) >= <=)))
  22.         (do ((var first (+ var step)))
  23.             ((ge? var past))
  24.           body ...)))
  25.     ((for (var first past) body ...)
  26.       (let* ((f first) (p past) (s (if (< first past) 1 -1)))
  27.         (for (var f p s) body ...)))
  28.     ((for (var past) body ...)
  29.       (let* ((p past)) (for (var 0 p) body ...)))))
  30.  
  31. (define (matrix-multiply a b)
  32.   (let ((ar (matrix-rows a)) (ac (matrix-cols a))
  33.         (br (matrix-rows b)) (bc (matrix-cols b)))
  34.     (if (not (= ac br))
  35.         (error 'matrix-multiply "incompatible matrices")
  36.         (let ((c (make-matrix ar bc 0)))
  37.           (for (i ar)
  38.             (for (j bc)
  39.               (for (k ac)
  40.                 (matrix-set! c i j
  41.                   (+ (matrix-ref c i j)
  42.                      (* (matrix-ref a i k)
  43.                         (matrix-ref b k j)))))))
  44.           c))))
  45.  
  46. (define (matrix-power m n)
  47.   (cond ((= n 1) m)
  48.         ((even? n)
  49.           (matrix-power
  50.             (matrix-multiply m m)
  51.             (/ n 2)))
  52.         (else (matrix-multiply m
  53.                 (matrix-power
  54.                   (matrix-multiply m m)
  55.                   (/ (- n 1) 2))))))
  56.  
  57. (define perrin-matrix '#(#(0 1 0) #(0 0 1) #(1 1 0)))
  58. (define perrin-init '#(#(3) #(0) #(2)))
  59.  
  60. (define (perrin n)
  61.   (if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2
  62.     (vector-ref (vector-ref
  63.       (matrix-multiply
  64.         (matrix-power perrin-matrix (- n 2))
  65.         perrin-init)
  66.       2) 0)))))
  67.  
  68. (display (perrin 200)) (newline)
  69.  
  70. (define (matrix-multiply-mod a b m)
  71.   (let ((ar (matrix-rows a)) (ac (matrix-cols a))
  72.         (br (matrix-rows b)) (bc (matrix-cols b)))
  73.     (if (not (= ac br))
  74.         (error 'matrix-multiply-mod "incompatible matrices")
  75.         (let ((c (make-matrix ar bc 0)))
  76.           (for (i ar)
  77.             (for (j bc)
  78.               (for (k ac)
  79.                 (matrix-set! c i j
  80.                   (modulo
  81.                     (+ (matrix-ref c i j)
  82.                        (* (matrix-ref a i k)
  83.                           (matrix-ref b k j))) m)))))
  84.           c))))
  85.  
  86. (define (matrix-power-mod b e m)
  87.   (cond ((= e 1) b)
  88.         ((even? e)
  89.           (matrix-power-mod
  90.             (matrix-multiply-mod b b m)
  91.             (/ e 2)
  92.             m))
  93.         (else
  94.           (matrix-multiply-mod
  95.             b
  96.             (matrix-power-mod
  97.               (matrix-multiply-mod b b m)
  98.               (/ (- e 1) 2)
  99.               m)
  100.             m))))
  101.  
  102. (define (perrin-mod n)
  103.   (if (= n 0) 3 (if (= n 1) 0 (if (= n 2) 2
  104.     (vector-ref (vector-ref
  105.       (matrix-multiply-mod
  106.         (matrix-power-mod perrin-matrix (- n 2) n)
  107.         perrin-init n)
  108.       2) 0)))))
  109.  
  110. (display (perrin-mod 200)) (newline)
  111. (display (modulo (perrin-mod 200) 200)) (newline)
  112.  
  113. (define (perrin-prime? n) (zero? (perrin-mod n)))
  114.  
  115. (display (perrin-prime? 271441)) (newline)
  116.  
  117. (define prime?
  118.   (let ((seed 3141592654))
  119.     (lambda (n)
  120.       (define (rand)
  121.         (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))
  122.         (+ (quotient (* seed (- n 2)) 4294967296) 2))
  123.       (define (expm b e m)
  124.         (define (times x y) (modulo (* x y) m))
  125.         (let loop ((b b) (e e) (r 1))
  126.           (if (zero? e) r
  127.             (loop (times b b) (quotient e 2)
  128.                   (if (odd? e) (times b r) r)))))
  129.       (define (spsp? n a)
  130.         (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
  131.             ((odd? d)
  132.               (let ((t (expm a d n)))
  133.                 (if (or (= t 1) (= t (- n 1))) #t
  134.                   (do ((s (- s 1) (- s 1))
  135.                        (t (expm t 2 n) (expm t 2 n)))
  136.                       ((or (zero? s) (= t (- n 1)))
  137.                         (positive? s))))))))
  138.       (if (not (integer? n))
  139.           (error 'prime? "must be integer")
  140.           (if (< n 2) #f
  141.             (do ((a (rand) (rand)) (k 25 (- k 1)))
  142.                 ((or (zero? k) (not (spsp? n a)))
  143.                   (zero? k))))))))
  144.  
  145. (display (prime? 271441)) (newline)
  146.  
  147. (define (perrin-pseudo-primes n)
  148.   (do ((i 2 (+ i 1))) ((< n i))
  149.     (when (and (perrin-prime? i)
  150.                (not (prime? i)))
  151.       (display i) (newline))))
  152.  
  153. ; about ten minutes -- linear
  154. ; > (perrin-pseudo-primes 1000000)
  155. ; 271441
  156. ; 904631
  157.  
  158. (define (perrin-pseudo-primes n)
  159.   (let loop ((p-2 3) (p-1 0) (p 2) (i 3))
  160.     (when (< i n)
  161.       (let ((p+1 (+ p-2 p-1)))
  162.         (when (and (zero? (modulo p+1 i))
  163.                    (not (prime? i)))
  164.           (display i) (newline))
  165.         (loop p-1 p p+1 (+ i 1))))))
  166.  
  167. ; about two minutes -- quadratic
  168. ; > (perrin-pseudo-primes 1000000)
  169. ; 271441
  170. ; 904631
Success #stdin #stdout 0.27s 9176KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
2658793989922287946990250
50
50
#t
#f
https://ideone.com/SOZgf8
language:
Scheme (guile 2.2.4)
created:
8 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!