Skip to content

Instantly share code, notes, and snippets.

@tabidots
Last active April 25, 2019 11:30
Show Gist options
  • Save tabidots/52e3ee425ca5930894f1efa20fb451f7 to your computer and use it in GitHub Desktop.
Save tabidots/52e3ee425ca5930894f1efa20fb451f7 to your computer and use it in GitHub Desktop.
[Draft of Pollard-Brent in Clojure] Attempt at writing the Pollard-Brent factorization algorithm in a functional style. Not performant though
;; 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)))))))))
# 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
(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))))))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment