fork download
copy 
  1. ; proth primes and sierpinski numbers
  2.  
  3. (define (expm b e m)
  4.   (define (m* x y) (modulo (* x y) m))
  5.   (cond ((zero? e) 1)
  6.         ((even? e) (expm (m* b b) (/ e 2) m))
  7.         (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
  8.  
  9. (define (jacobi a m)
  10.   (if (not (integer? a)) (error 'jacobi "must be integer")
  11.     (if (not (and (integer? m) (positive? m) (odd? m)))
  12.         (error 'jacobi "modulus must be odd positive integer")
  13.         (let loop1 ((a (modulo a m)) (m m) (t 1))
  14.           (if (zero? a) (if (= m 1) t 0)
  15.             (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
  16.               (let loop2 ((a a) (t t))
  17.                 (if (even? a) (loop2 (/ a 2) (* t z))
  18.                   (loop1 (modulo m a) a
  19.                          (if (and (= (modulo a 4) 3)
  20.                                   (= (modulo m 4) 3))
  21.                              (- t) t))))))))))
  22.  
  23. (define (proth? n)
  24.   (let loop ((a '(3 5 7 17)))
  25.     (if (null? a) #f
  26.       (if (not (= (jacobi (car a) n) -1)) (loop (cdr a))
  27.         (if (= (expm (car a) (/ (- n 1) 2) n) (- n 1)) #t
  28.           (loop (cdr a)))))))
  29.  
  30. (define (proth k limit)
  31.   (do ((k2n (* k 2) (* k2n 2)) (n 1 (+ n 1))) ((< limit n))
  32.     (when (proth? (+ k2n 1)) (display n) (newline))))
  33.  
  34. (proth 12909 1000)
Success #stdin #stdout 14.87s 10424KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1
2
5
7
11
14
17
18
22
30
39
58
62
77
85
91
99
109
125
131
149
171
185
215
223
253
267
285
287
299
310
337
347
391
639
641
685
767
771
781
887
https://ideone.com/e6FnQS
language:
Scheme (guile 2.2.4)
created:
6 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!