质因数分解

Pollard$\rho$质因数分解法

假设$N$是一个能分解为$p*q$的数,那么对$N$的质因数分解就可以转换成分别对$p$和$q$求解质因数,现在的问题是如何快速的找到$N$的一个因数$p$.
传统的试除法需要一个一个找到因子,对于$N$来说,如果$N$只有$p$和$q$两个非平凡因子,那么找到它们的概率为$\frac{2}{N-2}$,当`$N$非常大的时候,那么找到因子的概率显然很小,而Pollard$\rho$分解法的主要思想就是提高找到因子的概率.

生日悖论

问题:在$n$个人中至少两个人同一天生日的概率.
当$n>365$时,此时两人同一天生日的概率为1.
当$n\leq365$时,令$p(n)$表示$n$个人中至少两个人同一天生日的概率,那么$\overline p(n)$表示每个人的生日都不相同的概率.
那么

$$\overline p(n)=1\cdot(1-\frac{1}{365})\cdot(1-\frac{2}{365})\cdots(1-\frac{n-1}{365}) =\frac{365!}{365^n(365-n)!}$$
所以
$$p(n)=1-\overline p(n)$$
通过计算可以算出, $p(25)\approx0.5686$, $p(50)\approx0.9703$, $p(75)\approx1$,当 $n$不断增大$p(n)$的概率趋近于1.

利用生日悖论因数分解

运用Birthday Trick

选取$k$个数,判断是否存在$x_i-x_j$能够整除$N$,即$x_i-x_j=p$或$x_i-x_j=q$

随机函数

对于过大的整数,显然不能将所有的$x_i$都存下来进行比较,所以可以通过比较通过随机函数生成的序列中两个相邻的随机数

GCD优化

不检测$x_i-x_j$整除$N$,检测$gcd(x_i-x_j,N)>1$.
在实际过程中,$N$的因子可能不止只有$p$和$q$,如果找到了因子$p$与$q$的话,那么$p,2p,3p,\cdots,(q-1)*p$,$q,2q,3q,\cdots,(p-1)*q$同样满足为$N$的因子,所以只要判断$x_i$与$x_j$的公约数不为1即可.

判环

尽管引入了随机函数,但是这里的随机函数是伪随机函数,生成的随机数可能存在循环,因此要在过程中加入判断环是否存在的步骤,常用的floyd判环法:设置两个变量,一个变量通过随机函数变换一次,另一个变换两次,那么如果存在环的话,最后必定相遇.

代码

uint64_t gcd(uint64_t a, uint64_t b)
{
    return a == 0 ? b : gcd(b % a, a);
}
uint64_t PollardRho(uint64_t n, uint64_t c)
{
    uint64_t x = rand() % (n - 1) + 1;
    uint64_t y = x;
    uint64_t i = 0, k = 1,d=1;

    /*
    Floyd循环判断算法对储存整个迭代序列的空间要求很高，一般实现时都使用时间换空间的办法，同时计算x_k和x_2k来进行判断。
    Brent提出了另外一种效率更高的循环判断算法：每步只计算x_k，当k是2的方幂时，令y=x_k；
    每一步都拿当前的x_k和y计算d=gcd(x_k-y,n)。
    */
    while (d==1)
    {
        ++i;
        x = (mul_mod(x, x, n) + c) % n;//随机函数选择f(x)=x^2+c c这里选1 可以选择随机数
        d = gcd(max(x,y)-min(x,y),n);
        if (1 < d&&d < n)//找到了n的一个因数
            return d;
        if (d==n)//x与y相等 此时的因数为n (可能会存在刚开始x=y的情况导致找不出n的非平凡因子，尽管概率很小，仍然存在这种可能，所以要多次找)
            return n;
        if (k == i) //用k与i来检测环
        {
            y = x;
            k <<= 1;
        }
    }
}
void primeDecomposition(uint64_t n,map<uint64_t, uint64_t> &count)
{
    if (MillerRabin(n))//用MillerRabin法进行质数判断
    {
        if (count.find(n) == count.end())
            count[n] = 1;
        else
            count[n]++;
        return;
    }
    uint64_t p = n;
    if (n == 1) return;
    while (p==n) p = PollardRho(n, rand() % (n-1) + 1);
    primeDecomposition(p, count);
    primeDecomposition(n / p, count);
}