浅谈质因数分解

· · 个人记录

浅谈质因数分解

->part 1：

算数基本定理：

任何一个大于1的正整数都能唯一分解为有限个质数的乘积，可写作：

N=\prod_{i=1}^m p_i^ {c_i}

其中c_i都是正整数，p_i都是质数，且满足p_1<p_2<…<p_m

->part 2：

分解方法：

结合质数判定的“试除法”和质数筛选的“Eratosthenes 筛法”,我们可以扫描 2-\sqrt N的每个数x∈Z，若x能整除N，则从N中除掉所有质因子x，同时累计除去x的个数。

上述算法保证一个合数的因子一定在扫描到这个合数之前除掉了，所以这个过程中所得到的能整除N的一定是质数。时间复杂度为O(\sqrt N)。

值得注意的一点是，若N没有被任何2-\sqrt N的整数整除，则N是质数，无需分解。

试除法的时间复杂度在O(\sqrt N)而 Pollard Rho 算法的理论实践复杂度在O(\sqrt[4]{N})，这已经大大提高了程序效率。

学习 Pollard Rho 算法前，需要掌握几个知识点：

1. Miller-Rabin 算法

因此，要对同余相关知识有所了解，详见同余相关及中国剩余定理

若p为质数，则

∀a∈Z,a^p \equiv a\pmod{p}

费马小定理的证明可以参考《算法竞赛进阶指南》，读者也可以自己尝试证明或者百度。

上述定理对于同个p、多个底数 a，若不满足上式，则p为合数。

若p为质数， 且 x^2 \equiv 1\pmod{p}， 则x \equiv 1 \pmod{p} 或 x \equiv p-1 \pmod{p}

二次探测定理就用作提高程序准确性：

我们对于一个数p且满足： a^{p-1} \equiv 1\pmod{p}[费马小定理]

那么对于2|p-1，可以得到(x^\frac{p-1}{2})^2 \equiv 1 \pmod{p}，我们在费马小定理的基础上讨论x^\frac{p-1}{2} \mod{p} 的值(以下记为K)即可：

  1. 若K ≠ 1 且 K ≠p-1 ，则p是合数 return \ false;
  2. 若K=1，则继续递归x^\frac{p-1}{2};
  3. 若K=p-1,不能利用二次探测定理继续递归，说明目前无法验证p是合数，return \ true;

(以上部分为引用内容)

这样我们多取几个x就能够提高正确率。

以上内容是为了实现Miller-Rabin 算法，这是一个高效(O(logN))的判断素数的随机算法。

昨天两个毒瘤数据还是卡了我很久，感谢 @雪中舞 巨佬的帮助

然后拿两个质数 61 和 24251 以及三个随机数去check，这样就基本不会被卡了。

然而我三个随机数还是被卡了，因为某种意外，详见下面代码；于是我写了四个随机数，这样就稳稳AC了。

    //Miller_Rabin 算法
    inline ll qmul(ll x, ll y, ll mod) { //快速乘 
        return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
    }

    inline ll qpow(ll x, ll m, ll mod) { //快速幂 
        register ll ans = 1;
        while(m) {
            if(m & 1)
                ans = qmul(ans, x, mod) % mod;
            x = qmul(x, x, mod) % mod;
            m >>= 1;
        }
        return ans;
    }

    inline bool FMLT(ll x, ll p) { //费马小定理 
        return qpow(x, p - 1, p) == 1; 
    }

    inline bool Miller_Rabin(ll x, ll p) { 
        if(!FMLT(x, p)) return 0;
        register ll k = p - 1;
        for(;!(k & 1);) {
            k >>= 1;
            ll K = qpow(x, k, p);
            if(K != 1 && K != p - 1) return 0;
            if(K == p - 1) return 1;
        }
        return 1;
    }

    inline bool check(ll p) { // 判断质数 
        if(p == 1) return 0;
        ll t1 = rand(), t2 = rand(), t3 = rand(), t4 = rand();
        if(p == 61 || p == 24251 || p == t1 || p == t2 || p == t3 || p == t4) return 1;
        /*这行写成
        if(p == 2 || p == 3 || p == t1 || p == t2 || p == t3) return 1;     
        （原三个随机数）是可以AC的，但是上面那种写法才是理论正确的，所以多加了一个随机数吧tql
        */
        bool s1 = Miller_Rabin(61, p), s2 = Miller_Rabin(24251, p);
        bool s3 = Miller_Rabin(t1, p), s4 = Miller_Rabin(t2, p);
        bool s5 = Miller_Rabin(t3, p), s6 = Miller_Rabin(t4, p);
        return s1 && s2 & s3 && s4 && s5 && s6;

    }

难道你以为这就完了？你会发现这题的小数据就跟毒瘤一样，为什么呢？因为我们rand的随机数有可能大于当前的数p，这对答案是没有贡献的，那么我们需要点优化：

    ll t1 = rand() % (p - 3) + 2

