The Monte Carlo method is a computational simulation scheme, originally introduced by Stanisław Ulam, that solves a wide variety of problems arising in chemistry, economics, finance, mathematics, physics, etc. In this note, we discuss an application of the Monte Carlo method in factoring of integers. It is also called rho method and was introduced by J. M. Pollack.
The first step is to choose an easily evaluated map from $\mathbb{Z}_n$ to itself. A popular choice is $f(x)=x^2+1$. Next, one chooses an initial value $x=x_0$, and then computes the successive iterates of $f$: $x_1=f(x_0)$, $x_2=f(x_1)$, $x_3=f(x_2)$, etc. i.e. $x_{j+1}=f(x_j)$, $j=0,1,2,\cdots$. Compare the $x_j$’s, hoping to find two which are different modulo $n$ but the same modulo some divisor of $n$. Once we find such $x_i$, $k_k$, we have $(x_k-x_j,n)$ equal to a proper divisor of $n$
Example. Let us factor $91$ by choosing $f(x)=x^2+1$, $x_0=1$.
\begin{align*} x_1&=f(x_0)=2\\ x_2&=f(2)=5\\ x_3&=f(5)=26\\ &\vdots \end{align*}
Since $(x_3-x_2,n)=(21,91)=7$, 7 is a factor.
The method works by successively computing $x_k=f(x_{k-1})$ and comparing $x_k$ with the earlier $x_j$ until we find a pair satisfying $(x_k-x_j,n)=r>1$. But as $k$ gets larger, it becomes more time consuming to compute $(x_k-x_j,n)$ for each $j<k$. Note that once there is a pair $(k_0,j_0)$ such that $x_{k_0}\equiv x_{j_0} \mod r|n$, we have the same relation $x_k\equiv x_j\mod r$ for any pair $(j,k)$ such that $k-j=k_0-j_0$: Set $k=k_0+m$ and $j=j_0+m$, and apply $f$ to both sides of the congruence $x_{k_0}\equiv x_{j_0}\mod r$ repeatedly $m$ times.
The previous algorithm can be modified so that we need to calculate the gcd only once for each $k$. This significantly reduces the required computational burden. Here is the modified algorithm.
We successively compute the $x_k$. For each $k$, we proceed as follows. Suppose $k$ is an $(h+1)$-bit integer, i.e. $2^h\leq k<2^{h+1}$. Let $j$ be the largest $h$-bit integer: $j=2^h-1$. We compare $x_k$ with this particular $x_j$, i.e. compute $(x_k-x_j,n)$. If this gcd gives a nontrivial factor of $n$, stop. Otherwise continue on to $k+1$.
Example. $n=91$, $f(x)=x^2+1$, $x_0=1$.
\begin{align*} x_1&=f(1)=2\\ x_2&=f(2)=5;\ (x_2-x_1,n)=(5-2,91)=1\\ x_3&=f(5)=26;\ (x_3-x_1,n)=(24,91)=1\\ x_4&=f(26)=26^2+1\equiv 40\mod 91;\ (x_4-x_3,n)=(14,91)=7 \end{align*}
Example. Factor 4087 using $f(x)=x^2+x+1$ and $x_0=2$.
\begin{align*} x_1&=f(2)=7;\ (x_1-x_0,n)=(7-2,4087)=1\\ x_2&=f(7)=57;\ (x_2-x_1,n)=(57-7,4087)=1\\ x_3&=f(57)=3307;\ (x_3-x_1,n)=(3307-7,4087)=1\\ x_4&=f(3307)\equiv\mod 4087;\ (x_4-x_3,n)=(2745-3307,4087)=1\\ x_5&=f(2745)\equiv 1343\mod 4087; (x_5-x_3,n)=(1343-3307,4087)=1\\ x_6&=f(1343)\equiv 2626\mod 4087;\ (x_6-x_3, n)=(2626-3307,4087)=1\\ x_7&=f(2626)\equiv 3734\mod 4087;\ (x_7-x_3,n)=(3734-3307,4087)=61 \end{align*}
Hence, $61$ is a factor of $4087$ and $4087=61\cdot 67$.