[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Feb 16:
; home primes

(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 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 (prime? n)
  (letrec (
    (expm (lambda (b e m)
      (let ((times (lambda (x y) (modulo (* x y) m))))
        (cond ((zero? e) 1) ((even? e) (expm (times b b) (/ e 2) m))
        (else (times b (expm (times b b) (/ (- e 1) 2) m)))))))
    (digits (lambda (n)
      (let loop ((n n) (ds '()))
        (if (zero? n) ds (loop (quotient n 2) (cons (modulo n 2) ds))))))
    (isqrt (lambda (n)
      (let loop ((x n) (y (quotient (+ n 1) 2)))
        (if (<= 0 (- y x) 1) x (loop y (quotient (+ y (quotient n y)) 2))))))
    (square? (lambda (n) (let ((n2 (isqrt n))) (= n (* n2 n2)))))
    (jacobi (lambda (a n)
      (let loop ((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) (* (loop 2 n) (loop (/ a 2) n)))
              ((< n a) (loop (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (loop n a)))
              (else (loop n a))))))
    (inverse (lambda (x m)
      (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w m))
        (if (zero? w) (modulo a m)
          (let ((q (quotient g w)))
            (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w))))))))
    (strong-pseudo-prime? (lambda (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))))))))))
    (chain (lambda (m f g u v)
      (let loop ((ms (digits m)) (u u) (v v))
        (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)))))))
    (lucas-pseudo-prime? (lambda (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 (jacobi 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)))))))))))
    (if (not (integer? n)) (error 'prime? "must be integer")
      (if (< n 2) #f (if (even? n) (= n 2) (if (zero? (modulo n 3)) (= n 3)
        (and (strong-pseudo-prime? n 2)
             (strong-pseudo-prime? n 3)
             (lucas-pseudo-prime? n))))))))

(define (factors n)
  (letrec (
    (last-pair (lambda (xs) (if (pair? (cdr xs)) (last-pair (cdr xs)) xs)))
    (cycle (lambda xs (set-cdr! (last-pair xs) xs) xs))
    (cons* (lambda (x . xs)
      (if (null? xs) x (cons x (apply cons* (car xs) (cdr xs)))))))
    (if (not (integer? n)) (error 'factors "must be integer")
      (if (member n '(-1 0 1)) (cdr (assoc n '((-1 -1) (0 0) (1 1))))
        (if (negative? n) (cons -1 (factors (- n)))
          (let wheel ((n n) (p 2) (ps '()) (ws (cons* 1 2 2 (cycle 4 2 4 2 4 6 2 6))))
            (if (= n 1) (reverse ps)
              (if (< n (* p p)) (reverse (cons n ps))
                (if (zero? (modulo n p)) (wheel (/ n p) p (cons p ps) ws)
                  (if (< p 1000) (wheel n (+ p (car ws)) ps (cdr ws))
                    (let rho ((ps ps) (cs (list n)) (a 1) (m 10000000))
                      (letrec (
                        (f (lambda (y) (modulo (+ (* y y) a) (car cs))))
                        (g (lambda (p x y) (modulo (* p (abs (- x y))) (car cs)))))
                        (if (null? cs) (sort < ps)
                          (let ((n (car cs)))
                            (if (= n 1) (rho ps (cdr cs) a m)
                              (if (prime? n) (rho (cons n ps) (cdr cs) a m)
                                (let loop ((x 2) (y (f 2)) (j 1) (q 2) (p 1))
                                  (if (= j m) (cons (sort < ps) cs)
                                    (if (= j q)
                                        (let ((t (f y)))
                                          (loop y t (+ j 1) (* q 2) (g p y t)))
                                        (if (positive? (modulo j 100))
                                            (loop x (f y) (+ j 1) q (g p x y))
                                            (let ((d (gcd p n)))
                                              (if (= d 1) (loop x (f y) (+ j 1) q 1)
                                                (if (= d n) (rho ps cs (+ a 1) (- m j))
                                                  (rho ps (cons* d (/ n d) (cdr cs))
                                                    a (- m j)))))))))))))))))))))))))

(define (home-prime n)
  (let* ((fs (factors n))
         (hp (undigits (apply append (map digits fs)))))
    (if (prime? hp) hp (home-prime hp))))

(define (home-factors n)
  (let loop ((n n) (hps '()))
    (if (prime? n) (reverse (cons n hps))
      (let* ((fs (factors n))
             (hp (undigits (apply append (map digits fs)))))
        (loop hp (cons fs hps))))))

(display (home-prime 99)) (newline)
(display (home-factors 8)) (newline)

(define (td-least n)
  (let loop ((p 3))
    (cond ((< 1000000 p) #f)
          ((zero? (modulo n p)) p)
          (else (loop (+ p 2))))))

(define (euclid-mullin)
  (let loop ((n 1) (ems '(2)))
    (display n) (display " ")
    (display (car ems)) (newline)
    (let* ((next (+ (apply * ems) 1))
           (lpf (td-least next)))
      (if lpf (loop (+ n 1) (cons lpf ems))
        (let ((fs (factors next)))
          (cond ((= n 16)
                  (loop 17 (cons 30693651606209 ems)))
                ((= n 27)
                  (loop 28 (cons 643679794963466223081509857 ems)))
                ((number? (car fs))
                  (loop (+ n 1) (cons (car fs) ems)))
                ((pair? (car fs))
                  (loop (+ n 1) (cons (caar fs) ems)))
                (else (error 'euclid-mullin
                        "can't complete factorization")))))))) 

(euclid-mullin)


Output:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
71143
((2 2 2) (2 3 37) (3 19 41) (3 3 3 7 13 13) (3 11123771) (7 149 317 941) (229 31219729) (11 2084656339) (3 347 911 118189) (11 613 496501723) (97 130517 917327) (53 1832651281459) (3 3 3 11 139 653 3863 5107) 3331113965338635107)
1 2
2 3
3 7
4 43
5 13
6 53
7 5
8 6221671
9 38709183810571
10 139
11 2801
12 11
13 17
14 5471
15 52662739
16 23003

Timeout


Create a new paste based on this one


Comments: