[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Jan 20:
; pollard rho, revisited

(define (cons* first . rest)
  (let loop ((curr first) (rest rest))
    (if (null? rest) curr
        (cons curr (loop (car rest) (cdr rest))))))

(define (last-pair xs)
    (if (pair? (cdr xs))
        (last-pair (cdr xs))
        xs))

(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)
             (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))))))
           (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? "non-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 (wheel-factors n limit)
  (define (wheel . xs) (set-cdr! (last-pair xs) xs) xs)
  (let loop ((n n) (f 2) (fs '())
             (ws (cons* 1 2 2 (wheel 4 2 4 2 4 6 2 6))))
    (cond ((< limit f) (values n (reverse fs)))
          ((or (= n 1) (< n (* f f))) (values 1 (reverse (cons n fs))))
          ((zero? (modulo n f)) (loop (/ n f) f (cons f fs) ws))
          (else (loop n (+ f (car ws)) fs (cdr ws))))))

(define (pollard factor n)
  (let-values (((n fs) (wheel-factors n 1000)))
    (sort < (let fact ((n n) (fs fs))
      (cond ((= n 1) fs)
            ((prime? n) (cons n fs))
            (else (let ((f (factor n 1 100000000)))
                  (append fs (fact f '()) (fact (/ n f) '())))))))))

(define (floyd n c limit)
  (define (f x) (modulo (+ (* x x) c) n))
  (define (g p t h) (modulo (* p (abs (- t h))) n))
  (let loop1 ((j 1) (t 2) (h (f 2)) (x 2) (y (f 2)) (p 1))
    (if (= j limit) (error 'floyd "timeout")
      (if (= t h) (floyd n (+ c 1) (- limit j))
        (if (positive? (modulo j 100)) (loop1 (+ j 1) (f t) (f (f h)) x y (g p t h))
          (let ((d (gcd p n)))
            (if (= d 1) (loop1 (+ j 1) (f t) (f (f h)) t h 1)
              (if (< 1 d n) d
                (let loop2 ((k 1) (x x) (y y) (d (gcd (- x y) n)))
                  (if (= k 100) (floyd n (+ c 1) (- limit j))
                    (if (= d 1) (loop2 (+ k 1) (f x) (f (f y)) (gcd (- x y) n))
                      (if (= d n) (floyd n (+ c 1) (- limit j))
                        d))))))))))))

(define (brent n c limit)
  (define (f y) (modulo (+ (* y y) c) n))
  (define (g p x y) (modulo (* p (abs (- x y))) n))
  (let loop1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1))
    (if (= j limit) (error 'brent "timeout")
      (if (= x y) (brent n (+ c 1) (- limit j)) ; cycle
        (if (= j q) (let ((t (f y))) (loop1 y (f y) z (+ j 1) (* q 2) (g p y t)))
          (if (positive? (modulo j 100)) (loop1 x (f y) z (+ j 1) q (g p x y))
            (let ((d (gcd p n)))
              (if (= d 1) (loop1 x (f y) y (+ j 1) q (g p x y))
                (if (< 1 d n) d ; factor
                  (let loop2 ((k 1) (z (f z)))
                    (if (= k 100) (brent n (+ c 1) (- limit j))
                      (let ((d (gcd (- z x) n)))
                        (if (= d 1) (loop2 (+ k 1) (f z))
                          (if (= d n) (brent n (+ c 1) (- limit j))
                            d))))))))))))))

(define (fib n)
  (define (square x) (* x x))
  (cond ((= n 0) 0) ((= n 1) 1) ((= n 2) 1)
        ((odd? n) (let ((n2 (+ (quotient n 2) 1)))
                    (+ (square (fib n2)) (square (fib (- n2 1))))))
        (else (let ((n2 (quotient n 2)))
                (* (fib n2) (+ (* (fib (- n2 1)) 2) (fib n2)))))))

(time (do ((i 3 (+ i 1))) ((< 100 i))
  (display i) (display " ")
  (display (pollard floyd (fib i)))
  (newline)))

(time (do ((i 3 (+ i 1))) ((< 100 i))
  (display i) (display " ")
  (display (pollard brent (fib i)))
  (newline)))


Output:
3 (2)
4 (3)
5 (5)
6 (2 2 2)
7 (13)
8 (3 7)
9 (2 17)
10 (5 11)
11 (89)
12 (2 2 2 2 3 3)
13 (233)
14 (13 29)
15 (2 5 61)
16 (3 7 47)
17 (1597)
18 (2 2 2 17 19)
19 (37 113)
20 (3 5 11 41)
21 (2 13 421)
22 (89 199)
23 (28657)
24 (2 2 2 2 2 3 3 7 23)
25 (5 5 3001)
26 (233 521)
27 (2 17 53 109)
28 (3 13 29 281)
29 (514229)
30 (2 2 2 5 11 31 61)
31 (557 2417)
32 (3 7 47 2207)
33 (2 89 19801)
34 (1597 3571)
35 (5 13 141961)
36 (2 2 2 2 3 3 3 17 19 107)
37 (73 149 2221)
38 (37 113 9349)
39 (2 233 135721)
40 (3 5 7 11 41 2161)
41 (2789 59369)
42 (2 2 2 13 29 211 421)
43 (433494437)
44 (3 43 89 199 307)
45 (2 5 17 61 109441)
46 (139 461 28657)
47 (2971215073)
48 (2 2 2 2 2 2 3 3 7 23 47 1103)
49 (13 97 6168709)
50 (5 5 11 101 151 3001)
51 (2 1597 6376021)
52 (3 233 521 90481)
53 (953 55945741)
54 (2 2 2 17 19 53 109 5779)
55 (5 89 661 474541)
56 (3 7 7 13 29 281 14503)
57 (2 37 113 797 54833)
58 (59 19489 514229)
59 (353 2710260697)
60 (2 2 2 2 3 3 5 11 31 41 61 2521)
61 (4513 555003497)
62 (557 2417 3010349)
63 (2 13 17 421 35239681)
64 (3 7 47 1087 2207 4481)
65 (5 233 14736206161)
66 (2 2 2 89 199 9901 19801)
67 (269 116849 1429913)
68 (3 67 1597 3571 63443)
69 (2 137 829 18077 28657)
70 (5 11 13 29 71 911 141961)
71 (6673 46165371073)
72 (2 2 2 2 2 3 3 3 7 17 19 23 107 103681)
73 (9375829 86020717)
74 (73 149 2221 54018521)
75 (2 5 5 61 3001 230686501)
76 (3 37 113 9349 29134601)
77 (13 89 988681 4832521)
78 (2 2 2 79 233 521 859 135721)
79 (157 92180471494753)
80 (3 5 7 11 41 47 1601 2161 3041)
81 (2 17 53 109 2269 4373 19441)
82 (2789 59369 370248451)
83 (99194853094755497)
84 (2 2 2 2 3 3 13 29 83 211 281 421 1427)
85 (5 1597 9521 3415914041)
86 (6709 144481 433494437)
87 (2 173 514229 3821263937)
88 (3 7 43 89 199 263 307 881 967)
89 (1069 1665088321800481)
90 (2 2 2 5 11 17 19 31 61 181 541 109441)
91 (13 13 233 741469 159607993)
92 (3 139 461 4969 28657 275449)
93 (2 557 2417 4531100550901)
94 (2971215073 6643838879)
95 (5 37 113 761 29641 67735001)
96 (2 2 2 2 2 2 2 3 3 7 23 47 769 1103 2207 3167)
97 (193 389 3084989 361040209)
98 (13 29 97 6168709 599786069)
99 (2 17 89 197 19801 18546805133)
100 (3 5 5 11 41 101 151 401 3001 570601)
cpu time: 280 real time: 1413 gc time: 110
3 (2)
4 (3)
5 (5)
6 (2 2 2)
7 (13)
8 (3 7)
9 (2 17)
10 (5 11)
11 (89)
12 (2 2 2 2 3 3)
13 (233)
14 (13 29)
15 (2 5 61)
16 (3 7 47)
17 (1597)
18 (2 2 2 17 19)
19 (37 113)
20 (3 5 11 41)
21 (2 13 421)
22 (89 199)
23 (28657)
24 (2 2 2 2 2 3 3 7 23)
25 (5 5 3001)
26 (233 521)
27 (2 17 53 109)
28 (3 13 29 281)
29 (514229)
30 (2 2 2 5 11 31 61)
31 (557 2417)
32 (3 7 47 2207)
33 (2 89 19801)
34 (1597 3571)
35 (5 13 141961)
36 (2 2 2 2 3 3 3 17 19 107)
37 (73 149 2221)
38 (37 113 9349)
39 (2 233 135721)
40 (3 5 7 11 41 2161)
41 (2789 59369)
42 (2 2 2 13 29 211 421)
43 (433494437)
44 (3 43 89 199 307)
45 (2 5 17 61 109441)
46 (139 461 28657)
47 (2971215073)
48 (2 2 2 2 2 2 3 3 7 23 47 1103)
49 (13 97 6168709)
50 (5 5 11 101 151 3001)
51 (2 1597 6376021)
52 (3 233 521 90481)
53 (953 55945741)
54 (2 2 2 17 19 53 109 5779)
55 (5 89 661 474541)
56 (3 7 7 13 29 281 14503)
57 (2 37 113 797 54833)
58 (59 19489 514229)
59 (353 2710260697)
60 (2 2 2 2 3 3 5 11 31 41 61 2521)
61 (4513 555003497)
62 (557 2417 3010349)
63 (2 13 17 421 35239681)
64 (3 7 47 1087 2207 4481)
65 (5 233 14736206161)
66 (2 2 2 89 199 9901 19801)
67 (269 116849 1429913)
68 (3 67 1597 3571 63443)
69 (2 137 829 18077 28657)
70 (5 11 13 29 71 911 141961)
71 (6673 46165371073)
72 (2 2 2 2 2 3 3 3 7 17 19 23 107 103681)
73 (9375829 86020717)
74 (73 149 2221 54018521)
75 (2 5 5 61 3001 230686501)
76 (3 37 113 9349 29134601)
77 (13 89 988681 4832521)
78 (2 2 2 79 233 521 859 135721)
79 (157 92180471494753)
80 (3 5 7 11 41 47 1601 2161 3041)
81 (2 17 53 109 2269 4373 19441)
82 (2789 59369 370248451)
83 (99194853094755497)
84 (2 2 2 2 3 3 13 29 83 211 281 421 1427)
85 (5 1597 9521 3415914041)
86 (6709 144481 433494437)
87 (2 173 514229 3821263937)
88 (3 7 43 89 199 263 307 881 967)
89 (1069 1665088321800481)
90 (2 2 2 5 11 17 19 31 61 181 541 109441)
91 (13 13 233 741469 159607993)
92 (3 139 461 4969 28657 275449)
93 (2 557 2417 4531100550901)
94 (2971215073 6643838879)
95 (5 37 113 761 29641 67735001)
96 (2 2 2 2 2 2 2 3 3 7 23 47 769 1103 2207 3167)
97 (193 389 3084989 361040209)
98 (13 29 97 6168709 599786069)
99 (2 17 89 197 19801 18546805133)
100 (3 5 5 11 41 101 151 401 3001 570601)
cpu time: 180 real time: 1035 gc time: 70


Create a new paste based on this one


Comments: