Skip to content

Instantly share code, notes, and snippets.

@tabidots
Created April 15, 2019 04:33
Show Gist options
  • Save tabidots/9bd38be86cf66ef7ee74927b431152d5 to your computer and use it in GitHub Desktop.
Save tabidots/9bd38be86cf66ef7ee74927b431152d5 to your computer and use it in GitHub Desktop.
[sieve-factorization] Attempt to do prime factorization using JVM array of smallest prime factors
(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)))))))
(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))))))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment