add lenstra eliptic curve factoring minor tidy

Author: oscardssmith

src/Primes.jl

@@ -387,8 +387,20 @@ function iterate(f::FactorIterator{T}, state=(f.n, T(3))) where T
     if n <= 2^32 || isprime(n)
         return (n, 1), (T(1), n)
     end
+    # check for n=root^r
+    r = cld(ndigits(n, base=2), ndigits(Primes.PRIMES[end], base=2))
+    root = inthroot(n, r)
+    while r >= 2
+        if root^r == n
+            isprime(root) && return (root, r), (n, root+2)
+            
+        end
+        r -= 1
+        root = inthroot(n, r)
+    end
+
     should_widen = T <: BigInt || widemul(n - 1, n - 1) ≤ typemax(n)
-    p = should_widen ? pollardfactor(n) : pollardfactor(widen(n))
+    p = should_widen ? lenstrafactor(n) : lenstrafactor(widen(n))
     num_p = 0
     while true
         q, r = divrem(n, p)
@@ -520,55 +532,65 @@ julia> radical(2*2*3)
 """
 radical(n) = n==1 ? one(n) : prod(p for (p, num_p) in eachfactor(n))
 
-function pollardfactor(n::T) where {T<:Integer}
-    while true
-        c::T = rand(1:(n - 1))
-        G::T = 1
-        r::T = 1
-        y::T = rand(0:(n - 1))
-        m::T = 100
-        ys::T = 0
-        q::T = 1
-        x::T = 0
-        while c == n - 2
-            c = rand(1:(n - 1))
-        end
-        while G == 1
-            x = y
-            for i in 1:r
-                y = y^2 % n
-                y = (y + c) % n
-            end
-            k::T = 0
-            G = 1
-            while k < r && G == 1
-                ys = y
-                for i in 1:(m > (r - k) ? (r - k) : m)
-                    y = y^2 % n
-                    y = (y + c) % n
-                    q = (q * (x > y ? x - y : y - x)) % n
-                end
-                G = gcd(q, n)
-                k += m
+function lenstrafactor(n::T) where{T<:Integer}
+    # 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
+                return isprime(res) ? res : lenstrafactor(res)
             end
-            r *= 2
-        end
-        G == n && (G = 1)
-        while G == 1
-            ys = ys^2 % n
-            ys = (ys + c) % n
-            G = gcd(x > ys ? x - ys : ys - x, n)
         end
-        if G != n
-            G2 = div(n,G)
-            f = min(G, G2)
-            _gcd = gcd(G, G2)
-            if _gcd != 1
-                f = _gcd
+    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, y = T(0), T(1)
+    for B in small_primes
+        t = B^trunc(Int64, log(B, plimit))
+        xn, yn = T(0), T(0)
+        sx, sy = x, y
+
+        first = true
+        while t > 0
+            if isodd(t)
+                if first
+                    xn, yn = sx, sy
+                    first = false
+                else
+                    g, u, _ = gcdx(sx-xn, N)
+                    g > 1 && (g == N ? break : return g)
+                    u += N
+                    L  = (u*(sy-yn)) % N
+                    xs = (L*L - xn - sx) % N
+
+                    yn = (L*(xn - xs) - yn) % N
+                    xn = xs
+                end
             end
-            return isprime(f) ? f : pollardfactor(f)
+            g, u, _ = gcdx(2*sy, N)
+            g > 1 && (g == N ? break : return g)
+            u += N
+
+            L  = (u*(3*sx*sx + a)) % N
+            x2 = (L*L - 2*sx) % N
+
+            sy = (L*(sx - x2) - sy) % N
+            sx = x2
+
+            sy == 0 && return T(1)
+            t >>= 1
         end
+        x, y = xn, yn
     end
+    return T(1)
 end
 
 """