codepad

Language:

; factoring with bicycle chains

(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 (square x) (* x x))

(define (every? xs)
  (if (null? xs) #t
    (if (not (car xs)) #f
      (every? (cdr xs)))))

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

(define (cycle xs)
  (set-cdr! (last-pair xs) xs) xs)

(define (chain n len)
  (let ((a (+ (isqrt n) 1))
        (qs (make-vector len #f)) ; quadratic residue?
        (cs (make-vector len #f))) ; chain link has peg?
    (do ((i 0 (+ i 1))) ((= i len))
      (vector-set! qs (modulo (square i) len) #t))
    (do ((i 0 (+ i 1)))
        ((= i len) (cycle (vector->list cs)))
      (vector-set! cs i
        (vector-ref qs (modulo (- (square (+ a i)) n) len))))))

(define (bicycle n chains)
  (let loop ((a (+ (isqrt n) 1))
             (cs (map (lambda (len) (chain n len)) chains)))
    (if (not (every? (map car cs))) (loop (+ a 1) (map cdr cs))
      (let* ((b2 (- (square a) n)) (b (isqrt b2)))
        (if (= (square b) b2) (list (- a b) (+ a b))
          (loop (+ a 1) (map cdr cs)))))))

(define chains ; lehmer's 19-chain machine
  '(64 27 25 49 22 26 17 19 23 29 31 37 41 43 47 53 59 61 67))

(display (bicycle 5237178673 chains))

Private [?]

Run code