Category talk:Wren-ecm
Lenstra elliptic-curve factorization
This implementation is based on the Python code here which doesn't appear to be subject to any license.
Although ECM is said to be the third fastest factorization algorithm (after MPQS and GNFS), the Wren implementation below is usually much slower in practice than a combination of Pollard's Rho and wheel based factorization which is why I've used the latter in Wren-gmp. This may be because the ECM algorithm is far more complicated and memory intensive than Pollard's Rho and the cost of interpreting the additional code may outweigh any speed advantage it would otherwise have.
Pollard's Rho may nevertheless take a long time to factor certain numbers (generally those whose factors are all large) and this is where ECM comes in useful. However, the implementation uses a prime sieve and you may experience memory problems if you try to factorize numbers greater than 60 digits long after all small factors have been removed.
Where it is clear that simpler methods will find factor(s) as quick or quicker than ECM these are used instead.
The ECM class contains its own prime sieve and binary search methods to avoid importing these from other modules.
Finally, it should be noted that the algorithm depends on picking an elliptic curve at random and is therefore non-deterministic. This may result in execution time varying over quite a large range when factorizing a particular number. It also may result in the factors not always being found in the same order.
Source code
/* Module "ecm.wren" */
import "./gmp" for Mpz
import "random" for Random
var MaxCurves = 10000
var MaxRnd = 1 << 31
var MaxB1 = 43 * 1e7
var MaxB2 = 2 * 1e10
var Rand = Random.new()
class Ecm {
// Private method to compute Stage 1 and Stage 2 bounds.
static computeBounds_(n) {
var le = n.toString.count
var b1
var b2
if (le <= 30) {
b1 = 2000
b2 = 147396
} else if (le <= 40) {
b1 = 11000
b2 = 1873422
} else if (le <= 50) {
b1 = 50000
b2 = 12746592
} else if (le <= 60) {
b1 = 250000
b2 = 128992510
} else if (le <= 70) {
b1 = 1000000
b2 = 1045563762
} else if (le <= 80) {
b1 = 3000000
b2 = 5706890290
} else {
b1 = MaxB1
b2 = MaxB2
}
return [b1, b2]
}
// Private method to add two specified P and Q points (in Montgomery form).
// Assumes R = P - Q.
static pointAdd_(px, pz, qx, qz, rx, rz, n) {
var u = (px - pz) * (qx + qz)
var v = (px + pz) * (qx - qz)
var upv = u + v
var umv = u - v
var x = rz * upv.square
if (x >= n) x.mod(n)
var z = rx * umv.square
if (z >= n) z.mod(n)
return [x, z]
}
// Private method to double a point P (in Montgomery form).
static pointDouble_(px, pz, n, a24) {
var u = (px + pz).square
var v = (px - pz).square
var t = u - v
var x = u * v
if (x >= n) x.mod(n)
var z = t * (v + a24 * t)
if (z >= n) z.mod(n)
return [x, z]
}
// Private method to multiply a specified point P (in Montgomery form)
// by a specified scalar.
static scalarMultiply_(k, px, pz, n, a24) {
var sk = k.toString(2)
var lk = sk.count
var qx = px.copy()
var qz = pz.copy()
var pr = pointDouble_(px, pz, n, a24)
var rx = pr[0]
var rz = pr[1]
for (i in 1...lk) {
if (sk[i] == "1") {
pr = pointAdd_(rx, rz, qx, qz, px, pz, n)
qx = pr[0]
qz = pr[1]
pr = pointDouble_(rx, rz, n, a24)
rx = pr[0]
rz = pr[1]
} else {
pr = pointAdd_(qx, qz, rx, rz, px, pz, n)
rx = pr[0]
rz = pr[1]
pr = pointDouble_(qx, qz, n, a24)
qx = pr[0]
qz = pr[1]
}
}
return [qx, qz]
}
// Private method which sieves for primes up to and including 'n'
// and returns a list of the primes found.
static primeSieve_(n) {
if (n < 2) return []
var primes = [2]
var k = ((n-3)/2).floor + 1
var marked = List.filled(k, true)
var limit = ((n.sqrt.floor - 3)/2).floor + 1
limit = limit.max(0)
for (i in 0...limit) {
if (marked[i]) {
var p = 2*i + 3
var s = ((p*p - 3)/2).floor
var j = s
while (j < k) {
marked[j] = false
j = j + p
}
}
}
for (i in 0...k) {
if (marked[i]) primes.add(2*i + 3)
}
return primes
}
// Private method which uses binary seach to find the index
// of the first occurrence of an element in a list.
// Returns list.count (not -1) if the element is not present.
static findFirst_(a, t) {
var n = a.count
var l = 0
var r = n - 1
while (l != r) {
var m = ((l + r)/2).ceil
if (a[m] > t) {
r = m - 1
} else {
l = m
}
}
if (a[l] == t) return l
return n
}
// Private worker method for Lenstra's two-state ECM algorithm.
static ecm_(n) {
var pr = computeBounds_(n)
var b1 = pr[0]
var b2 = pr[1]
var dd = b2.sqrt.floor
var beta = List.filled(dd + 1, null)
for (i in 0...beta.count) beta[i] = Mpz.new()
var s = List.filled(2*dd + 2, null)
for (i in 0...s.count) s[i] = Mpz.new()
// stage 1 and stage 2 precomputations
var curves = 0
var logB1 = b1.log
var primes = primeSieve_(b2 - 1)
var numPrimes = primes.count
var idxB1 = findFirst_(primes, b1)
if (idxB1 == -1) idxB1 = numPrimes
// compute a B1-powersmooth integer 'k'
var k = Mpz.one
for (i in 0...idxB1) {
var p = primes[i]
var bp = Mpz.from(p)
var t = (logB1 / p.log).floor
var bt = bp.pow(t)
k.mul(bt)
}
var g = Mpz.one
while ((g == 1 || g == n) && curves <= MaxCurves) {
curves = curves + 1
var st = Rand.int(6, MaxRnd)
var sigma = Mpz.from(st)
// generate a new random curve in Montgomery form with Suyama's parameterization
var u = (sigma * sigma).sub(5).mod(n)
var v = (sigma * 4).mod(n)
var vmu = v - u
var a = vmu.cube
var t = (u * 3).add(v)
a.mul(t)
t = (u * u).mul(u).mul(v).mul(4)
a.div(t).sub(2).mod(n)
var a24 = (a + 2).div(4)
// stage 1
var px = (u * u).mul(u)
t = (v * v).mul(v)
px.div(t).mod(n)
var pz = Mpz.one
pr = scalarMultiply_(k, px, pz, n, a24)
var qx = pr[0]
var qz = pr[1]
g.gcd(n, qz)
// if stage 1 is successful, return a non-trivial factor else
// move on to stage 2
if (g != 1 && g != n) return g
// stage 2
pr = pointDouble_(qx, qz, n, a24)
s[1] = pr[0]
s[2] = pr[1]
pr = pointDouble_(s[1], s[2], n, a24)
s[3] = pr[0]
s[4] = pr[1]
beta[1].mul(s[1], s[2]).mod(n)
beta[2].mul(s[3], s[4]).mod(n)
for (d in 3..dd) {
var d2 = 2 * d
pr = pointAdd_(s[d2-3], s[d2-2], s[1], s[2], s[d2-5], s[d2-4], n)
s[d2-1] = pr[0]
s[d2] = pr[1]
beta[d].mul(s[d2-1], s[d2]).mod(n)
}
g.set(1)
var b = Mpz.from(b1 - 1)
pr = scalarMultiply_(b, qx, qz, n, a24)
var rx = pr[0]
var rz = pr[1]
t = Mpz.from(dd) * 2
t.sub(b, t)
pr = scalarMultiply_(t, qx, qz, n, a24)
var tx = pr[0]
var tz = pr[1]
var q = idxB1
var step = 2 * dd
var r = b1 - 1
while (r < b2) {
var alpha = (rx * rz).mod(n)
var limit = r + step
while (q < numPrimes && primes[q] <= limit) {
var d = ((primes[q] - r) / 2).truncate
var t = rx - s[2 * d - 1]
var f = (rz + s[2 * d]).mul(t).sub(alpha).add(beta[d])
g.mul(f).mod(n)
q = q + 1
}
var trx = rx.copy()
var trz = rz.copy()
pr = pointAdd_(rx, rz, s[2*dd-1], s[2*dd], tx, tz, n)
rx = pr[0]
rz = pr[1]
tx.set(trx)
tz.set(trz)
r = r + step
}
g.gcd(n, g)
}
// no non-trivial factor found, return zero
if (curves > maxCurves) return Mpz.zero
return g
}
// Returns a factor (Mpz) of 'm' (an Mpz object or an integral Num) using the
// ECM algorithm. Returns Mpz.zero in the event of failure.
static ecm(m) {
if (m < 2) return Mpz.zero
if (m is Num) m = Mpz.from(m)
if (m.probPrime(15) > 0) return m.copy()
if (m.isSquare) return m.copy().sqrt
for (p in primeSieve_(1e5)) {
if (m.isDivisibleUi(p)) return Mpz.from(p)
}
return ecm_(m)
}
// Private worker method (recursive) to obtain the prime factors of a number (Mpz).
static primeFactors_(m, trialDivs) {
if (m.probPrime(15) > 0) return [m.copy()]
var n = m.copy()
var factors = []
var threshold = 1e11 // from which using EM may be advantageous
if (trialDivs) {
for (p in primeSieve_(1e5)) {
while (n.isDivisibleUi(p)) {
factors.add(Mpz.from(p))
n.div(p)
}
}
}
while (n > 1) {
if (n.probPrime(15) > 0) {
factors.add(n)
break
}
if (n >= threshold) {
var d = ecm_(n)
if (d != 0) {
factors.addAll(primeFactors_(d, false))
n.div(d)
} else {
factors.addAll(Mpz.primeFactors_(n, false))
break
}
} else {
factors.addAll(Mpz.primeFactorsWheel_(n))
break
}
}
factors.sort()
return factors
}
// Returns a list of the primes factors (Mpz) of 'm' (an Mpz object or an integral Num)
// using ECM or a simpler method as appropriate.
static primeFactors(m) {
if (m < 2) return []
if (m is Num) m = Mpz.from(m)
return primeFactors_(m, true)
}
}