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
<lang ecmascript>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) }
}</lang>