Question

I am trying to write a Scheme macro that loops over the prime numbers. Here is a simple version of the macro:

(define-syntax do-primes
  (syntax-rules ()
    ((do-primes (p lo hi) (binding ...) (test res ...) exp ...)
      (do ((p (next-prime (- lo 1)) (next-prime p)) binding ...)
          ((or test (< hi p)) res ...) exp ...))
    ((do-primes (p lo) (binding ...) (test res ...) exp ...)
      (do ((p (next-prime (- lo 1)) (next-prime p)) binding ...)
          (test res ...) exp ...))
    ((do-primes (p) (binding ...) (test res ...) exp ...)
      (do ((p 2 (next-prime p)) binding ...) (test res ...) exp ...))))

The do-primes macro expands to a do with three possible syntaxes: if the first argument to do-primes is (p lo hi) then the do loops over the primes from lo to hi unless iteration is stopped early by the termination clause, if the first argument to do-primes is (p lo) then the do loops over the primes starting from lo and continuing until the termination clause stops the iteration, and if the first argument to do-primes is (p) then the do loops over the primes starting from 2 and continuing until the termination clause stops the iteration. Here are some sample uses of the do-primes macro:

; these lines display the primes less than 25
(do-primes (p 2 25) () (#f) (display p) (newline))
(do-primes (p 2) () ((< 25 p)) (display p) (newline))
(do-primes (p) () ((< 25 p)) (display p) (newline))

; these lines return a list of the primes less than 25
(do-primes (p 2 25) ((ps (list) (cons p ps))) (#f (reverse ps)))
(do-primes (p 2) ((ps (list) (cons p ps))) ((< 25 p) (reverse ps)))
(do-primes (p) ((ps (list) (cons p ps))) ((< 25 p) (reverse ps)))

; these lines return the sum of the primes less than 25
(do-primes (p 2 25) ((s 0 (+ s p))) (#f s))
(do-primes (p 2) ((s 0 (+ s p))) ((< 25 p) s))
(do-primes (p) ((s 0 (+ s p))) ((< 25 p) s))

What I want to do is write a version of the do-primes macro that uses a local version of the next-prime function; I want to do that because I can make the next-prime function run faster than my general-purpose next-prime function because I know the environment in which it will be called. I tried to write the macro like this:

(define-syntax do-primes
  (let ()
    (define (prime? n)
      (if (< n 2) #f
        (let loop ((f 2))
          (if (< n (* f f)) #t
            (if (zero? (modulo n f)) #f
              (loop (+ f 1)))))))
    (define (next-prime n)
      (let loop ((n (+ n 1)))
        (if (prime? n) n (loop (+ n 1)))))
  (lambda (x) (syntax-case x ()
    ((do-primes (p lo hi) (binding ...) (test res ...) exp ...)
      (syntax
        (do ((p (next-prime (- lo 1)) (next-prime p)) binding ...)
            ((or test (< hi p)) res ...) exp ...)))
    ((do-primes (p lo) (binding ...) (test res ...) exp ...)
      (syntax
        (do ((p (next-prime (- lo 1)) (next-prime p)) binding ...)
            (test res ...) exp ...)))
    ((do-primes (p) (binding ...) (test res ...) exp ...)
      (syntax
        (do ((p 2 (next-prime p)) binding ...) (test res ...) exp ...))))))))

(Ignore the prime? and next-prime functions, which are there only for illustration. The real version of the do-primes macro will use a segmented sieve for small primes and switch to a Baillie-Wagstaff pseudoprime test for larger primes.) But that doesn't work; I get an error message telling me that that I am trying "to reference out-of-phase identifier next-prime." I understand the problem. But my macrology wizardry is insufficient to solve it.

Can someone show me how to write the do-primes macro?

EDIT: Here is the final macro:

(define-syntax do-primes (syntax-rules () ; syntax for iterating over primes

  ; (do-primes (p lo hi) ((var init next) ...) (pred? result ...) expr ...)

  ; Macro do-primes provides syntax for iterating over primes. It expands to
  ; a do-loop with variable p bound in the same scope as the rest of the (var
  ; init next) variables, as if it were defined as (do ((p (primes lo hi) (cdr
  ; p)) (var init next) ...) (pred result ...) expr ...). Variables lo and hi
  ; are inclusive; for instance, given (p 2 11), p will take on the values 2,
  ; 3, 5, 7 and 11. If hi is omitted the iteration continues until stopped by
  ; pred?. If lo is also omitted iteration starts from 2. Some examples:

  ;   three ways to display the primes less than 25
  ;   (do-primes (p 2 25) () (#f) (display p) (newline))
  ;   (do-primes (p 2) () ((< 25 p)) (display p) (newline))
  ;   (do-primes (p) () ((< 25 p)) (display p) (newline))

  ;   three ways to return a list of the primes less than 25
  ;   (do-primes (p 2 25) ((ps (list) (cons p ps))) (#f (reverse ps)))
  ;   (do-primes (p 2) ((ps (list) (cons p ps))) ((< 25 p) (reverse ps)))
  ;   (do-primes (p) ((ps (list) (cons p ps))) ((< 25 p) (reverse ps)))

  ;   three ways to return the sum of the primes less than 25
  ;   (do-primes (p 2 25) ((s 0 (+ s p))) (#f s))
  ;   (do-primes (p 2) ((s 0 (+ s p))) ((< 25 p) s))
  ;   (do-primes (p) ((s 0 (+ s p))) ((< 25 p) s))

  ;   functions to count primes and return the nth prime (from P[1] = 2)
  ;   (define (prime-pi n) (do-primes (p) ((k 0 (+ k 1))) ((< n p) k)))
  ;   (define (nth-prime n) (do-primes (p) ((n n (- n 1))) ((= n 1) p)))

  ; The algorithm used to generate primes is a segmented Sieve of Eratosthenes
  ; up to 2^32. For larger primes, a segmented sieve runs over the sieving
  ; primes up to 2^16 to produce prime candidates, then a Baillie-Wagstaff
  ; pseudoprimality test is performed to confirm the number is prime.

  ; If functions primes, expm, jacobi, strong-pseudoprime?, lucas, selfridge
  ; and lucas-pseudoprime? exist in the outer environment, they can be removed
  ; from the macro.

  ((do-primes (p lo hi) (binding ...) (test result ...) expr ...)
    (do-primes (p lo) (binding ...) ((or test (< hi p)) result ...) expr ...))

  ((do-primes (pp low) (binding ...) (test result ...) expr ...)

    (let* ((limit (expt 2 16)) (delta 50000) (limit2 (* limit limit))
           (sieve (make-vector delta #t)) (ps #f) (qs #f) (bottom 0) (pos 0))

      (define (primes n) ; sieve of eratosthenes
        (let ((sieve (make-vector n #t)))
          (let loop ((p 2) (ps (list)))
            (cond ((= n p) (reverse ps))
                  ((vector-ref sieve p)
                    (do ((i (* p p) (+ i p))) ((<= n i))
                      (vector-set! sieve i #f))
                    (loop (+ p 1) (cons p ps)))
                  (else (loop (+ p 1) ps))))))

      (define (expm b e m) ; modular exponentiation
        (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 (jacobi a m) ; jacobi symbol
        (let loop1 ((a (modulo a m)) (m m) (t 1))
          (if (zero? a) (if (= m 1) t 0)
            (let ((z (if (member (modulo m 8) (list 3 5)) -1 1)))
              (let loop2 ((a a) (t t))
                (if (even? a) (loop2 (/ a 2) (* t z))
                  (loop1 (modulo m a) a
                         (if (and (= (modulo a 4) 3)
                                  (= (modulo m 4) 3))
                             (- t) t))))))))

      (define (strong-pseudoprime? n a) ; strong pseudoprime base a
        (let loop ((r 0) (s (- n 1)))
          (if (even? s) (loop (+ r 1) (/ s 2))
            (if (= (expm a s n) 1) #t
              (let loop ((r r) (s s))
                (cond ((zero? r) #f)
                      ((= (expm a s n) (- n 1)) #t)
                      (else (loop (- r 1) (* s 2)))))))))

      (define (lucas p q m n) ; lucas sequences u[n] and v[n] and q^n (mod m)
        (define (even e o) (if (even? n) e o))
        (define (mod n) (if (zero? m) n (modulo n m)))
        (let ((d (- (* p p) (* 4 q))))
          (let loop ((un 1) (vn p) (qn q) (n (quotient n 2))
                     (u (even 0 1)) (v (even 2 p)) (k (even 1 q)))
            (if (zero? n) (values u v k)
              (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn))))
                    (q2 (mod (* qn qn))) (n2 (quotient n 2)))
                (if (even? n) (loop u2 v2 q2 n2 u v k)
                  (let* ((uu (+ (* u v2) (* u2 v)))
                         (vv (+ (* v v2) (* d u u2)))
                         (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu))
                         (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv))
                         (uu (mod (/ uu 2))) (vv (mod (/ vv 2))))
                    (loop u2 v2 q2 n2 uu vv (* k q2)))))))))

      (define (selfridge n) ; initialize lucas sequence
        (let loop ((d-abs 5) (sign 1))
          (let ((d (* d-abs sign)))
            (cond ((< 1 (gcd d n)) (values d 0 0))
                  ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4)))
                  (else (loop (+ d-abs 2) (- sign)))))))

      (define (lucas-pseudoprime? n) ; standard lucas pseudoprime
        (call-with-values
          (lambda () (selfridge n))
          (lambda (d p q)
            (if (zero? p) (= n d)
              (call-with-values
                (lambda () (lucas p q n (+ n 1)))
                (lambda (u v qkd) (zero? u)))))))

      (define (init lo) ; initialize sieve, return first prime
        (set! bottom (if (< lo 3) 2 (if (odd? lo) (- lo 1) lo)))
        (set! ps (cdr (primes limit))) (set! pos 0)
        (set! qs (map (lambda (p) (modulo (/ (+ bottom p 1) -2) p)) ps))
        (do ((p ps (cdr p)) (q qs (cdr q))) ((null? p))
          (do ((i (+ (car p) (car q)) (+ i (car p)))) ((<= delta i))
            (vector-set! sieve i #f)))
        (if (< lo 3) 2 (next)))

      (define (advance) ; advance to next segment
        (set! bottom (+ bottom delta delta)) (set! pos 0)
        (do ((i 0 (+ i 1))) ((= i delta)) (vector-set! sieve i #t))
        (set! qs (map (lambda (p q) (modulo (- q delta) p)) ps qs))
        (do ((p ps (cdr p)) (q qs (cdr q))) ((null? p))
          (do ((i (car q) (+ i (car p)))) ((<= delta i))
            (vector-set! sieve i #f))))

      (define (next) ; next prime after current prime
        (when (= pos delta) (advance))
        (let ((p (+ bottom pos pos 1)))
          (if (and (vector-ref sieve pos) (or (< p limit2)
                (and (strong-pseudoprime? p 2) (lucas-pseudoprime? p))))
              (begin (set! pos (+ pos 1)) p)
              (begin (set! pos (+ pos 1)) (next)))))

      (do ((pp (init low) (next)) binding ...) (test result ...) expr ...)))

  ((do-primes (p) (binding ...) (test result ...) expr ...)
    (do-primes (p 2) (binding ...) (test result ...) expr ...))))
Was it helpful?

Solution

To get correct phasing, your next-prime needs to defined within the macro output. Here's one way to go about it (tested with Racket):

(define-syntax do-primes
  (syntax-rules ()
    ((do-primes (p lo hi) (binding ...) (test res ...) exp ...)
     (do-primes (p lo) (binding ...) ((or test (< hi p)) res ...) exp ...))
    ((do-primes (p lo) (binding ...) (test res ...) exp ...)
     (let ()
       (define (prime? n)
         ...)
       (define (next-prime n)
         ...)
       (do ((p (next-prime (- lo 1)) (next-prime p)) binding ...)
           (test res ...)
         exp ...)))
    ((do-primes (p) (binding ...) (test res ...) exp ...)
     (do-primes (p 2) (binding ...) (test res ...) exp ...))))

This way, this defines the prime? and next-prime in the most local scope possible, while not having tons of duplicate code in your macro definition (since the 1- and 3-argument forms are simply rewritten to use the 2-argument form).

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top