[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/xSBfwIMn    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Oct 30:
; emirps

(define-syntax fold-of
  (syntax-rules (range in is)
    ((_ "z" f b e) (set! b (f b e)))
    ((_ "z" f b e (v range fst pst stp) c ...)
      (let* ((x fst) (p pst) (s stp)
             (le? (if (positive? s) <= >=)))
        (do ((v x (+ v s))) ((le? p v) b)
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (v range fst pst) c ...)
      (let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
        (fold-of "z" f b e (v range x p s) c ...)))
    ((_ "z" f b e (v range pst) c ...)
      (fold-of "z" f b e (v range 0 pst) c ...))
    ((_ "z" f b e (x in xs) c ...)
      (do ((t xs (cdr t))) ((null? t) b)
        (let ((x (car t)))
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (x is y) c ...)
      (let ((x y)) (fold-of "z" f b e c ...)))
    ((_ "z" f b e p? c ...)
      (if p? (fold-of "z" f b e c ...)))
    ((_ f i e c ...)
      (let ((b i)) (fold-of "z" f b e c ...)))))

(define-syntax list-of (syntax-rules ()
  ((_ arg ...) (reverse (fold-of
    (lambda (d a) (cons a d)) '() arg ...)))))

(define (digits n . args)
  (let ((b (if (null? args) 10 (car args))))
    (let loop ((n n) (d '()))
      (if (zero? n) d
          (loop (quotient n b)
                (cons (modulo n b) d))))))

(define (undigits ds . args)
  (let ((b (if (null? args) 10 (car args))))
    (let loop ((ds ds) (n 0))
      (if (null? ds) n
          (loop (cdr ds) (+ (* n b) (car ds)))))))

(define (prime? n)
  (define (expm b e m)
    (define (m* x y) (modulo (* x y) m))
    (cond ((zero? e) 1)
          ((even? e) (expm (m* b b) (/ e 2) m))
          (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
  (define (digits n . args)
    (let ((b (if (null? args) 10 (car args))))
      (let loop ((n n) (d '()))
        (if (zero? n) d
            (loop (quotient n b)
                  (cons (modulo n b) d))))))
  (define (isqrt n)
    (let loop ((x n) (y (quotient (+ n 1) 2)))
      (if (<= 0 (- y x) 1) x
        (loop y (quotient (+ y (quotient n y)) 2)))))
  (define (square? n)
    (let ((n2 (isqrt n)))
      (= n (* n2 n2))))
  (define (jacobi a n)
    (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
        (error 'jacobi "modulus must be positive odd integer")
        (let jacobi ((a a) (n n))
          (cond ((= a 0) 0)
                ((= a 1) 1)
                ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
                ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
                ((< n a) (jacobi (modulo a n) n))
                ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
                (else (jacobi n a))))))
  (define legendre jacobi)
  (define (inverse x n)
    (let loop ((x (modulo x n)) (a 1))
      (cond ((zero? x) (error 'inverse "division by zero"))
            ((= x 1) a)
            (else (let ((q (- (quotient n x))))
                    (loop (+ n (* q x)) (modulo (* q a) n)))))))
  (define (miller? n 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 (chain m f g x0 x1)
    (let loop ((ms (digits m 2)) (u x0) (v x1))
      (cond ((null? ms) (values u v))
            ((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))
            (else (loop (cdr ms) (g u v) (f v))))))
  (define (lucas? n)
    (let loop ((a 11) (b 7))
      (let ((d (- (* a a) (* 4 b))))
        (cond ((square? d) (loop (+ a 2) (+ b 1)))
              ((not (= (gcd n (* 2 a b d)) 1))
                (loop (+ a 2) (+ b 2)))
              (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))
                           (m (quotient (- n (legendre d n)) 2))
                           (f (lambda (u) (modulo (- (* u u) 2) n)))
                           (g (lambda (u v) (modulo (- (* u v) x1) n))))
                      (let-values (((xm xm1) (chain m f g 2 x1)))
                        (zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))
  (cond ((or (not (integer? n)) (< n 2))
          (error 'prime? "must be integer greater than one"))
        ((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))
        (else (and (miller? n 2) (miller? n 3) (lucas? n)))))

(define (emirp? n)
  (let ((rev (undigits (reverse (digits n)))))
    (and (prime? n) (prime? rev))))

(display (list-of n (n range 10 1000) (emirp? n)))


Output:
1
(11 13 17 31 37 71 73 79 97 101 107 113 131 149 151 157 167 179 181 191 199 311 313 337 347 353 359 373 383 389 701 709 727 733 739 743 751 757 761 769 787 797 907 919 929 937 941 953 967 971 983 991)


Create a new paste based on this one


Comments: