add lenstra eliptic curve factoring

src/Primes.jl

@@ -249,11 +249,10 @@ Base.isempty(f::FactorIterator) = f.n == 1
 # Uses a variety of algorithms depending on the size of n to find a factor.
 #     https://en.algorithmica.org/hpc/algorithms/factorization
 # Cache of small factors for small numbers (n < 2^16)
-# Trial division of small (< 2^16) precomputed primes
-# Pollard rho's algorithm with Richard P. Brent optimizations
+# Trial division of small (< 2^16) precomputed primes and
+# Lenstra elliptic curve algorithm
 #     https://en.wikipedia.org/wiki/Trial_division
-#     https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
-#     http://maths-people.anu.edu.au/~brent/pub/pub051.html
+#     https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
 #
 
 """
@@ -316,8 +315,25 @@ 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
-    should_widen = !Base.hastypemax(T) || widemul(n - 1, n - 1) ≤ typemax(n)
-    p = should_widen ? pollardfactor(n) : pollardfactor(widen(n))
+    # check for n=root^r
+    r = cld(ndigits(n, base=2), ndigits(N_SMALL_FACTORS, base=2))
+    root = iroot(n, r)
+    while r >= 2
+        if root^r == n
+            if isprime(root)
+                return (root, r), (1, root+2)
+            else
+                f = first(eachfactor(root))
+                return f, (n÷f, root+2)
+            end
+        end
+        r -= 1
+        root = iroot(n, r)
+    end
+
+    # lenstrafactor's modular arithmetic needs intermediates up to 2n^2
+    fits = !Base.hastypemax(T) || widemul(n, n) ≤ typemax(n) ÷ 2
+    p = fits ? lenstrafactor(n) : lenstrafactor(widen(n))
     num_p = 0
     while true
         q, r = divrem(n, p)
@@ -327,14 +343,157 @@ function iterate(f::FactorIterator{T}, state=(f.n, T(3))) where T
     end
 end
 
+# the curve arithmetic uses plain +/- with negative intermediates, so run
+# unsigned inputs in the equally-sized signed domain (the widening check
+# already caps n^2 ≤ typemax(n)÷2, so n fits the signed type)
+lenstrafactor(n::Unsigned) = lenstrafactor(signed(typeof(n))(n))
+
+function lenstrafactor(n::T; use_stage2::Bool=true) where{T<:Integer}
+    # bounds and runs per bound taken from
+    # https://www.rieselprime.de/ziki/Elliptic_curve_method
+    # the 200 tier is a cheap pre-pass that catches most factors ≤ ~20 bits
+    B1s = Int[200, 2e3, 11e3, 5e4, 25e4, 1e6, 3e6, 11e6,
+                43e6, 11e7, 26e7, 85e7, 29e8, 76e8, 25e9]
+    # published counts assume GMP-ECM's much larger B2; scaled 1.5x for B2 = 25*B1
+    runs = Int[6, 38, 135, 450, 1050, 2700, 7650, 2700, 15900,
+              29000, 73500, 186000, 315000, 510000, 15*10^5, 15*10^6]
+    for (B1, run) in zip(B1s, runs)
+        small_primes = primes(B1)
+        for σ in 6:2*run+6
+            res = lenstra_stage_1(n, σ, small_primes, B1; use_stage2)
+            if res != 1
+                return isprime(res) ? res : lenstrafactor(res; use_stage2)
+            end
+        end
+    end
+    throw(ArgumentError("This number is too big to be factored with this algorithm effectively"))
+end
+
+# x-only arithmetic on Montgomery curves B*y^2 = x^3 + A*x^2 + x, projective X:Z.
+# mod-reduced values live in [0, N), but intermediate sums/differences may reach
+# (-N, 2N), so this arithmetic needs signed T and |products| ≤ 2N^2 to fit in T —
+# both guaranteed by the callers (widening check + unsigned-to-signed conversion).
+# _addmod keeps sums that get squared in [0, N); squaring an unreduced sum would
+# need 4N^2 of headroom, i.e. cost half a bit of the widening threshold.
+@inline _addmod(a, b, N) = (c = a + b; c ≥ N ? c - N : c)
+
+# doubling, with a24 = (A+2)/4
+@inline function _xdbl(X, Z, a24, N)
+    s = _addmod(X, Z, N)
+    d = X - Z
+    s2 = mod(s*s, N)
+    d2 = mod(d*d, N)
+    t = s2 - d2     # 4XZ
+    return mod(s2*d2, N), mod(t*(d2 + mod(a24*t, N)), N)
+end
+
+# differential addition: P+Q given their difference D = P-Q
+@inline function _xadd(XP, ZP, XQ, ZQ, XD, ZD, N)
+    u = mod((XP - ZP)*(XQ + ZQ), N)
+    v = mod((XP + ZP)*(XQ - ZQ), N)
+    s = _addmod(u, v, N)
+    d = u - v
+    return mod(ZD*mod(s*s, N), N), mod(XD*mod(d*d, N), N)
+end
+
+# Montgomery ladder: k*(X:Z) for k ≥ 1
+@inline function _xladder(k, X, Z, a24, N)
+    k == 1 && return X, Z
+    X0, Z0 = X, Z
+    X1, Z1 = _xdbl(X, Z, a24, N)
+    for i in 8*sizeof(k)-leading_zeros(k)-2:-1:0
+        if (k >> i) & 1 == 0
+            X1, Z1 = _xadd(X0, Z0, X1, Z1, X, Z, N)
+            X0, Z0 = _xdbl(X0, Z0, a24, N)
+        else
+            X0, Z0 = _xadd(X1, Z1, X0, Z0, X, Z, N)
+            X1, Z1 = _xdbl(X1, Z1, a24, N)
+        end
+    end
+    return X0, Z0
+end
+
+# One ECM stage-1 attempt on the Suyama curve with parameter σ ≥ 6.
+# Returns a nontrivial factor of N, or 1 if this curve fails.
+function lenstra_stage_1(N::T, σ, small_primes, plimit; use_stage2::Bool=true) where T <: Integer
+    # Suyama parametrization: guarantees torsion 12 and a known starting point
+    u = mod(T(σ)*T(σ) - 5, N)
+    v = mod(4*T(σ), N)
+    X = mod(mod(u*u, N)*u, N)   # u^3
+    Z = mod(mod(v*v, N)*v, N)   # v^3
+    # a24 = (v-u)^3*(3u+v) / (16*u^3*v); a non-invertible denominator is itself a factor
+    den = mod(mod(T(16)*X, N)*v, N)
+    g, di, _ = gcdx(den, N)
+    g > 1 && return g == N ? T(1) : g
+    w = v - u
+    num = mod(mod(mod(w*w, N)*w, N)*mod(3u + v, N), N)
+    a24 = mod(num*mod(di, N), N)
+    for B in small_primes
+        # multiply in the largest power of B that stays ≤ plimit
+        t = B*B > plimit ? B : B^trunc(Int64, log(B, plimit))
+        X, Z = _xladder(t, X, Z, a24, N)
+        Z == 0 && return T(1)  # hit the identity mod N: curve failed
+    end
+    g = gcd(Z, N)
+    1 < g < N && return g
+    g == N && return T(1)
+    return use_stage2 ? _lenstra_stage2(N, X, Z, a24, plimit) : T(1)
+end
+
+# Stage 2 (standard continuation): look for a single prime q ∈ (B1, B2] with
+# q*Q ≡ identity mod p. Writing q = m*D ± j, that means x(mD*Q) = x(j*Q), i.e.
+# p | X_m*Z_j - X_j*Z_m; accumulate those cross products and gcd once at the end.
+# Tests every mD ± j with gcd(j, D) = 1 rather than sieving primes out of (B1, B2] —
+# a ~3x overhead that avoids enumerating primes up to B2.
+# B2 = 25*B1 puts stage-2 cost at ~half of stage 1, which is throughput-optimal for
+# this linear-cost continuation: marginal returns decay like 1/(q*ln q) while the
+# cost per unit B2 is flat, so the optimum sits below cost balance.
+function _lenstra_stage2(N::T, X, Z, a24, B1) where T <: Integer
+    B2 = 25 * B1
+    D = 210
+    # baby steps: j*Q for odd j ≤ D/2, chained by Q_{j+2} = Q_j + 2Q (diff Q_{j-2})
+    nb = (D ÷ 2 + 1) ÷ 2
+    Xb = Vector{T}(undef, nb)
+    Zb = Vector{T}(undef, nb)
+    Xb[1], Zb[1] = X, Z
+    X2, Z2 = _xdbl(X, Z, a24, N)
+    Xp, Zp = X, Z
+    Xj, Zj = _xadd(X2, Z2, X, Z, X, Z, N)  # 3Q
+    Xb[2], Zb[2] = Xj, Zj
+    for idx in 3:nb
+        Xn, Zn = _xadd(Xj, Zj, X2, Z2, Xp, Zp, N)
+        Xp, Zp = Xj, Zj
+        Xj, Zj = Xn, Zn
+        Xb[idx], Zb[idx] = Xj, Zj
+    end
+    # giant steps: R = m*D*Q, stepping by D*Q
+    XD, ZD = _xladder(D, X, Z, a24, N)
+    m0 = fld(B1, D)
+    XRp, ZRp = _xladder(m0 - 1, XD, ZD, a24, N)
+    XR, ZR = _xladder(m0, XD, ZD, a24, N)
+    acc = T(1)
+    for m in m0:cld(B2, D)
+        for idx in 1:nb
+            j = 2idx - 1
+            gcd(j, D) == 1 || continue
+            acc = mod(acc*(mod(XR*Zb[idx], N) - mod(Xb[idx]*ZR, N)), N)
+        end
+        acc == 0 && return T(1)  # N divides the product: factor unrecoverable, next curve
+        XRn, ZRn = _xadd(XR, ZR, XD, ZD, XRp, ZRp, N)
+        XRp, ZRp = XR, ZR
+        XR, ZR = XRn, ZRn
+    end
+    g = gcd(acc, N)
+    return 1 < g < N ? g : T(1)
+end
+
 function factor!(n::T, h::AbstractDict{K,Int}) where {T<:Integer,K<:Integer}
     for (p, num_p) in eachfactor(n)
         increment!(h, num_p, p)
     end
     return h
 end
 
-
 """
     factor(n::Integer) -> Primes.Factorization
 
