[Draft of Pollard-Brent in Clojure] Attempt at writing the Pollard-Brent factorization algorithm in a functional style. Not performant though

brent.clj

;; Hard-coded values for c, m, x here to debug the code and compare performance vs. Python implementation
;; Performance is VERY poor (~4 min). This is the fifth rewrite, and the most functional 
;; (imperative attempts did not fare much better, and were incredibly ugly).
;; However, with the same seed values, this does perform the algorithm correctly 
;; (that is, identically to the fast Python implementation below), finding a factor at some point in the loop
;; between r = 1048576 and r = 2097152.

;; With reference to Brent's paper "An Improved Monte Carlo Factorization Algorithm
;; https://maths-people.anu.edu.au/~brent/pd/rpb051i.pdf

(require '[clojure.core.reducers :as r])
(require '[clojure.math.numeric-tower :refer [abs])
         
;; For the other functions, see helpers.clj below

(let [n 5687450887724258126853054080419196271544306051]
  (if (even? n) 2
    (let [c  1384618644036103387726606019273596888099126033 ;(rand-num 1 (dec n))
          f  (fn [x] (-> (+ c (mod-pow x 2 n)) (mod n)))
          m  163544 ;(rand-num 1 (dec n))
          x  117210326061389040030818806148762920404935678]  ;(rand-num 1 (dec n))
      (loop [r 1 q 1 x x]                               ;; x  = Tortoise value
        (let [y   (f x)                                 ;; y  = Backup value if algo fails
              ys  (->> (iterate f y)                    ;; ys = Sequence of hare values
                       (r/take r)
                       (r/foldcat))
              zs  (r/foldcat (r/map #(abs (- x %)) ys)) ;; zs = Differences between tortoise & hare values
              qs  (reductions (fn [res z]               ;; qs = Cumulative product of differences
                                (mod-mul res z n))      ;;      all the way from the beginning
                              q zs)
              g   (or (->> (take-nth m qs)              ;; m = Step size. We are only interested in
                           (r/map #(gcd % n))           ;;     computing gcd(q, n) at the steps that are
                           (r/filter #(> % 1))          ;;     multiples of the step size.
                           (r/foldcat)                  ;;
                           first)                       ;;     If none of them are 1, leave g = 1 so the
                      1)]                               ;;     loop continues.
          (cond
            (< 1 g n) g                                 ;; Success
            (= g n)   (let [second-try (->> (iterate f y)
                                            (r/map #(gcd (- x %) n))
                                            (r/filter #(> % 1))
                                            (r/foldcat)
                                            first)]
                        (when (< 1 second-try n) second-try))
            :else     (do (println r)
                        (recur (*' 2 r)
                               (last qs)
                               (last ys)))))))))

brent.py

# Source: https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
# This implementation factors the test input in ~3s!

from math import gcd
import random

def brent(N): # Test with N = 5687450887724258126853054080419196271544306051
        if N%2==0:
                return 2
        # y,c,m = random.randint(1, N-1),random.randint(1, N-1),random.randint(1, N-1)
        y,c,m = 117210326061389040030818806148762920404935678, 1384618644036103387726606019273596888099126033, 163544
        g,r,q = 1,1,1
        while g==1:             
                x = y
                for i in range(r):
                        y = ((y*y)%N+c)%N
                k = 0
                while (k<r and g==1):
                        ys = y
                        for i in range(min(m,r-k)):
                                y = ((y*y)%N+c)%N
                                q = q*(abs(x-y))%N
                        g = gcd(q,N)
                        k = k + m
                r = r*2
        if g==N:
                while True:
                        ys = ((ys*ys)%N+c)%N
                        g = gcd(abs(x-ys),N)
                        if g>1:
                                break
         
        return g   

helpers.clj

(import 'java.util.concurrent.ThreadLocalRandom)

(defn rand-num
  "Random number generator of arbitrary-precision integers."
  ([max]
   (rand-num 0 max))
  ([min max]
   (let [src (ThreadLocalRandom/current)]
     (cond
       ;; Clojure translation of https://stackoverflow.com/a/2290089
       (> max Long/MAX_VALUE)    (->> #(BigInteger. (.bitLength (bigint max)) src)
                                      repeatedly
                                      (filter #(<= min % max))
                                      first)
       ;; For longs instead of ints, use ThreadLocalRandom.current().nextLong(start, end)
       ;; Clojure High Performance Programming, p. 75
       (> max Integer/MAX_VALUE) (.nextLong src min max)
       :else                     (.nextInt src min max)))))

(defn mod-mul
  ;; Translated from https://www.geeksforgeeks.org/multiply-large-integers-under-large-modulo/
  "Quickly calculates a * b % m. Useful when a and b are very large integers."
  [a b m]
  (if (or (= m 1) (= a m)) 0
    (loop [a (mod a m) b b res 0]
      (if (zero? b) res
        (recur (-> (+' a a) (mod m))
               (.shiftRight (biginteger b) 1)
               (if (odd? b)
                 (-> (+' res a) (mod m))
                 res))))))

(defn mod-pow
  ;; Adapted from https://en.wikipedia.org/wiki/Modular_exponentiation
  "Quickly calculates a ^ b % m. Useful when a and b are very large integers."
  [base exp m]
  (if (or (= m 1) (= base m)) 0
    (loop [base (mod base m) e exp res 1]
      (if (zero? e) res
        (recur (-> (*' base base) (mod m))
               (.shiftRight (biginteger e) 1)
               (if (odd? e)
                 (-> (*' res base) (mod m))
                 res))))))