[euler437]

euler437.clj

(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 problem-437
  [lim]
  (let [sieve           (spf-sieve (inc lim))
        primes          (->> sieve
                             (keep-indexed (fn [i x] (when (= i x) i)))
                             (drop 2))
        get-spf*        (fn [n] (aget sieve n))
        get-spf         (memoize get-spf*)
        prime-factors   (fn [n]
                          (loop [n n d (get-spf n) fs (sorted-set)]
                            (if (= 1 d) fs
                              (let [f (/ n d)]
                                (recur f (get-spf f) (conj fs d))))))
        primitive-root? (fn [x p]
                          (let [phi (dec p) pfs (prime-factors phi)]
                            (not-any? #(= 1 (ma/mod-pow x (/ phi %) p)) pfs)))
        has-fpr?        (fn [p]
                          (when-let [roots (qr/mod-quadratic-zeroes {2 1, 1 -1, 0 -1} p)]
                            (some #(primitive-root? % p) roots)))]
    (->> (filter has-fpr? primes)
         (cons 5)
         (reduce +'))))