[sieve-factorization] Attempt to do prime factorization using JVM array of smallest prime factors

factorization.clj

(ns numthy.factorization
  (:require [clojure.math.numeric-tower :as tower]
            [numthy.helpers :as h]))

(defn spf-sieve
  [n]
  ;; Adapted from https://stackoverflow.com/a/48137788/4210855
  (let [^ints sieve (int-array n (range n))
        upper-bound (h/isqrt n)]
    (loop [p 2] ;; p's are the bases
      (if (> p upper-bound) sieve
        (do
          (when (= p (aget sieve p))
            (loop [i (* 2 p)] ;; i's are the multiples of p
              (when (< i n)
                (aset sieve i (min (aget sieve i) p))
                (recur (+ i p)))))
          (recur (inc p)))))))

(defn primes-to
  [n]
  (->> (spf-sieve n)
       (keep-indexed (fn [i x] (when (= i x) i)))
       (drop 2)
       (into [])))

(time
  (let [n     659119 ;5394992290384 exceeds Java int limit
        sieve (spf-sieve (inc n))]
    (loop [n n d (aget sieve n) fs (sorted-set)]
      (if (= 1 d) fs
        (let [f (/ n d)]
          (recur f (aget sieve f) (conj fs d)))))))

helpers.clj

(ns numthy.helpers)

(defn isqrt
  "floor(√n). When incremented, provides an upper bound for factorization."
  ;; Java interop is super fast but not accurate for n > 1E24 (approx) due to
  ;; floating-point rounding. Uses a slightly slower but pinpoint-precise method for n > 1E24.
  [n]
  (if (< n 1E24)
    (-> (Math/sqrt n) bigint)
    ;; https://cs.stackexchange.com/a/30383
    (let [half-bit-length (quot (.bitLength (bigint n)) 2)]
      (loop [a (tower/expt 2 half-bit-length)
             b a
             c (*' a a)]
        (cond
          (zero? b) a
          (> c n)   (recur (-' a b) (quot b 2) (+ c (*' -2 a b) (*' b b)))
          :else     (recur (+' a b) (quot b 2) (+ c (*' 2 a b) (*' b b))))))))