这句话把随机的范围缩小到了[2,p-1)。然后小数据就可以过了

但是，这还不是结尾，我们知道，4个随机数把合数判成质数的概率是\frac{1}{4}，我们需要缩小这个概率，那么就需要增加实验次数，于是我们可以把check和Miller \ Rabin写在一起，并增加一个“次数”repeat，在一定的时间复杂度内，我们进行多次实验。

下面是优化后的核心部分(repaet开到23貌似可以过了，当然，为了保险，可以开到50)

    inline bool Miller_Rabin(ll p, int repeat) { 

        if(p == 2 || p == 3) return true;
        if(p % 2 == 0 || p == 1) return false;
        //特判
        register ll k = p - 1;
        register int cnt = 0;
        while(!(k & 1)) ++cnt, k >>= 1;

        srand((unsigned)time(NULL));//种子
        for(register int i = 0; i < repeat; i++) { 

            register ll test = rand() % (p - 3) + 2;
            register ll x = qpow(test, k, p), y = 0;

            for(register int j = 0; j < cnt; j++) { //二次探测逆过程
                y = qmul(x, x, p);
                if(y == 1 && x != 1 && x != (p - 1)) 
                    return false;
                x = y;
            }
            if(y != 1) return false; //费马小定理 
        }
        return true;
    }

所有的改进思路感谢巨佬@雪中舞

2. 快速幂 | 快速乘(相信这个我不用说了)

    typedef long long ll;

    ll quick_mul(ll x, ll y, ll mod) {
        return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;  
    }

    ll quick_pow(ll x, ll m, ll mod) {
        ll ans = 1;
        while(m) {
            if(m & 1) 
                ans = ans * x % mod;
            x = x * x % mod;  
            m >>= 1; 
        }
        return ans;
    }

3. gcd(这个我也不用说了吧)

两种办法：更相减损术 | 辗转相除法

下面给出辗转相除法的代码：

    ll gcd(ll x,ll y) {
        ll temp;
        while(y) {
            temp = x % y;
            x = y;
            y = temp;   
        }
        return x;
    }

当然，这是最简单的gcd算法，我们还可以进行二进制优化： 详见约数相关

下面为Pollard-Rho 算法实现：

这里还要引入一个“生日悖论”，利用这个玄学的悖论，我们就可以快速质因数分解了。然后我们可以利用rand()函数进行随机，在一直找到一个不符合的数进入死循环时利用Floyd的判断环的方式，一个“正常倍速”，一个“二倍速”，两者“第一次相遇时”，说明“走完了一圈”，即进入死循环。

//参考代码 
    void Pollard_Rho(int x) 
    {
        if(test(x)) {//素数测试
            ans=max(x, ans);
            return; 
        }
        ll t1 = rand() % (x-1) + 1;
        ll t2 = (qmul(t2, t2, x)+b) % x;
        ll b = rand() % (x-1) + 1;

        for(;t1 != t2;) {//Flord判断是否进入死循环 
            ll t = gcd(abs(t1-t2), x);
            if(t != 1 && t != x) 
            {
                Pollard_Rho(t),Pollard_Rho(x/t);    
            }
            t1 = (qmul(t1, t1, x)+b) % x;
            t2 = (qmul(t2, t2, x)+b) % x;
            t2 = (qmul(t2, t2, x)+b) % x;//“两倍速” 
        }
    }

上述代码时间复杂度还没达到理想时间复杂度O(\sqrt [4]N)，主要瓶颈在于频繁的求gcd。

那么这里有个优化的玄学常数127，每取127样本才进行1次gcd，看起来是挺划算的？

下面是优化的核心部分代码(这真的是我手打的)：

    inline void Pollard_Rho(ll x) {
        if(check(x)) { //检查质数 
            ans = max(x, ans);
            return;
        }
        ll s1 = rand() % (x - 1) + 1;
        ll m = rand() % (x - 1) + 1, p = 1;
        ll s2 = (qmul(s1, s1, x) + m) % x;
        for(register ll i = 1; s1 != s2; i++) { 
        //边界：不进入死循环 即S1 != S2（Floyd判断环） 
            p = qmul(p, abs(s1 - s2), x);
            if(p == 0) {
                ll K = gcd(abs(s1 - s2), x);
                if(K != 1 && K != x) 
                    Pollard_Rho(K), Pollard_Rho(x / K);
                return; 
            }
            if(i % 127 == 0) {
                p = gcd(p, x);
                if(p != 1 && p != x) {
                    Pollard_Rho(p), Pollard_Rho(x / p);
                    return;
                }
                p = 1;
            }
            s1 = (qmul(s1, s1, x) + m) % x;
            s2 = (qmul(s2, s2, x) + m) % x;
            s2 = (qmul(s2, s2, x) + m) % x; //跑两倍 
        }
        //以下部分：最后有可能一个环长不到127 
        p = gcd(p, x);
        if(p != 1 && p != x) {
            Pollard_Rho(p), Pollard_Rho(x / p);
            return;
        }
    }

