离散对数的波拉德ρ算法

**离散对数的波拉德ρ算法**是1978年所发明解决离散对数问题的算法。

算法的目标是求 \gamma 使得 \alpha ^ \gamma = \beta，其中 \beta 属于一个由 \alpha 生成的 循环群 G。该算法寻找 a, b,  A, B 使得 \alpha^a \beta^b = \alpha^A \beta^B。若他们基于的群是一个 n 阶的循环群，则 \gamma 是方程 (B-b)\gamma = (a-A) \pmod{n} 的其中一个解。

为求得 a, b,  A, B， 该算法使用 Floyd判圈算法 在数列 x_i = \alpha^{a_i} \beta^{b_i} 中寻找一个环。 假设映射 f: x_i \mapsto x_{i+1} 是近似于随机的，则有可能在约 \sqrt{\frac{\pi n}{2}} 步后发现一个环。可使用一下规则来生成一个此类映射：将 G 分割为三个不相交的子集 S_0，S_1
，S_2 ，且其所含元素数量大致相等，如果 x_i \in S_0 则将 a 和 b 加倍； 如果 x_i \in S_1 则将 a 自增； 如果 x_i \in S_2 则将 b 自增。

## 算法
使  是一个  阶的 循环群, 且有 \alpha, \beta\in G, 以及一个分割 G = S_0\cup S_1\cup S_2, 定义映射 f:G\to G 

:
f(x) = \begin{cases}
\beta x & x\in S_0\\
x^2 & x\in S_1\\
\alpha x & x\in S_2
\end{cases}


并据以下方式定义映射 g:G\times\mathbb{Z}\to\mathbb{Z} 和 h:G\times\mathbb{Z}\to\mathbb{Z} 

:\begin{align}
g(x,n) &= \begin{cases}
n & x\in S_0\\
2n \pmod p & x\in S_1\\
n+1 \pmod p & x\in S_2
\end{cases}
\\
h(x,n) &= \begin{cases}
n+1 \pmod p & x\in S_0\\
2n \pmod p & x\in S_1\\
n & x\in S_2
\end{cases}
\end{align}

  **输入** *a*: a 是 *G* 的生成元, *b*: *G* 的一个元素
  **输出** 整数 *x* 使得 *a<sup>x</sup>* = *b*, 或者失败
 
  初始化 *a*0 &larr; 0, *b*0 &larr; 0, *x*0 &larr; 1 &isin; *G*, 
 
  *i* &larr; 1
  **loop**
      *x<sub>i</sub>* &larr; *f*(*x**i*-1), 
      *a<sub>i</sub>* &larr; *g*(*x**i*-1, *a**i*-1), 
      *b<sub>i</sub>* &larr; *h*(*x**i*-1, *b**i*-1)
 
      *x*2*i* &larr; *f*(*f*(*x*2*i*-2)), 
      *a*2*i* &larr; *g*(*f*(*x*2*i*-2), *g*(*x*2*i*-2, *a*2*i*-2)), 
      *b*2*i* &larr; *h*(*f*(*x*2*i*-2), *h*(*x*2*i*-2, *b*2*i*-2))
 
      **if** *x<sub>i</sub>* = *x*2*i* **then**
          *r* &larr; *b<sub>i</sub>* - *b*2*i*
          **if** r = 0 **return failure**
          *x* &larr; *r*−1(*a*2*i* - *a<sub>i</sub>*) mod *p*
          **return** *x*
      **else** # *x<sub>i</sub>* &ne; *x*2*i*
          *i* &larr; *i*+1, 
          **break loop**
      **end if**
   **end loop**

## 举例
例如一个由 2 模 N=1019 生成的群（群的阶是n=1018，2是生成元，生成群的元素模1019同余)。这个算法可由以下 C++ 程序实现。


```cpp
#include <stdio.h>
 
 const int n = 1018, N = n + 1;  /* N = 1019 -- 素数       */
 const int alpha = 2;            /* 生成元                 */
 const int beta = 5;             /* 2^{10} = 1024 = 5 (N) */
 
 void new_xab( int& x, int& a, int& b ) {
   switch( x%3 ) {
   case 0: x = x*x     % N;  a =  a*2  % n;  b =  b*2  % n;  break;
   case 1: x = x*alpha % N;  a = (a+1) % n;                  break;
   case 2: x = x*beta  % N;                  b = (b+1) % n;  break;
   }
 }
 
 int main(void) {
   int x=1, a=0, b=0;
   int X=x, A=a, B=b;
   for(int i = 1; i < n; ++i ) {
     new_xab( x, a, b );
     new_xab( X, A, B );
     new_xab( X, A, B );
     printf( "%3d  %4d %3d %3d  %4d %3d %3d\n", i, x, a, b, X, A, B );
     if( x == X ) break;
   }
   return 0;
 }
```


结果如下 (已截断):

  i     x   a   b     X   A   B
 ------------------------------
  1     2   1   0    10   1   1
  2    10   1   1   100   2   2
  3    20   2   1  1000   3   3
  4   100   2   2   425   8   6
  5   200   3   2   436  16  14
  6  1000   3   3   284  17  15
  7   981   4   3   986  17  17
  8   425   8   6   194  17  19
 ..............................
 48   224 680 376    86 299 412
 49   101 680 377   860 300 413
 50   505 680 378   101 300 415
 51  1010 681 378  1010 301 416

可见 2^{681} 5^{378} = 1010 = 2^{301} 5^{416} \pmod{1019} 以及 (416-378)\gamma = 681-301 \pmod{1018}。
正如预期，其中 \gamma_1=10 是一个解。由于 n=1018 不是素数，因此存在另一个解 \gamma_2=519，使得 2^{519} = 1014 = -5\pmod{1019} 成立。

## 复杂度
时间复杂度近似于 \mathcal{O}(\sqrt{n})。如果配合使用 ，则整体时间复杂度近似于 \mathcal{O}(\sqrt{p})，其中 p 是 n 的最大质因数。

## 参考文献
*
*