Large Integer Factorization Algorithm and Its Practice
By Max  20230802Large Integer Factorization
In 1977, RSA (Rivest–Shamir–Adleman), a cryptosystem based on factorizing large integers was born. It was named after its three creators, Ron Rivest, Adi Shamir, and Leonard Adleman, using their surnames' initials.
With the significant contributions of the RSA algorithm to publickey cryptography and its profound impact worldwide, the three creators were then awarded the Turing Award in 2002.
Until today, the integer factorization problem remains one of the most crucial subjects in almost all aspects of publickey cryptography security.
The security of the RSA algorithm depends on large number factorization, which is
to factorize a large integer into a product of several large prime factors.
What's interesting is that, although this is a "wellknown" challenge that has been practically demonstrated, it is still unknown whether it is a problem at $P$ class or $NP$ class.
Some typical large integer factorization algorithms targeting small factors are:
 Pollard's rho algorithm
 Pollard’s p1 algorithm
 ECM algorithm (known as Lenstra elliptic curve factorization or ellipticcurve factorization method)
In this article, we will first introduce Pollard's p1 algorithm and the ECM algorithm. Additionally, we will present a specific implementation of the ECM algorithm called GMPECM.
Pollard’s p1 Algorithm
Powersmooth: If all the prime powers $p^v$ of the number $m$ are Bsmooth, that is $p^v \le B$ , then $m!$ is Bsmooth.
Pollard’s p1 algorithm is an integer factorization algorithm using Fermat's little theorem, invented by John Pollard in 1974. According to Fermat's little theorem, $a^{K(p1)} \equiv 1 \mod p$ , if $p$ is a factor of $n$ which is to be factorized, then we can express it as $pn$ . For $p$ , if $d = a^{p1}  1 \equiv 0 \mod p$ , then $d$ is a multiple of $p$ , therefore
$p = gcd(d, n) = gcd(a^{p1}  1, n)$
For data selection, the Pollard p1 algorithm employs random numbers. If no result is obtained through random selection, it then generates $d$ from the prime numbers (less than $B$ (Bpowersmooth)) until the prime number list is exhausted.
STEPS
Input: $n$ (composite number)
Output: $n$ or a failed nontrivial factor

Select a smoothness bound $B$ .

Definition $M = \prod_{q \le B} {q^{\left \lfloor log_q{B} \right \rfloor }}$ , $q$ is a prime number.

Randomly select a number $a$ that is coprime with $n$ .

Calculation
$g = gcd(a^M  1, n)$
 If $1 < g < n$ , return $g$ .
 If $g=1$ , choose a larger value for $B$ and return to Step 2, or return failure.
 If $g=n$ , choose a smaller value for $B$ and return to Step 2, or return failure.
ECM Algorithm
In Pollard's p1 method, it requires that $n$ has a factor $p$ , and that $p1$ has relatively small factors, allowing us to find a factor of $n$ using $p=gcd(a^{p1}  1, n)$ . However, when the factors of $p−1$ are large, Pollard's p1 method is less likely to successfully factorize $n$ .
To improve this method, the ECM (Elliptic Curve Method) introduces random elliptic curves to transform the original multiplicative group into a group consisting of points on the elliptic curve.
The elliptic curve method is a numbertheoretic algorithm that utilizes the properties and arithmetic rules of points on an elliptic curve to search for factors that satisfy certain conditions.
Overview of ECM Algorithm
Assuming $p$ is a prime factor of $n$ , considering the elliptic curve $E_{a,b} \mod p$ , and according to Hasse's theorem, the order $g$ of the elliptic curve satisfies the following condition: $g(p+1)<2 \sqrt{p}$ . When the parameters $a$ and $b$ of the elliptic curve vary, the order of the elliptic curve will fall within the range $[ p + 1  2\sqrt{p}, p + 1+ 2 \sqrt{p}]$ .
When using ECM to find factors of $n$ , there are three steps:
 Select a $B$ and calculate $m = lcm(1,2,\dots ,B)$ .
 Select a random curve out of $Z/NZ$ and a random point $P$ on the curve.
 Calculate $mP$ . During the calculation, if the addition of points on the corresponding elliptic curve cannot be calculated, then a factor of $n$ has been found. Otherwise, you can select another curve for calculation, or return failure.
The principle behind the ECM algorithm finding factors of $n$ lies in the fact that when a randomly chosen elliptic curve order $g$ is Bsmooth, during the calculation of $mP$ , if we encounter $gcd(v,p) = p$ or $gcd(v,q) = q$ , it will result in a nontrivial factor of $gcd(v,n)$ .
GMPECM: An Implementation of ECM Algorithm
The GMPECM implementation makes optimizations based on the original ECM method.
The overall algorithm process in GMPECM has been optimized. In GMPECM, the calculation of factors is divided into two steps, each step selecting two values $(B1, B2)$ .
The algorithm is divided into two stages, where calculations are performed using $B1$ and $B2$ , respectively. The calculations in stage 1 are similar to the ECM algorithm.
Stage 1
First, calculate the product of the point $P$ on the elliptic curve $E_{a,b}$ :
$Q = \Pi_{\pi \le B_2} \pi^{\left \lfloor (log B_1)/(log \pi) \right \rfloor} P_0$
In Stage 1, GMPECM has optimized the selection and calculation of the elliptic curve.
Random Elliptic Curve Generation
During the factorization, when the order of the randomly generated elliptic curve varies within the range $[ p + 1  2\sqrt{p}, p + 1+ 2 \sqrt{p}]$ , and the random elliptic curve is B1smooth, it can factorize $n$ .
In GMPECM, certain methods are employed for elliptic curve generation to ensure that the order of the generated elliptic curve satisfies certain properties, for example:
 Suyama's form: Ensure that the order $g$ of the elliptic curve is divisible by 12.
 Montgomery's form: Ensure that the order $g$ of the elliptic curve is divisible by 3.
