fork download
copy 
  1. ; gauss's criterion
  2.  
  3. (define (range . args)
  4.   (case (length args)
  5.     ((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
  6.     ((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
  7.     ((3) (let ((le? (if (negative? (caddr args)) >= <=)))
  8.            (let loop ((x(car args)) (xs '()))
  9.              (if (le? (cadr args) x)
  10.                  (reverse xs)
  11.                  (loop (+ x (caddr args)) (cons x xs))))))
  12.     (else (error 'range "unrecognized arguments"))))
  13.  
  14. (define-syntax assert
  15.   (syntax-rules ()
  16.     ((assert expr result)
  17.       (if (not (equal? expr result))
  18.           (for-each display `(
  19.             #\newline "failed assertion:" #\newline
  20.             expr #\newline "expected: " ,result
  21.             #\newline "returned: " ,expr #\newline))))))
  22.  
  23. (define (gauss-criterion b p)
  24.   (define (gauss 2k-1) (modulo (* 2k-1 b) p))
  25.   (expt -1 (length (filter even? (filter positive? (map gauss (range 1 p 2)))))))
  26.  
  27. (define (jacobi a m)
  28.   (if (not (integer? a)) (error 'jacobi "must be integer")
  29.     (if (not (and (integer? m) (positive? m) (odd? m)))
  30.         (error 'jacobi "modulus must be odd positive integer")
  31.         (let loop1 ((a (modulo a m)) (m m) (t 1))
  32.           (if (zero? a) (if (= m 1) t 0)
  33.             (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
  34.               (let loop2 ((a a) (t t))
  35.                 (if (even? a) (loop2 (/ a 2) (* t z))
  36.                   (loop1 (modulo m a) a
  37.                          (if (and (= (modulo a 4) 3)
  38.                                   (= (modulo m 4) 3))
  39.                              (- t) t))))))))))
  40.  
  41. (define (primes n) ; list of primes not exceeding n
  42.   (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
  43.     (let loop ((i 0) (p 3) (ps (list 2)))
  44.       (cond ((< n (* p p))
  45.               (do ((i i (+ i 1)) (p p (+ p 2))
  46.                    (ps ps (if (vector-ref bits i) (cons p ps) ps)))
  47.                   ((= i len) (reverse ps))))
  48.             ((vector-ref bits i)
  49.               (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
  50.                   ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
  51.                 (vector-set! bits j #f)))
  52.             (else (loop (+ i 1) (+ p 2) ps))))))
  53.  
  54. (do ((ps (cdr (primes 100)) (cdr ps))) ((null? ps))
  55.   (do ((b 1 (+ b 1))) ((= (car ps) b))
  56.     (assert (gauss-criterion b (car ps))
  57.             (jacobi b (car ps))))); your code goes here
Success #stdin #stdout 0.58s 9728KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/ONFr06
language:
Scheme (guile 2.2.4)
created:
5 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!