From 63fe1a0c62b7ee84f17d1815903041cedda9e89c Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Sun, 19 Dec 2021 03:01:20 -0500
Subject: [PATCH] faster and more modular

---
 src/Primes.jl | 43 ++++++++++++++++++++++++++-----------------
 1 file changed, 26 insertions(+), 17 deletions(-)

diff --git a/src/Primes.jl b/src/Primes.jl
index c5bdd9c..7c61930 100644
--- a/src/Primes.jl
+++ b/src/Primes.jl
@@ -270,7 +270,7 @@ function factor!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
         end
     end
     isprime(n) && (h[n]=1; return h)
-    T <: BigInt || widemul(n - 1, n - 1) ≤ typemax(n) ? lenstrafactors!(n, h) : lenstrafactors!(widen(n), h)
+    lenstrafactors!(widen(n), h)
 end
 
 
@@ -394,36 +394,46 @@ function inthroot(n::Integer, r::Int)
     return ans
 end
 
-function lenstrafactors!(n::T, h::AbstractDict{K,Int}, plimit=10000) where{T<:Integer,K<:Integer}
-    isprime(n) && (h[n] = get(h, n, 0)+1; return h)
+function primepowerfactor!(n::T, h::AbstractDict{K,Int}) where{T<:Integer,K<:Integer}
     r = ceil(Int, inv(log(n, Primes.PRIMES[end])))+1
     root = inthroot(n, r)
     while r >= 2
-        if root^r == n
-            isprime(root) ? h[root] = get(h, root, 0) + 1 : lenstrafactors!(root, h)
-            return lenstrafactors!(div(n, root), h)
+        if root^r == n && isprime(root)
+            h[root] = get(h, root, 0) + r
+            return true
         end
         r -= 1
         root = inthroot(n, r)
     end
-    small_primes = primes(plimit)
-    for a in zero(T):(n>>1)
-        for p in (-1, 1)
-            res = lenstra_get_factor(n, p*a, small_primes, plimit)
+    return false
+end
+
+function lenstrafactors!(n::T, h::AbstractDict{K,Int}) where{T<:Integer,K<:Integer}
+    isprime(n) && (h[n] = get(h, n, 0)+1; return h)
+    primepowerfactor!(n::T, h) && return h
+    # bounds and runs per bound taken from
+    # https://www.rieselprime.de/ziki/Elliptic_curve_method
+    B1s = Int[2e3, 11e3, 5e4, 25e4, 1e6, 3e6, 11e6,
+                43e6, 11e7, 26e7, 85e7, 29e8, 76e8, 25e9]
+    runs = (74, 221, 453, 984, 2541, 4949, 8266, 20158,
+              47173, 77666, 42057, 69471, 102212, 188056, 265557)
+    for (B1, run) in zip(B1s, runs)
+        small_primes = primes(B1)
+        for a in -run:run
+            res = lenstra_get_factor(n, a, small_primes, B1)
             if res != 1
                 isprime(res) ? h[res] = get(h, res, 0) + 1 : lenstrafactors!(res, h)
-                res2 = div(n,res)
-                isprime(res2) ? h[res2] = get(h, res2, 0) + 1 : lenstrafactors!(res2, h)
-                return h
+                n = div(n,res)
+                primepowerfactor!(n::T, h) && return h
+                isprime(n) && (h[n] = get(h, n, 0) + 1; return h)
             end
         end
     end
+    throw(ArgumentError("This number is too big to be factored with this algorith effectively"))
 end
 
 function lenstra_get_factor(N::T, a, small_primes, plimit) where T <: Integer
-    x = T(0)
-    y = T(1)
-
+    x, y = T(0), T(1)
     for B in small_primes
         t = B^trunc(Int64, log(B, plimit))
         xn, yn = T(0), T(0)
@@ -463,7 +473,6 @@ function lenstra_get_factor(N::T, a, small_primes, plimit) where T <: Integer
     end
     return T(1)
 end
-
 """
     ismersenneprime(M::Integer; [check::Bool = true]) -> Bool