@@ -449,58 +608,6 @@ 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
-            end
-            yield()
-            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
-            return isprime(f) ? f : pollardfactor(f)
-        end
-    end
-end
-
 """
     ismersenneprime(M::Integer; [check::Bool = true]) -> Bool
 
@@ -894,19 +1001,13 @@ prevprimes(start::T, n::Integer) where {T<:Integer} =
 
 """
     divisors(n::Integer) -> Vector
-
 Return a vector with the positive divisors of `n`.
-
 For a nonzero integer `n` with prime factorization `n = p₁^k₁ ⋯ pₘ^kₘ`, `divisors(n)`
 returns a vector of length (k₁ + 1)⋯(kₘ + 1) containing the divisors of `n` in
 lexicographic (rather than numerical) order.
-
 `divisors(-n)` is equivalent to `divisors(n)`.
-
 For convenience, `divisors(0)` returns `[]`.
-
 # Example
-
 ```jldoctest; filter = r"(\\s+#.*)?"
 julia> divisors(60)
 12-element Vector{Int64}:
@@ -922,14 +1023,12 @@ julia> divisors(60)
  15         # 5 * 3
  30         # 5 * 3 * 2
  60         # 5 * 3 * 2 * 2
-
 julia> divisors(-10)
 4-element Vector{Int64}:
   1
   2
   5
  10
-
 julia> divisors(0)
 Int64[]
 ```
@@ -947,7 +1046,6 @@ end
 
 """
     divisors(f::Factorization) -> Vector
-
 Return a vector with the positive divisors of the number whose factorization is `f`.
 Divisors are sorted lexicographically, rather than numerically.
 """

test/runtests.jl

@@ -186,6 +186,12 @@ end
 # test second branch, after all small primes in list have been searched
 @test factor(10009 * Int128(1000000000000037)) == Dict(10009 => 1, 1000000000000037 => 1)
 
+# 62-bit semiprimes exercise lenstrafactor's widened (Int128) modular arithmetic
+@test factor(2147483659 * 2147484679) == Dict(2147483659 => 1, 2147484679 => 1)
+@test factor(2147539211 * 3221125483) == Dict(2147539211 => 1, 3221125483 => 1)
+# unsigned inputs route through lenstrafactor's signed-domain conversion
+@test factor(UInt64(2147483659) * UInt64(2147484679)) == Dict(2147483659 => 1, 2147484679 => 1)
+
 for n = 1:100
     m = 1
     for (p,k) in factor(n)