[perfect-powers] Testing for perfect powers #number-theory

perfect_powers.clj

(defn babylonian-root
  "High-precision BigDecimal nth-root using the Babylonian algorithm,
  with a close initial approximation for ridiculously fast convergence."
  [n root]
  (let [eps 0.000000000000000000000000000000000000000001M]
    (loop [t (bigdec (tower/expt n (/ root)))] ;; rough initial approx
      (let [ts  (repeat (- root 1) t)
            nxt (with-precision 100 (-> (/ n (reduce *' ts))
                                        (+ (reduce +' ts))
                                        (/ root)))]
        (if (< (.abs (- nxt t)) eps) t
          (recur nxt))))))

(defn nth-root-is-integer? [x n]
  (let [nth-root (babylonian-root x n)
        floor    (.setScale nth-root 0 BigDecimal/ROUND_HALF_EVEN)
        exp      (bigint (reduce *' (repeat n floor)))]
    (= x exp)))

;; Some test values
;; 6221821273427820544 = 484^7
;; 92709463147897837085761925410587 = 3^67

(defn perfect-power? [n]
  (let [max-power (Math/log n) ;  https://cr.yp.to/papers/powers-ams.pdf
        powers    (cons 2 (filter odd? (range 3 (inc max-power))))]
        ;; Strictly speaking, this should be prime powers only, but
        ;; that would require additional overhead by filtering primes
        ;; which seems circuitous since the purpose of this function is already
        ;; the first step in an algorithm to test the primality of a gigantic
        ;; integer. Plus, removing the prime? filter allows this function to
        ;; work without specifying a custom prime? function.
        ;;
        ;; Otherwise, use (filter prime? (range 2 (inc max-power)))
    (some (partial nth-root-is-integer? n) powers)))