# Elliptic curve method

The **elliptic curve method** (sometimes called **Lenstra elliptic curve factorization**, commonly abbreviated as **ECM**) is a factorization method which computes a large multiple of a point on a random elliptic curve modulo the number to be factored. It is currently the best algorithm known, among those whose complexity depends mainly on the size of the factor found.

This complexity can be represented by:

- [math]O(\exp{\sqrt{(\log p \,\log \log p)(1+O(1)}})[/math]

where p is the smallest factor.

## Contents

## Maths

### Simple Explanation

This is a simple explanation of the Elliptic Curve Method. It is not a mathematical proof of the method, nor is it an explanation of why it works; it is simply an explanation of what the method is.

The Elliptic Curve Method is a way of determining factors of a composite number. This method cannot be used when it is not known in advance that the number is composite, so it cannot be used as a primality test.

First we need to understand some simple mathematical terms. A prime number is a number greater than 1 that is only divisible by itself and 1. A composite number is a number that has divisors that are neither itself, nor 1. A highly composite number is a number that has lots and lots of divisors. We can make highly composite numbers simply by multiplying lots of different prime numbers together. So that 2 x 3 x 5 x 7 x 11 x 13 etc will be a highly composite number. But that is only the first six prime numbers multiplied together. If we multiply the first 2000 prime numbers together we get a huge number that has lots and lots of divisors. The significance of this to the Elliptic Curve Method is that a huge amount of other numbers will have factors in common with our highly composite number.

Modular arithmetic can be thought of as a system of arithmetic in which we only concern ourselves with the remainder after dividing by a specific number called the modulus. So that if we take 17 as our modulus, when we multiply 47 x 23, instead of saying that the answer is 1081, we say it leaves a remainder of 10 after division by the modulus. We can then abbreviate this to 47 x 23 = 10 (mod 17). The significance of this to the Elliptic Curve Method is that all of our maths is done with the number we are trying to factor as the modulus.

If we plug different values for *x* into an algebraic expression such as [math]y\,=\,3x^2\,+\,10x\,-\,9[/math], we get out different values for *y*. If we plug in *x* = 2, we get out *y* = 23. Plug in *x* = 3, we get out *y* = 48. If we plug in enough different values for *x* and plot the resultant *y* values on a graph, when we join the dots on the graph we get a curve. For the expression given above we would get a quadratic curve because the expression [math]y\,=\,3x^2\,+\,10x\,-\,9[/math] is a quadratic expression (because the highest power of *x* is 2).

If we change the expression to [math]y^2\,\equiv \,x^3\,+\,ax\,+\,b[/math] (where *a*, *b*, *x* and *y* are integers like 3 or -12, and the sign [math]\equiv[/math] means that the difference between the left hand side and the right hand side is a multiple of the number to be factored) we get a different kind of curve called an elliptic curve. For example, if we want to factor the number 77, the point (*x*, *y*) = (12, 25) is located in the elliptic curve [math]y^2\,\equiv \,x^3\,+\,17x\,+\,2[/math] because the difference is: [math]25^2-(12^3+17*12+2) = -1309[/math], a multiple of 77.

So now we have an algebraic expression that generates an elliptic curve. We have a highly composite number that has loads and loads of divisors, and we have some arithmetic being done modulo the number being factored. How do these three go together ?

We start by calculating a point on the curve. (Naturally the computer programmes that do this do not actually draw graphs, they simply calculate the *x* and *y* values for a particular point on the curve). Now, by doing some very beautiful mathematics explained below, we can calculate the position of other points on the same curve, the *multiples* of the original point. This process is called "running a curve" as in "I have run 400 curves" and in terms of the computer program doing this it is one iteration of the algorithm. So we compute the *product* of the first point times the highly composite number, and then calculate the greatest common divisor of the coordinate *x* of the found point and the number to be factored. If the greatest common divisor is greater than 1 and less than the number to be factored, congratulations! we have found a proper factor. Otherwise we have to run more curves, using a different starting point.

As described above it probably sounds a bit of a hit-or-miss affair, but there are controls built into the programme. By carefully determining the starting values for the programme it is possible to limit our search to factors of no more than some size measured in bits (a bit is the smallest unit of computer memory). After running the requisite number of curves (detailed on a table on this page) at one particular bit level, we can then start to look for factors at a higher bit level. Even at the higher bit level it is still possible to find factors of the smaller size. You can think of this as being a tall man and his short wife with an umbrella. If the wife is holding the umbrella, when we look under the umbrella we will only find her. But if her husband is holding the umbrella when we look underneath it we may find either of them, or even both. By running sufficient curves at each bit level we can eliminate the possibility of their being any factors hiding under the umbrella.

### More detailed explanation

#### Introduction

The elliptic curve stated above (Weierstrass form) is not the best in terms of computations. The Montgomery form is used instead:

- [math]By^2z \equiv x(x^2+Axz+z^2)[/math] (mod [math]N[/math])
**(1)**

which is an elliptic curve when [math]B\neq 0[/math] and [math]A\neq \pm 2[/math].

Notice that the curve (really a surface in the 3-dimensional space) above always contains the point [math]x=0[/math], [math]y=1[/math], [math]z=0[/math]. This point is known as the neutral element **O**. Also the equation does not change by multiplying all coordinates by the same value [math]k\neq 0[/math].

In this paragraph let's suppose that we work modulo a prime [math]p[/math] (this will be an a-priori unknown prime factor of [math]N[/math]). Given the set of points for which **(1)** holds, we can define the addition operation. So we can write **P** + **Q** = **R** where **P**, **Q** and **R** are points in this set. We can write **P** + **P** = **S** as 2**P** = **S**. Given a point **P**, we could compute 2**P**, 3**P**, 4**P**, and so on. For some value [math]g[/math]. we will have [math]g[/math]**P** = **O**. This value [math]g[/math] is known as the group order and depends on the values of [math]p[/math], [math]A[/math] and [math]B[/math]. Its value is near [math]p[/math], or more exactly between [math]p-2sqrt{p}[/math] and [math]p+2sqrt{p}[/math]. It should be obvious that [math]2g[/math]**P** = **O**, [math]3g[/math]**P** = **O**, and so on. For some values of **P**, it is possible to find [math]f[/math]**P** = **O** where [math]f[/math] is a divisor of the group order, but this does not change the factorization method (in fact it works better).

It is interesting to note that the elliptic curve method does not use the [math]y[/math] coordinate, thus making additions and duplications faster.

#### Step 1

Now if we work modulo the number to be factored, [math]N[/math], if we start with the point **P**, and somehow we know a multiple of the group order, say [math]kg[/math], we will find that the coordinates [math]x[/math] and [math]z[/math] of [math]kg[/math]**P** will be multiples of the unknown prime factor. So we compute the greatest common divisor between the coordinate [math]x[/math] and the number to be factored [math]N[/math] in order to reveal the prime. Notice that it is possible that [math]kg[/math] is also multiple of the group factor for another prime factor, especially if [math]kg[/math] is very high. In that case the GCD will be a composite factor of [math]N[/math]. In the worst case the GCD will be [math]N[/math].

How do we find a multiple of the group order that is not known in advance? Using highly composite numbers. Given a bound B_{1}, we multiply the original point **P** by all prime numbers less than B_{1} (each prime number is raised to a power such that the result is also about B_{1}). If the group order is a composite number whose factors are all less than B_{1} we are done. Otherwise we have to select another curve and another initial point **P** or continue with *step 2*.

#### Step 2

The step 2 is useful when the group order has a prime factor [math]q[/math] between the bounds B_{1} and B_{2} and all the other prime factors less than B_{1}. Since this is almost always the case when ECM is successful, this step 2 speeds a lot the elliptic curve method.

Let **Q** be the point computed in the *step 1* (the previous paragraph), so [math]q[/math]**Q** = **O** (considering arithmetic [math]\pmod p[/math]). Let [math]C[/math] be the multiple of 6 before B_{1}. The idea is to compute [math]C[/math]**Q** and 6**Q**. Then using the addition formula we can compute ([math]C+6[/math])**Q**, ([math]C+12[/math])**Q**, ([math]C+18[/math])**Q**, ..., up to B_{2}**Q**. If [math]q = C+1+6k[/math] the point computed [math](C+6k)[/math]**Q** will be equal to **-Q** and if [math]q = C-1+6k[/math] the point computed will be equal to **Q**. Since the coordinates [math]x[/math] of the points **Q** and **-Q** are congruent mod p, this coordinate will be congruent to the coordinate [math]x[/math] of [math](C+6k)[/math]**Q** (mod p), so we just multiply all [math]X_{(C+6r)Q} - X_Q \,\pmod N[/math]. This value will be a multiple of [math]p[/math], so taking the GCD of the product and [math]N[/math] the factor [math]p[/math] will be revealed.

An improvement to step 2 can be done if there is enough memory to store intermediate values. As a simple example let's consider the modulus 30. The numbers coprime to it are: 1, -1, 7, -7, 11, -11, 13 and -13. So we can precompute a table of four values: **Q**, 7**Q**, 11**Q** and 13**Q** (of course the first one does not require any calculation), and the values 30**Q** and [math]C[/math]**Q** as done in the previous paragraph (in this case [math]C[/math] is the multiple of 30 before B_{1}). Then we multiply all [math](X_{(C+30k)Q} - X_Q)[/math] [math](X_{(C+30k)Q} - X_{7Q})[/math][math](X_{(C+30k)Q} - X_{11Q})[/math][math](X_{(C+30k)Q} - X_{13Q})\,\pmod N[/math] and take the GCD of the product and [math]N[/math] in order to reveal the factor.

The modulus 30 can be replaced by 210, or 2310 if there is enough memory to hold all values of [math]t[/math]**Q** where the number [math]t[/math] is positive, less than half the modulus and coprime to it.

Another optimization can be done when the modular multiplications are expensive: when both [math]C+dk+t[/math] and [math]C+dk-t[/math] are composite, do not multiply the value [math]X_{(C+dk)Q} - X_{tQ}[/math].

For high values of B_{2} another method, known as the birthday paradox continuation, is faster.

### Formulas for addition and duplication

Since the algorithm does not use the [math]y[/math] coordinate we represent a point **P** as (P_{x} : : P_{z}). Notice that the second coordinate is not shown deliberately because it is not used.

If **P** = -**Q** then P_{x} = Q_{x} and P_{z} = Q_{z} (they only differ in the [math]y[/math] coordinate).

Let **P** = (P_{x} : : P_{z}) and **Q** = (Q_{x} : : Q_{z}) where P_{x} / P_{z} is not equal to Q_{x} / Q_{z} (otherwise they are the same point or negatives). Also let **P** + **Q** = (X_{+} : : Z_{+}) and **P** - **Q** = (X_{-} : : Z_{-}).

[math]\large\frac {X_+}{Z_+} = \frac {Z_-(U + V)^2}{X_-(U - V)^2}[/math]

where [math]U = (P_x - P_z)(Q_x + Q_z)[/math] and [math]V = (P_x + P_z)(Q_x - Q_z)[/math]

Since the previous formula does not work when P_{x} / P_{z} = Q_{x} / Q_{z} we need another one to compute 2**P**.

Let **P** = (P_{x} : : P_{z}) and **Q** = 2**P** = (Q_{x} : : Q_{z}).

We get:

[math]\frac {Q_x}{Q_z} = \frac {(P_x + P_z)^2(P_x - P_z)^2}{T((P_x - P_z)^2 + ST)}[/math]

where: [math]S = \frac {A+2}4[/math] which can be precomputed when the elliptic curve is selected and

[math]T = (P_x + P_z)^2 - (P_x - P_z)^2[/math].

All the arithmetic above must be done modulo [math]N[/math].

Only modular additions, subtractions and multiplications are needed. No modular inversions are required (they are very expensive in execution time).

In order to compute [math]S[/math] without modular inversions, we get the first number of the sequence [math]A+2[/math],

[math]A+2+N[/math], [math]A+2+2N[/math] and [math]A+2+3N[/math] which is multiple of 4 (which can be seen by inspecting the two least significant bits). Then just divide by 4 that number.

### Multiplication

In order to perform multiplication of a point **P** by a natural number [math]M[/math] (which is the main operation in step 1), there are two approaches using the formulas shown above.

Montgomery's PRAC is an almost optimal algorithm, but this is too complicated to explain it here.

Another algorithm is the same used in exponentiation. We represent the value [math]M[/math] in binary form. Then erase the first "1" (all numbers when represented in binary start with the digit "1"). Then for every "1" write the string "DA" and for every "0" write "D".

For example for the number 21 which is 10101 when written in binary, you will have the string DDADDA. The letter D means duplication and A means addition.

So if we have to compute 21**P** using the recipe above we will get the sequence: 2**P**, 4**P**, 5**P**, 10**P**, 20**P**, 21**P**.

A problem with this approach is that the formula for addition shown above for **Q** + **R** needs the coordinates of **Q** - **R**. This can be overcomed by using another point. When we compute k**P** we will be also computing (k+1)**P**.

Then if we start with the points (k**P**, (k+1)**P**), when the bit equals zero we compute (2k**P**, (2k+1)**P**) and when the bit equals one we compute ((2k+1)**P**, 2(k+1)**P**). In both cases we need one duplication and one addition requiring only 10 modular multiplications.

The PRAC algorithm is 30% faster in average than the algorithm shown above to multiply a point by a natural number.

### Selecting curve and starting point

Using the Montgomery form of elliptic curve, the group order is always multiple of four. Suyama found a method where the group order is multiple of 12, thus increasing the probability that this number has only small prime factors.

First we have to select a number sigma ([math]\sigma [/math]) which is not equal to 0, -1, 1 or 5. Then with the following formulas we can find the coordinates [math]x[/math] and [math]z[/math] and the elliptic curve parameter [math]A[/math]. All computations must be done modulo [math]N[/math].

- [math]u = \sigma ^2 - 5[/math]
- [math]v = 4\sigma[/math]
- [math]x = u^3[/math]
- [math]z = v^3[/math]
- [math]A = \frac{(v-u)^3 (3u+v)}{4u^3v} - 2[/math]

## Implementations

- GMP-ECM
- Prime95
- A Web application, written by Dario Alpern, to factorize using ECM and SIQS.

## Choosing the best parameters for ECM

It's still an open question how to choose the best parameters for ECM. Although, a few conventions have been made, including the B1 parameter for factors up to 70 digits (as you can see, that is a good limit because ECM on 70-digit factors is almost impossible with current hardware). Many programs today use the default B2=100*B1:

Digits | B1(B2=B1*100) | Curves |
---|---|---|

15 | 2,000 | 25 |

20 | 11,000 | 90 |

25 | 50,000 | 300 |

30 | 250,000 | 700 |

35 | 1,000,000 | 1,800 |

40 | 3,000,000 | 5,100 |

45 | 11,000,000 | 10,600 |

50 | 43,000,000 | 19,300 |

55 | 110,000,000 | 49,000 |

60 | 260,000,000 | 124,000 |

65 | 850,000,000 | 210,000 |

70 | 2,900,000,000 | 340,000 |

GMP-ECM also has default B2 values to use, which are commonly taken by most factorers. Some say that the B2 value doesn't affect the chances of finding a factor much, but as a rule of thumb it's common to spend as many time in stage 2 as in stage 1, and with the progress made on algorithms, this value should be much higher than B1*100. As we can see, using a higher B2 does affect the chances of finding a factor, according to GMP-ECM's probability function. A red square means N/A:

Digits | Value of B1 | GMP-ECM B2 default | Curves |
---|---|---|---|

15 | 2,000 | 147,396 | |

20 | 11,000 | 1,873,422 | 86 |

25 | 50,000 | 12,746,592 | 214 |

30 | 250,000 | 128,992,510 | 430 |

35 | 1,000,000 | 1,045,563,762 | 910 |

40 | 3,000,000 | 5,706,890,290 | 2,351 |

45 | 11,000,000 | 35,133,391,030 | 4,482 |

50 | 43,000,000 | 240,490,660,426 | 7,557 |

55 | 110,000,000 | 776,278,396,540 | 17,884 |

60 | 260,000,000 | 3,178,559,884,516 | 42,057 |

65 | 850,000,000 | 15,892,628,251,516 | 69,471 |

70 | 2,900,000,000 | 105,101,237,217,912 | 102,212 |

75 | 7,600,000,000 | 425,332,376,469,022 | 188,056 |

80 | 25,000,000,000 | 2,551,982,328,195,322 | 265,557 |

## Ongoing factorization effort with ECM

There is an ongoing effort by GIMPS to factorize composite Mersenne and Fermat numbers using ECM. As of May 2018, the lowest composite Mersenne number without any known factors is 2^{1277}-1. [1]

The largest factor found using ECM so far has 83 decimal digits and was discovered on 2013-09-07 by R. Propper.[2]