; sieving a polynomial (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 (primes n) ; list of primes not exceeding n (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t))) (let loop ((i 0) (p 3) (ps (list 2))) (cond ((< n (* p p)) (do ((i i (+ i 1)) (p p (+ p 2)) (ps ps (if (vector-ref bits i) (cons p ps) ps))) ((= i len) (reverse ps)))) ((vector-ref bits i) (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p))) ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps))) (vector-set! bits j #f))) (else (loop (+ i 1) (+ p 2) ps)))))) (define (factors-1mod4 n) (let* ((len (quotient (- n 1) 4)) (factors (make-vector (+ len 1) (list)))) (define (sieve p pp) ; p = prime, pp = prime-power (let* ((k (if (= (modulo pp 4) 1) pp (* 3 pp))) (k (quotient (- k 1) 4))) (do ((i k (+ i pp))) ((< len i)) (vector-set! factors i (cons p (vector-ref factors i)))))) (do ((ps (cdr (primes (isqrt n))) (cdr ps))) ((null? ps)) (do ((pp (car ps) (* pp (car ps)))) ((< n pp)) (sieve (car ps) pp))) (do ((k 0 (+ k 1))) ((< len k) factors) (let* ((n (+ (* 4 k) 1)) (fs (vector-ref factors k)) (f (/ n (apply * fs)))) (vector-set! factors k (cond ((null? fs) (list n)) ((= f 1) (reverse fs)) (else (reverse (cons f fs))))))))) (define fs1mod4 (factors-1mod4 10001)) (do ((k 1 (+ k 1))) ((= (vector-length fs1mod4) k)) (display (+ (* 4 k) 1)) (display #\tab) (display (vector-ref fs1mod4 k)) (newline))