; two powering predicates

(define (ilog b n)
  (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b))
    (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi))
      (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi))
        (if (<= (- hi lo) 1) (if (= b^hi n) hi lo)
          (let* ((mid (quotient (+ lo hi) 2))
                 (b^mid (* b^lo (expt b (- mid lo)))))
            (cond ((< n b^mid) (loop2 lo b^lo mid b^mid))
                  ((< b^mid n) (loop2 mid b^mid hi b^hi))
                  (else mid))))))))

(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 (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 (isqrt n)
  (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 (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-power? n)
  (cond ((not (integer? n)) (error 'prime-power? "must be integer"))
        ((< n 1) #f) ((= n 1) 1) ((prime? n) n)
        ((even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f))
        (else (let loop ((a 2) (n n))
                (let* ((b (expm a n n)) (p (gcd (- b a) n)))
                  (cond ((or (= p 1) (= a b)) #f)
                        ((prime? p) (if (= (expt p (ilog p n)) n) p #f))
                        (else (loop (+ a 1) p))))))))

(display (isqrt 24)) (newline)
(display (isqrt 25)) (newline)

(display (square? 24)) (newline)
(display (square? 25)) (newline)

(display (prime-power? 8)) (newline)
(display (prime-power? 45)) (newline)
(display (prime-power? 561)) (newline)
(display (prime-power? 12167)) (newline)
(display (prime-power? 8388609)) (newline)

4
5
#f
5
2
#f
#f
23
#f