; shanks's squfof
(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?
(let ((q11 (make-vector 11 #f))
(q63 (make-vector 63 #f))
(q64 (make-vector 64 #f))
(q65 (make-vector 65 #f)))
(do ((k 0 (+ k 1))) ((< 5 k))
(vector-set! q11 (modulo (* k k) 11) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q63 (modulo (* k k) 63) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q64 (modulo (* k k) 64) #t))
(do ((k 0 (+ k 1))) ((< 32 k))
(vector-set! q65 (modulo (* k k) 65) #t))
(lambda (n)
(if (not (integer? n)) (error 'square? "must be integer")
(if (< n 1) #f
(if (not (vector-ref q64 (modulo n 64))) #f
(let ((r (modulo n 45045)))
(if (not (vector-ref q63 (modulo r 63))) #f
(if (not (vector-ref q65 (modulo r 65))) #f
(if (not (vector-ref q11 (modulo r 11))) #f
(let ((q (isqrt n)))
(if (= (* q q) n) q #f))))))))))))
(define (prime? n)
(define (expm 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)))))
(define (digits n . args)
(let ((b (if (null? args) 10 (car args))))
(let loop ((n n) (d '()))
(if (zero? n) d
(loop (quotient n b)
(cons (modulo n b) d))))))
(define (isqrt n)
(let loop ((x n) (y (quotient (+ n 1) 2)))
(if (<= 0 (- y x) 1) x
(loop y (quotient (+ y (quotient n y)) 2)))))
(define (square? n)
(let ((n2 (isqrt n)))
(= n (* n2 n2))))
(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 legendre jacobi)
(define (inverse x n)
(let loop ((x (modulo x n)) (a 1))
(cond ((zero? x) (error 'inverse "division by zero"))
((= x 1) a)
(else (let ((q (- (quotient n x))))
(loop (+ n (* q x)) (modulo (* q a) n)))))))
(define (miller? 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)))))))))
(define (chain m f g x0 x1)
(let loop ((ms (digits m 2)) (u x0) (v x1))
(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))))))
(define (lucas? 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 (legendre 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)))))))))
(cond ((or (not (integer? n)) (< n 2))
(error 'prime? "must be integer greater than one"))
((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))
(else (and (miller? n 2) (miller? n 3) (lucas? n)))))
(define (rho f)
(let* ((a (car f)) (b (cadr f)) (c (caddr f))
(d (- (* b b) (* 4 a c)))
(s (isqrt d)) (l (abs c))
(r (if (< l (quotient s 2))
(- (* 2 l (quotient (+ b s) (* 2 l))) b)
(- (* 2 l) b))))
(list c r (quotient (- (* r r) d) (* 4 c)))))
(define (reduce f)
(let* ((a (car f)) (b (cadr f)) (c (caddr f))
(d (- (* b b) (* 4 a c))) (s (isqrt d)))
(if (< (abs (- s (* 2 (abs a)))) b s) f (reduce (rho f)))))
(define (squfof n)
(cond ((prime? n) n) ((even? n) 2) ((square? n) => identity)
(else (let* ((n4 (= (modulo n 4) 1)) (big-d (if n4 n (* 4 n)))
(m (if n4 1 2)) (little-d (isqrt big-d))
(b (if n4
(+ (* 2 (quotient (- little-d 1) 2)) 1)
(* 2 (quotient little-d 2))))
(delta (quotient (- (* b b) big-d) 4))
(f (list 1 b delta))
(l (isqrt little-d)) (bound (* 4 l))
(q (let* ((x (abs delta)) (g (quotient x (gcd x m))))
(if (<= g l) (list g) '()))))
(let loop ((i 2) (f (rho f)) (q q))
(cond ((< bound i) #f)
((square? (caddr f)) =>
(lambda (c)
(if (not (member c q))
(let loop ((g (reduce (list (- c) (cadr f)
(* -1 c (car f))))))
(let ((b-prime (cadr g)) (new-g (rho g)))
(if (= (cadr new-g) b-prime)
(let ((c (caddr g)))
(abs (if (odd? c) c (/ c 2))))
(loop new-g))))
(if (= c 1) #f
(let* ((x (abs (caddr f)))
(g (quotient x (gcd x m))))
(loop (+ i 1) (rho f)
(if (and (<= g l) (not (member g q)))
(cons g q)
q)))))))
(else (let* ((x (abs (caddr f)))
(g (quotient x (gcd x m))))
(loop (+ i 1) (rho f)
(if (and (<= g l) (not (member g q)))
(cons g q)
q))))))))))
(display (squfof 633003781))