```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 ``` ```; dixon's factorization algorithm (define (cons* first . rest) (let loop ((curr first) (rest rest)) (if (null? rest) curr (cons curr (loop (car rest) (cdr rest)))))) (define (isqrt n) (if (not (and (positive? n) (integer? n))) (error 'isqrt "must be positive integer") (let loop ((x n)) (let ((y (quotient (+ x (quotient n x)) 2))) (if (< y x) (loop y) x))))) (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 (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 (make-factor-base n k bound lim) (let loop ((i (- lim 1)) (ps (cdr (primes bound))) (fb '(2))) (cond ((or (zero? i) (null? ps)) (reverse fb)) ((not (negative? (jacobi (* k n) (car ps)))) (loop (- i 1) (cdr ps) (cons (car ps) fb))) (else (loop i (cdr ps) fb))))) (define (smooth n q fb) (let loop ((q q) (fb fb) (fs '()) (i 0)) (cond ((= q 1) (reverse fs)) ((null? fb) #f) ((zero? (modulo q (car fb))) (if (and (pair? fs) (= (car fb) (car fs))) (loop (/ q (car fb)) fb (cdr fs) (+ i 1)) (loop (/ q (car fb)) fb (cons (car fb) fs) (+ i 1)))) (else (loop q (cdr fb) fs (+ i 1)))))) ;(define (residue n k fb lim) ; (let* ((kn (* k n)) (g (isqrt kn)) (a-3 0) (a-2 1) ; (q-2 kn) (r-2 g) (g+p-1 g) (q-1 1) ; (t-1 (quotient g+p-1 q-1)) (r-1 (- g+p-1 (* q-1 t-1))) ; (q0 (+ q-2 (* t-1 (- r-1 r-2))))) ; (let loop ((a-2 1) (a-1 g) (q-1 1) (r-1 0) (g+p (+ g g)) (q q0) (i 1) (qs '()) (lim lim)) ; ; (for-each display `(,i " " ,g+p " " ,q " " ,a-1 #\newline)) ; (if (or (= q 1) (zero? lim)) (reverse qs) ; (let* ((t (quotient g+p q)) ; (r (- g+p (* q t))) ; (a (modulo (+ (* t a-1) a-2) n)) ; (q+1 (+ q-1 (* t (- r r-1)))) ; (fs (smooth n q fb))) ; (cond ((null? fs) ; (let ((d (gcd (- a-1 (isqrt q)) n))) ; (if (< 1 d n) d ; (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1))))) ; (fs (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) (cons (cons* i q a-1 fs) qs) (- lim 1))) ; (else (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1))))))))) (define (residue n k fb lim) (let loop ((r (isqrt n)) (qs '())) (if (zero? lim) (reverse qs) (let* ((r2n (modulo (* r r) n)) (fs (smooth n r2n fb))) (cond ((null? fs) (let ((d (gcd (- r (isqrt r2n)) n))) (if (< 1 d n) d (loop (+ r 1) qs)))) (fs (loop (+ r 1) (cons (cons* 0 r2n r fs) qs))) (else (loop (+ r 1) qs))))))) ; (fibonacci 73) => 9375829 * 86020717 (define fb (make-factor-base 806515533049393 1 60000 270)) (define residues (residue 806515533049393 1 fb 2000000)) (define (make-row rs fb len) (let ((vec (make-vector (+ len 1) 0))) (vector-set! vec 0 (if (odd? (car rs)) 1 0)) (let loop ((i 1) (rs (cdddr rs)) (fb fb)) ; (for-each display `(,i " " ,rs " " ,fb " " #\newline)) (cond ((or (null? rs) (null? fb)) vec) ((= (car fb) (car rs)) (vector-set! vec i 1) (loop (+ i 1) (cdr rs) (cdr fb))) (else (loop (+ i 1) rs (cdr fb))))))) (define (make-ev fb len residues) (let loop ((rss residues) (es '()) (as '()) (qs '())) (if (null? rss) (values (list->vector (reverse es)) (list->vector (reverse as)) (list->vector (reverse qs))) (loop (cdr rss) (cons (make-row (car rss) fb len) es) (cons (caddar rss) as) (cons (cadar rss) qs))))) (define (make-hv ev-len) (let ((hv (make-vector ev-len))) (do ((i 0 (+ i 1))) ((= i ev-len) hv) (let ((h (make-vector ev-len 0))) (vector-set! h i 1) (vector-set! hv i h))))) (define (right-most-one rv start) (let loop ((i start)) (cond ((negative? i) i) ((= (vector-ref rv i) 1) i) (else (loop (- i 1)))))) (define (vector-add v1 v2) (let* ((len (vector-length v1)) (v3 (make-vector len))) (do ((i 0 (+ i 1))) ((= i len) v3) (vector-set! v3 i (if (= (vector-ref v1 i) (vector-ref v2 i)) 0 1))))) (define (vector-sum vec) (do ((i 0 (+ i 1)) (sum 0 (+ sum (vector-ref vec i)))) ((= i (vector-length vec)) sum))) (define (a-minus-q n vec ev-len av qv) (let loop ((i 0) (a-prod 1) (q-prod 1)) (cond ((= i ev-len) (modulo (- a-prod (isqrt q-prod)) n)) ((= (vector-ref vec i) 1) (loop (+ i 1) (modulo (* a-prod (vector-ref av i)) n) (* q-prod (vector-ref qv i)))) (else (loop (+ i 1) a-prod q-prod))))) (define (answer n fb residues) (let-values (((ev av qv) (make-ev fb (length fb) residues))) (let* ((fb-len (length fb)) (ev-len (vector-length ev)) (hv (make-hv ev-len)) (ptr (make-vector ev-len))) (do ((i 0 (+ i 1))) ((= i ev-len)) (vector-set! ptr i (right-most-one (vector-ref ev i) fb-len))) ; reduce ev and hv (do ((j fb-len (- j 1))) ((negative? j)) (let loop ((i 0)) (when (< i ev-len) (if (= (vector-ref ptr i) j) (do ((m (+ i 1) (+ m 1))) ((= m ev-len)) (when (= (vector-ref ptr m) j) (vector-set! ev m (vector-add (vector-ref ev i) (vector-ref ev m))) (vector-set! hv m (vector-add (vector-ref hv i) (vector-ref hv m))) (vector-set! ptr m (right-most-one (vector-ref ev m) j)))) (loop (+ i 1)))))) ; examine s-sets (let loop ((i 0)) (cond ((= i ev-len) #f) ((zero? (vector-sum (vector-ref ev i))) (let ((d (gcd (a-minus-q n (vector-ref hv i) ev-len av qv) n))) (if (< 1 d n) d (loop (+ i 1))))) (else (loop (+ i 1)))))))) (answer 806515533049393 fb residues) ```