对此有一定理解的读者可以尝试完成模板题：

P4718 【模板】Pollard-Rho算法 ，这也是上一章质数筛法的必做题目，下面给出参考代码：

Pollard-Rho$ 算法 $Code: 
    #include<iostream>
    #include<cstdlib>
    #include<cstdio>
    #include<cmath>
    #include<ctime>

    using namespace std;

    typedef long long ll;
    ll ans,n;

    inline ll qmul(ll x, ll y, ll mod) {  
        return (x*y-(ll)((long double)x/mod*y)*mod+mod)%mod;
    }

    inline ll qpow(ll x, ll m, ll mod) { 
        register ll ans = 1;
        while(m) {
            if(m & 1)
                ans = qmul(ans, x, mod) % mod;
            x = qmul(x, x, mod) % mod;
            m >>= 1;
        }
        return ans;
    }

    inline bool FMLT(ll x, ll p) { 
        return qpow(x, p - 1, p) == 1; 
    }

    inline bool Miller_Rabin(ll p, int repeat) { 
        if(p == 2 || p == 3) return true;
        if(p % 2 == 0 || p == 1) return false;
        register ll k = p - 1;
        register int cnt = 0;
        while(!(k & 1)) ++cnt, k >>= 1;
        srand((unsigned)time(NULL));
        for(register int i = 0; i < repeat; i++) {
            register ll test = rand() % (p - 3) + 2;
            register ll x = qpow(test, k, p), y = 0;
            for(register int j = 0; j < cnt; j++) {
                y = qmul(x, x, p);
                if(y == 1 && x != 1 && x != (p - 1)) 
                    return false;
                x = y;
            }
            if(y != 1) return false;
        }
        return true;
    }

    inline ll gcd(ll x, ll y) {
        register ll temp;
        while(y) {
            temp = x % y;
            x = y;
            y = temp;
        }
        return x;
    }

    inline void Pollard_Rho(ll x) {
        if(Miller_Rabin(x, 50)) { 
            ans = max(x, ans);
            return;
        }
        register ll s1 = rand() % (x - 1) + 1;
        register ll m = rand() % (x - 1) + 1, p = 1;
        register ll s2 = (qmul(s1, s1, x) + m) % x;
        for(register ll i = 1; s1 != s2; i++) { 
            p = qmul(p, abs(s1 - s2), x);
            if(p == 0) {
                register ll K = gcd(abs(s1 - s2), x);
                if(K != 1 && K != x) 
                    Pollard_Rho(K), Pollard_Rho(x / K);
                return; 
            }
            if(i % 127 == 0) {
                p = gcd(p, x);
                if(p != 1 && p != x) {
                    Pollard_Rho(p), Pollard_Rho(x / p);
                    return;
                }
                p = 1;
            }
            s1 = (qmul(s1, s1, x) + m) % x;
            s2 = (qmul(s2, s2, x) + m) % x;
            s2 = (qmul(s2, s2, x) + m) % x; 
        }
        p = gcd(p, x);
        if(p != 1 && p != x) {
            Pollard_Rho(p), Pollard_Rho(x / p);
            return;
        }
    }

    namespace IN_OUT {
        template <class T>
        inline void read(T &x) {
            static char c;
            static bool op;
            while(!isdigit(c = getchar()) && c != '-');
            x = (op = c == '-')? 0: c-'0';
            while(isdigit(c = getchar()))
                x = x * 10 + c - '0';
            if(op) x = ~x + 1; 
        }

        template <class T>
        inline void output(T x) {
            register char buf[30], *tail = buf;
            if(x == 0) putchar('0');
            else {
                if(x < 0) putchar('-'), x = ~x + 1;
                for(; x; x /= 10) *++tail = x % 10 + '0';
                for(; tail != buf; --tail)
                    putchar(*tail);
            }
            putchar('\n');
        }
    }

    using namespace IN_OUT;

    int main() {
        read(n);
        for(register int i = 1; i <= n; i++) {
            register ll x; 
            read(x);
            if(Miller_Rabin(x, 50))
                printf("Prime\n");

            else {
                ans=0;
                while(ans == 0)
                    Pollard_Rho(x);
                output(ans);
            }
        }
        return 0;
    }

下面是模题的评测记录，加了快读快输和一堆register \ ， \ repeat的值为50。

评测记录

※章末练习

  1. P2441 角色属性树

  2. P1069 细胞分裂

  3. P4358 [CQOI2016]密钥破解

  4. P4718 【模板】Pollard-Rho算法

END PS:

Pollard Rho 算法简介

【快速因数分解】Pollard's Rho 算法

PollardRho-快速分解质因数

质因数分解的题目较少，欢迎补充