Proof of Pollard's Rho Algorithm
Algorithm
Let’s say that an element $g$ of a group $G$ has order $N = q_{1}^{r_{1}} q_{2}^{r_{2}} \cdots q_{t}^{r_{t}}$. Then, the discrete logarithm problem $g^{x} = h$ can be solved in no more than $\displaystyle O \left( \sum_{i=1}^{t} S_{q_{i}^{r_{i}}} + \log N \right)$ steps according to the following algorithm.
Step 1.
Compute $\displaystyle g_{i} : = g^{N / q_{i}^{r_{i}}}$ and $\displaystyle h_{i} := h^{N / q_{i}^{r_{i}}}$.
Step 2.
Find the solution $y_{i}$ to the discrete logarithm problem $g_{i}^{y} = h_{i}$ using Shanks’ algorithm.
Step 3.
Find $1 \le x \le N$ that satisfies $\begin{cases} x \equiv y_{1} \pmod{ q_{1}^{r_{1}} } \\ \qquad \vdots \\ x \equiv y_{t} \pmod{ q_{t}^{r_{t}} } \end{cases}$ using the Chinese Remainder Theorem.
- $S_{q_{i}^{r_{i}}}$ is the time it takes for Shanks’ algorithm, which is about $\displaystyle O \left( S_{q_{i}^{r_{i}}} \right) \approx O \left( q_{i}^{r_{i} / 2} \right) $.
Shanks’ algorithm is an attack method that can be used when the order is known. The Pollard-Rho algorithm utilizes the fact that if $p$ is smooth, then the factorization of $(p-1)$ becomes easier, making it easier to find $\displaystyle g_{i} = g^{N / q_{i}^{r_{i}}}$. The order of the created $g^{i}$ is trivially $q_{i}^{r_{i}}$, removing the constraint of using Shanks’ algorithm. This involves reducing the original problem into smaller problems, solving them individually, and then obtaining the answer using the Chinese Remainder Theorem.
Proof
Part 1. Existence
That $x$ is a solution to $g^{x} = h$ means that for all $i = 1, \cdots , t$, there exists $z_{i}$ that satisfies $x = y_{i} + q_{i}^{r_{i}} z_{i}$.
- Part 1-1. $\displaystyle {{N} \over { q_{i}^{r_{i}} } } x \equiv {{N} \over { q_{i}^{r_{i}} } } \log_{g} ( h ) \pmod{N} $
By the definition of $N$ and $g_{i} , y_{i} , h_{i}$, $$ \begin{align*} \left( g^{x} \right)^{N / q_{i}^{r_{i}} } =& \left( g^{y_{i} + q_{i}^{r_{i}} z_{i}} \right)^{N / q_{i}^{r_{i}} } \\ =& \left( g^{ N / q_{i}^{r_{i}} } \right)^{ y_{i} } \cdot g^{N z_{i} } \\ =& \left( g^{ N / q_{i}^{r_{i}} } \right)^{ y_{i} } \\ =& g_{i}^{ y_{i} } \\ =& h_{i} \\ =& h^{N / q_{i}^{r_{i}} } \end{align*} $$ is obtained. To organize, $$ \left( g^{x} \right)^{N / q_{i}^{r_{i}} } = h^{N / q_{i}^{r_{i}} } $$ and taking the logarithm gives $$ {{N} \over { q_{i}^{r_{i}} } } x \equiv {{N} \over { q_{i}^{r_{i}} } } \log_{g} ( h ) \pmod{N} $$ - Part 1-2. $\displaystyle \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} = 1$
It is trivial that $\displaystyle \gcd \left( {{N} \over { q_{1}^{r_{1}} } } , \cdots , {{N} \over { q_{t}^{r_{t}} } } \right) = 1$ holds according to the definition of $N$, and by the Extended Euclidean Theorem, $$ {{N} \over { q_{1}^{r_{1}} } } c_{1} + \cdots + {{N} \over { q_{t}^{r_{t}} } } c_{t} = 1 $$ there exists $c_{1} , \cdots , c_{t} \in \mathbb{Z}$. - Part 1-3. $x = \log_{g} (h) \pmod{N}$
From Part 1-1, $$ {{N} \over { q_{i}^{r_{i}} } } x \equiv {{N} \over { q_{i}^{r_{i}} } } \log_{g} ( h ) \pmod{N} $$ Multiplying both sides of by $c_{i}$ gives $$ {{N} \over { q_{i}^{r_{i}} } } c_{i} x \equiv {{N} \over { q_{i}^{r_{i}} } } c_{i} \log_{g} ( h ) \pmod{N} $$ Adding this from $i=1$ to $t$ gives $$ \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} x \equiv \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} \log_{g} ( h ) \pmod{N} $$ $x$ and $log_{g} (h)$ are independent of index $i$, so $$ x \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} \equiv \log_{g} ( h ) \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} \pmod{N} $$ From Part 1-2, since $\displaystyle \sum_{i=1}^{t} {{N} \over { q_{i}^{r_{i}} } } c_{i} = 1$, it follows $$ x \equiv \log_{g} (h) \pmod{N} $$
Part 2. Time Complexity
Since Step 2 is repeated for every $i = 1, \cdots , t$, it takes $\displaystyle O \left( S_{q_{i}^{r_{i}}} \right) $. Moreover, because using the Chinese Remainder Theorem in Step 3 takes time $O ( \log N )$, $\displaystyle O \left( \sum_{i=1}^{t} S_{q_{i}^{r_{i}}} + \log N \right)$
■
Code
The following is the implementation of the Pollard-Rho algorithm in R. Factorization code, Code for finding the order, Shanks’ algorithm, Exponentiation by squaring, and Code using the Chinese Remainder Theorem were used.
prime = read.table("../attachment
/cfile8.uf@25411C3C5968BBE322F0D4.txt"); prime = prime[,1]
factorize<-function(p)
{
q=p
factors<-numeric(0)
i=1; j=1
while(q!=1)
{
if(q%%prime[i]) {i=i+1}
else
{
q<-q/prime[i]
factors[j]<-prime[i]
i=1
j=j+1
}
}
return(factors)
}
order<-function(g,p,h=1) #Calculate a order of g in modulo p
{
qe<-table(factorize(p-1))
qe<-rbind(as.numeric(names(qe)),qe)
divisor<-qe[1,1]^(0:qe[2,1])
if((length(qe)/2)==1) {return(qe[1,1]^qe[2,1])}
for(i in 2:(length(qe)/2)) {divisor=c(divisor%*%t(qe[1,i]^(0:qe[2,i])))}
for(i in divisor) {if((FPM(g,i,p))%%p==1) break;}
return(i)
}
FPM<-function(base,power,mod) #It is equal to (base^power)%%mod
{
i<-0
if (power<0) {
while((base*i)%%mod != 1) {i=i+1}
base<-i
power<-(-power)}
if (power==0) {return(1)}
if (power==1) {return(base%%mod)}
n<-0
while(power>=2^n) {n=n+1}
A<-rep(1,n)
A[1]=base
for(i in 1:(n-1)) {A[i+1]=(A[i]^2)%%mod}
for(i in n:1) {
if(power>=2^(i-1)) {power=power-2^(i-1)}
else {A[i]=1} }
for(i in 2:n) {A[1]=(A[1]*A[i])%%mod}
return(A[1])
}
shanks<-function(g,h,p)
{
N<-order(g,p)
n<-1+floor(sqrt(N))
gn<-FPM(g,-n,p) #gn := g^{-n}
x<-p
List_1<-numeric(n+1)
List_1[1]=1
for(i in 1:n) {List_1[i+1]=(List_1[i]*g)%%p}
List_2<-numeric(n+1)
List_2[1]=h
for(i in 1:n) {List_2[i+1]=(List_2[i]*gn)%%p}
for(i in 0:n+1) {
for(j in 0:n+1) {
if (List_1[i]==List_2[j]) {x[length(x)+1]<-((i-1)+(j-1)*n)}
}
}
return(min(x))
}
CRA<-function(S) #Algorithm of chinese remainder theorem
{
r<-S[,1] # matrix S express below sysyem.
mod<-S[,2] # x = r[1] (mod mod[1])
n<-length(r) # x = r[2] (mod mod[2])
# x = r[3] (mod mod[3])
A<-seq(r[1],to=mod[1]*mod[2],by=mod[1])
for(i in 2:n)
{
B=seq(r[i],to=mod[1]*mod[i],by=mod[i])
r[1]=min(A[A %in% B])
mod[1]=mod[1]*mod[i]
if (i<n) {A=seq(r[1],to=mod[1]*mod[i+1],by=mod[1])}
}
return(r[1])
}
PHA<-function(g,h,p){
N<-order(g,p)
m_<-table(factorize(N))
m_<-rbind(as.numeric(names(m_)),m_)
m_<-c(data.frame(m_)[1,]^data.frame(m_)[2,])
y_<-numeric(0)
for(i in 1:length(m_)){
g_i<-FPM(g,N/m_[i],p)
h_i<-FPM(h,N/m_[i],p)
y_[i]<-shanks(g_i,h_i,p)
}
return(CRA(cbind(y_,m_)))
}
PHA(7,166,433)
FPM(7,47,433)
PHA(10,243278,746497)
FPM(10,223755,746497)
PHA(2,39183497,41022299)
FPM(2,33703314,41022299)
The result of running the above code is as follows, and it has been checked that it works correctly by verifying with exponentiation by squaring.
