codepad

Language:

; the factorization of f7, part 1

(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) ; new version
  (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)))
          (if fs
              (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) (cons (cons* i q a-1 fs) qs) (- lim 1))
              (loop a-1 a q r (- (* 2 g) r) q+1 (+ i 1) qs (- lim 1))))))))

(for-each
  (lambda (x) (display x) (newline))
  (residue 13290059 1 '(2 5 31 41 43 53 113) 44))

Private [?]

Run code