Elliptic Curve Operation
In GMPECM, the Montgomery curve is used, and the multiplication or addition of points follows the following rules:
 $P(x_P :: z_P), Q(x_Q :: z_Q) \rightarrow P + Q$
 Calculate: $x_{P+Q} = 4z_{PQ} * (x_P x_Q  z_P z_Q)^2$
 Calculate: $z_{P+Q} = 4x_{PQ} * (x_P z_Q  z_P x_Q)^2$
 $P(x_P :: z_P) \rightarrow 2P$
 Calculate: $x_{2P} = (x_p^2  z_P^2)^2$
 Calculate: $z_{2P} = (4x_P z_P)[(x_p  z_P)^2 + d(4x_P z_P)]$ , $d = (a+2) / 4$ , $a$ is a parameter in the elliptic curve.
Stage 2
In the first stage of GMPECM, the main task is to calculate the point $Q$ on the elliptic curve $E$ . If it fails, meaning that $gcd(n, zQ) = 1$ , then the search range is expanded to $[B1，B2]$ in the second stage. In this new range, there exists a prime number $\pi$ , $\pi Q = O_E \mod p$ .
The main idea of Stage 2 is the meetinthemiddle strategy:
 In Stage 2, we calculate two points: $\sigma Q = ( x_\sigma : y_\sigma)$ and $\tau Q = ( x_\tau : y_\tau)$ .
 If $\sigma Q + \tau Q = O_E \mod p$ , it implies that $x_\sigma = x_\tau \mod p$ . Therefore, calculating $gcd(x_\sigma  x_\tau, n)$ will obtain the factor $p$ .
Algorithm Complexity
The ECM method is expected to take time of $O(L(p)^{\sqrt{2} + o(1)} M(\log n))$ to find a factor $p$ of a number $n$ , where $L(p) = e^{\sqrt{\log p \log \log p}}$ .
Here, $M(\log n)$ represents the complexity of multiplication modulo $n$ .
In the actual calculation process, due to the algorithm's randomness, the time to find the factor $p$ of $n$ may vary. Additionally, the ability of the algorithm to factorize and the time required for factorization also depend on the values of $B1$ and $B2$ .
For the calculations, you can refer to the recommended values of $B1$ and $B2$ as specified in the GMPECM documentation.
Example
Take a large integer $N$ (2048 bit) as an example:y
N = 0xA2862FB145AC7CE580E0BFF3CEFF42646050F8C611D36A3026E6EFB433B5CDD9EEFBA893E3F25C23A4951BA20992680162934CBDB6876CC791738C140C6EA6EC82938F7E18C5F0760367CF20883DCAB1D6885E222BBC7A1B5D434FD9FC148E387FA9AD69CE708FDDDE3E963F9AB2AF2046A37D7DBA21BC92E6F2A631C3E7770C77C95FAD6F1DC60D9F0645AEF2CC5D03D2151E35443FE0F90F1E2EAE58489C4450C8281EDCC58DDB409C306797C616954ACABCE4881E40ABDD2689A9F7596F29CD5CB42C752DA9306E7DB87CAC8957E3CA165421CF9E1B7220759A10588B386E33FD8E762E92C9F79D50150EF5FCA5F411DE23B8DFFE47A95D48ADDCF4797565
Use the tool which uses GMPECM to factorize:
ecm c 2 25e4 1.3e8 mpzmod threads: 2 mod: 1
A2862FB145AC7CE580E0BFF3CEFF42646050F8C611D36A3026E6EFB433B5CDD9EEFBA893E3F25C23A4951BA20992680162934CBDB6876CC791738C140C6EA6EC82938F7E18C5F0760367CF20883DCAB1D6885E222BBC7A1B5D434FD9FC148E387FA9AD69CE708FDDDE3E963F9AB2AF2046A37D7DBA21BC92E6F2A631C3E7770C77C95FAD6F1DC60D9F0645AEF2CC5D03D2151E35443FE0F90F1E2EAE58489C4450C8281EDCC58DDB409C306797C616954ACABCE4881E40ABDD2689A9F7596F29CD5CB42C752DA9306E7DB87CAC8957E3CA165421CF9E1B7220759A10588B386E33FD8E762E92C9F79D50150EF5FCA5F411DE23B8DFFE47A95D48ADDCF4797565
********** Factor found in step 2: 1021791499165844943393503 21 52761ms
Then, we successfully find the prime factor 1021791499165844943393503.