Skip to content

Instantly share code, notes, and snippets.

@tabidots
Last active April 2, 2019 14:46
Show Gist options
  • Save tabidots/aaef61a9e3418b74d18e5ef88604e997 to your computer and use it in GitHub Desktop.
Save tabidots/aaef61a9e3418b74d18e5ef88604e997 to your computer and use it in GitHub Desktop.
[perfect-powers] Testing for perfect powers #number-theory
(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)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment