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>/* 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)
   }

}</lang>