浅谈质因数分解

浅谈质因数分解

->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\)是质数,无需分解。

  • Pollard-Rho 算法

试除法的时间复杂度在\(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-Pho\)算法 参考于:

Pollard Rho 算法简介

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

PollardRho-快速分解质因数

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

上一篇:OpenFOAM 可压缩湍流库深度解析


下一篇:密码学RSA解密之Pollard_rho分解