[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Oct 31:
; emirps

(define (take-while pred? xs)
  (let loop ((xs xs) (ys '()))
    (if (or (null? xs) (not (pred? (car xs))))
        (reverse ys)
        (loop (cdr xs) (cons (car xs) ys)))))

(define (filter pred? xs)
  (let loop ((xs xs) (ys '()))
    (cond ((null? xs) (reverse ys))
          ((pred? (car xs))
            (loop (cdr xs) (cons (car xs) ys)))
          (else (loop (cdr xs) ys)))))

(define (range . args)
  (case (length args)
    ((1) (range 0 (car args) (if (negative? (car args)) -1 1)))
    ((2) (range (car args) (cadr args) (if (< (car args) (cadr args)) 1 -1)))
    ((3) (let ((le? (if (negative? (caddr args)) >= <=)))
           (let loop ((x(car args)) (xs '()))
             (if (le? (cadr args) x)
                 (reverse xs)
                 (loop (+ x (caddr args)) (cons x xs))))))
    (else (error 'range "unrecognized arguments"))))

(define sort #f)
(define merge #f)
(let ()
  (define dosort
    (lambda (pred? ls n)
      (if (= n 1)
          (list (car ls))
          (let ((i (quotient n 2)))
            (domerge pred?
                     (dosort pred? ls i)
                     (dosort pred? (list-tail ls i) (- n i)))))))
  (define domerge
    (lambda (pred? l1 l2)
      (cond
        ((null? l1) l2)
        ((null? l2) l1)
        ((pred? (car l2) (car l1))
         (cons (car l2) (domerge pred? l1 (cdr l2))))
        (else (cons (car l1) (domerge pred? (cdr l1) l2))))))
  (set! sort
    (lambda (pred? l)
      (if (null? l) l (dosort pred? l (length l)))))
  (set! merge
    (lambda (pred? l1 l2)
      (domerge pred? l1 l2))))

(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 (log10 x) (/ (log x) (log 10)))

(define (complement f) (lambda xs (not (apply f xs))))

(define (right-section proc . args)
  (lambda xs (apply proc (reverse (append (reverse args) (reverse xs))))))

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse ps))
                        ((vector-ref v j)
                          (loop (+ j 1) (cons (+ j j 3) ps)))
                        (else (loop (+ j 1) ps)))))
              ((vector-ref v i)
                (let loop ((j startj))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

(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 (rev n) (undigits (reverse (digits n))))
(define (palin? n) (= n (rev n)))
(define (emirp? n) (and (prime? n) (prime? (rev n)) (not (palin? n))))
(define (emirps1 n) (filter emirp? (range 2 n)))

(define (next-power-of-ten n)
  (let loop ((i 1) (ten^i 10))
    (if (<= n ten^i) (expt 10 i)
      (loop (+ i 1) (* ten^i 10)))))

(define (emirps2 n)
  (let ((n10 (next-power-of-ten n)))
    (let loop ((ps (filter (complement palin?) (primes n10))) (es '()))
      (cond ((null? ps) (take-while (right-section < n) (sort < es)))
            ((member (rev (car ps)) (cdr ps))
              (loop (cdr ps) (cons (car ps) (cons (rev (car ps)) es))))
            (else (loop (cdr ps) es))))))

(define (emirps3 n)
  (let* ((n10 (next-power-of-ten n))
         (ps (filter (complement palin?) (primes n10)))
         (qs (sort < (append ps (map rev ps)))))
    (let loop ((qs qs) (es '()))
      (cond ((or (< n (car qs)) (null? (cdr qs))) (reverse es))
            ((= (car qs) (cadr qs))
              (loop (cddr qs) (cons (car qs) es)))
            (else (loop (cdr qs) es))))))

(display (emirps1 500)) (newline)
(display (emirps2 500)) (newline)
(display (emirps3 500)) (newline)


Output:
1
2
3
(13 17 31 37 71 73 79 97 107 113 149 157 167 179 199 311 337 347 359 389)
(13 17 31 37 71 73 79 97 107 113 149 157 167 179 199 311 337 347 359 389)
(13 17 31 37 71 73 79 97 107 113 149 157 167 179 199 311 337 347 359 389)


Create a new paste based on this one


Comments: