• [科技]$Miller\_Rabin$ 和 $Pollard\_Rho$ 及各种玄学优化


    [科技](Miller\_Rabin)(Pollard\_Rho) 及各种玄学优化

    [科技] (Miller\_Rabin)(Pollard\_Rho)

    先讲(Miller\_Rabin)吧,(Miller\_Rabin)是用来检验素数的高效算法。

    我们先要知道两个定理

    1. 费马小定理:当(p)为质数时,(x^{p - 1} equiv 1 (mod p))。但这只是一个充分条件,但不是必要条件。即就算(x)(p)互质,那么(x^{p - 1})不一定在模(p)意义下同余于(1)
    2. (Fermat)定理:若(p)为奇素数,且(0 < x < p),那么(x ^ 2 equiv 1 (mod p))的解为(x = 1)(x = p - 1)。这个还是比较好证明的,移项可得(x ^ 2 - 1 equiv 0),即((x - 1)(x + 1) equiv 0),可得(p | (x - 1)(x + 1))。而(p)是质数,那么可知(x = 1)(x = p - 1)

    首先我们根据费马小定理就可以排除大量的合数了,即如果存在(x),不满足(x ^ {p - 1} equiv 1),则(p)就不是质数了。但是之前也说过,这个并不是一个必要条件,所以也存在满足费马小定理的合数,这种数叫做(Carmichael)数。最小的(Carmichael)数是(561),对于任意的数,都有(x ^ {560} equiv 1),而(561 = 3 * 11 * 17)(Carmichael)数使得我们无法使用费马小定理来判断一个数是否为素数。

    然后(Miller\_Rabin)算法就是优化了这一点,尝试使用多次判定来提高正确性。当我们判断了一个数满足了费马小定理之后,我们继续判断(x ^ { (p - 1) / 2} equiv 1)是否正确。于是我们就得出了一个算法:我们将(p - 1)表示成(2 ^ s * t)的形式,然后选取几个比较小的质数(pri),从(pri ^ t)开始,依次判断(pri ^ t, pri ^ {2t},pri ^ {4t} dots)等等是否满足这个同余式。如果对于所有取的素数都满足时,就可以确定这个数大概率是个素数了。经过测试发现,对于(int)范围以内的数,取出小于(30)的素数来验证就可以保证完全正确的判断出是否为素数,在(long long)范围内的错误的概率也是可以忽略不计的。

    再来讲讲(Phollard\_Rho),他是基于(Miller\_Rabin)的一种分解质因数的方法,同时也是基于随机化的方法,期望复杂度为(O(n ^ {frac{1}{4}})),但实际情况下的复杂度也是个玄学,似乎表现的非常优秀。

    假设我们需要对于(x)分解质因数,我们先通过(Miller\_Rabin)素性测试判断给定的(x)是否为素数。然后开始(Pollard\_Rho)进行分解。

    我们先要了解的是(Pollard\_Rho)是将某合数(x)分解为两个非平凡因子(a, b)(1)(x)为平凡因子)的算法。如果需要对(x)分解质因数,那么递归下去做就行了。

    先思考一下平时的(O(sqrt n))分解质因数的方法,即枚举(2)(sqrt n)的质因数,判断是否能被整除。实际上一个数的质因数是(log)级别的,所以我们其中的枚举很多都是无用枚举。考虑如何来优化这个枚举方法。

    有一个悖论叫做生日悖论,意思是在一个人数为(23)人的班级中,存在两个人生日相同的概率接近(0.5),而在人数大约为(60)人的班级中,存在两个人生日相同的概率已经非常接近(1)了。这虽然与实际非常不相符,但是简单计算一下发现的确是这样。以(23)人为例,我们计算没有两个人生日相同的概率(P = frac{365}{365} * frac{364}{365} * frac{363}{365} * dots * frac{365 - 22}{365}),此时的(P)大约为(0.49)

    举一个简单的例子就是假设我们需要在([1, 100])之间选一个数,那么选到(39)的概率为(1 \%),但是如果我们选择两个数(a, b),使得(| a - b | = 39),那么概率就变成了(2 \%)

    这给了我们启发,我们在选择测试因子时,也采用这种组合随机采样方法。我们随机两个数(a, b),判断(gcd( | a - b |, n))是否大于1,如果大于(1),则我们找到了一个因子,这样会大大提高我们找到因子的概率。

    那么我们该如何生成这些随机数呢?(Pollard\_Rho)通过伪随机数的方法来提供待测因子,即选定一个起始数字(x_0),以及一个常数(c),通过(x_i = x_{i - 1} ^ 2 + c (mod n))的方法来生成一系列的随机因子。至于为什么选用这个函数呢?似乎是这个函数在带入复数之后迭代得出的点集叫曼德勃罗集,这个集合又与混沌系统有关,即这个集合中的一系列数字非周期又不收敛,这使得这个函数生成的随机数非常优秀。并且这些生成的随机数,每一个都完全依赖于前一个数字,那么在迭代了一定次数之后,一定会出现循环。所以当我们出现循环的时候,就重新设置(x_0)(c),继续检测,这是不用系统的(rand)函数的一个原因,因为系统的随机并不一定会出现循环。而这个循环也构成了一个像( ho)一样的形状,这也是其名称后缀(Rho)的来源。

    然后我们就可以设计一个基于倍增的随机算法了,每次在( ho)的路径上取([2 ^ {k - 1}, 2 ^ k])这段区间,然后对于所有的(i in [2 ^ {k - 1}, 2 ^ k]),判断(gcd(x_i - x_{2^{k - 1}}, n))是否大于1,如果大于1,则找到了一个因子,返回即可。这样既可以不在( ho)的环上停留太久,也可以减少(gcd)的次数。

    于是我们得出了(Pollard\_Rho)的大致算法:

    1. 对于需要检验的数(n),用(Miller\_Rabin)进行素性测试,如果是素数,则直接返回即可,否则进入第二步。
    2. 随机生成数(x_0, c),然后开始倍增的判断(x_i - x_l)(n)是否有公因数,如果找到了公因数,返回即可。若已经出现了循环,则返回(n),表示该次探查失败,并重复2步骤。

    这样我们就可以非常快速的写出(Pollard\_Rho)的算法了,这里放出例题

    Code:

    #pragma GCC optimize (2, "inline", "Ofast")
    #include <bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    int T;
    ll n, ans = 0;
    
    ll QMul(ll x, ll y, ll Md) {
      ll ans = 0;
      for(; y; y >>= 1, x = (x + x) % Md) if(y & 1) ans = (ans + x) % Md;
      return ans;
    }
    
    ll Qpow(ll x, ll y, ll Md) {
      ll ans = 1;
      for(; y; y >>= 1, x = QMul(x, x, Md)) if(y & 1) ans = QMul(ans, x, Md);
      return ans;
    }
    
    bool Miller_Rabin(ll x) {
      if(x == 2) return 1;
      if(!(x & 1) || x == 1) return 0;
      ll s = 0, t = x - 1;
      while(!(t & 1)) {
        t >>= 1;
        s ++;
      }
      for(int i = 0; i < 10 && pri[i] < x; i++) {
        ll a = pri[i], b = Qpow(a, t, x);
        ll k;
        for(int j = 1; j <= s; j++) {
          k = QMul(b, b, x);
          if(k == 1 && b != 1 && b != x - 1) return 0;
          b = k;
        }
        if(b != 1) return 0;
      }
      return 1;
    }
    
    ll Gcd(ll a, ll b) {
      return !b ? a : Gcd(b, a % b);
    }
    
    ll Pollard_Rho(ll n, ll c) {
      ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
      int i = 1, k = 2;
      while(1) {
        i++;
        x = QMul(x, x, n);
        x = (x + c) % n;
        ll G = Gcd(y - x, n);
        if(G > 1 && G < n) return G;
        if(y == x) return n;
        if(i == k) {
          y = x;
          k <<= 1;
        }
      }
    }
    
    
    void Find(ll x, ll c) {
      if(x == 1 || x < ans) return ;
      if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
      ll p = x;
      while(p == x) p = Pollard_Rho(p, c--);
      while(x % p == 0) x /= p;
      Find(p, c); Find(x, c);
    }
    
    int main() {
      srand(2333);
      scanf("%d", &T);
      while(T--) {
        scanf("%lld", &n);
        ans = 0;
        bool f = Miller_Rabin(n);
        if(f) {
          puts("Prime");
          continue;
        }
        Find(n, 1000000);
        printf("%lld
    ", ans);
      }
      return 0;
    }
    

    好的……本科技讲解到此结束
    (nmdwsm!!!)
    为什么交到(luogu)的模板题上(T)飞了!!本机(0.7s) (luogu)上居然(T)了,难以置信……

    于是我们开始漫漫的卡常优化之旅……

    首先注意到我们在乘法的时候为了防止(long long)相乘导致溢出的情况,我们使用了龟速乘……实际上它确实是龟速,我们在乘法时强转成(\_\_int128),就可以避免溢出的情况了。

    Code:

    #pragma GCC optimize (2, "inline", "Ofast")
    #include <bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    int T;
    ll n, ans = 0;
    
    ll Qpow(ll x, ll y, ll Md) {
      ll ans = 1;
      for(; y; y >>= 1, x = x * x % Md) if(y & 1) ans = ans * x % Md;
      return ans;
    }
    
    bool Miller_Rabin(ll x) {
      if(x == 2) return 1;
      if(!(x & 1) || x == 1) return 0;
      ll s = 0, t = x - 1;
      while(!(t & 1)) {
        t >>= 1;
        s ++;
      }
      for(int i = 0; i < 10 && pri[i] < x; i++) {
        ll a = pri[i], b = Qpow(a, t, x);
        ll k;
        for(int j = 1; j <= s; j++) {
          k = (__int128)b * b % x;
          if(k == 1 && b != 1 && b != x - 1) return 0;
          b = k;
        }
        if(b != 1) return 0;
      }
      return 1;
    }
    
    ll Gcd(ll a, ll b) {
      return !b ? a : Gcd(b, a % b);
    }
    
    ll Pollard_Rho(ll n, ll c) {
      ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
      int i = 1, k = 2;
      while(1) {
        i++;
        x = (__int128)x * x % n;
        x = (x + c) % n;
        ll G = Gcd(y - x, n);
        if(G > 1 && G < n) return G;
        if(y == x) return n;
        if(i == k) {
          y = x;
          k <<= 1;
        }
      }
    }
    
    
    void Find(ll x, ll c) {
      if(x == 1 || x < ans) return ;
      if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
      ll p = x;
      while(p == x) p = Pollard_Rho(p, c--);
      while(x % p == 0) x /= p;
      Find(p, c); Find(x, c);
    }
    
    int main() {
      srand(2333);
      scanf("%d", &T);
      while(T--) {
        scanf("%lld", &n);
        ans = 0;
        bool f = Miller_Rabin(n);
        if(f) {
          puts("Prime");
          continue;
        }
        Find(n, 1000000);
        printf("%lld
    ", ans);
      }
      return 0;
    }
    

    于是我们修改之后再次交一发:

    (nmdwsm!!!)为什么这样就快了那么多啊……真是

    然后发现最后一个点我们死活跑不过去md真毒瘤。然后我们继续考虑优化算法。发现整个算法的复杂度的瓶颈在于(Pollard\_Rho)(gcd)上,虽然我们用了倍增算法减少了(gcd)的次数,但是这个复杂度仍然是无法承受的。但是需要注意到一件事情,我们只需要知道(x)的一个因子即可,并不需要做出规定一定是某个。然后一个显然的结论就是,根据欧几里得算法,(gcd(a, b) > 1 o gcd(a * c \% b, b) > 1)。所以我们可以把一次倍增中所有的随机因子都乘起来之后,一起做一次(gcd),这样原本是互质的,现在仍然互质,原本存在公因数的,现在仍然存在。这样就可以大大减少(gcd)的次数,唯一需要处理的边界条件就是所有随机因子的乘积与(n)(gcd)之后,(gcd)(n),这个时候虽然存在公因数,但是却不能直接返回(n)(因为(n)是我们用来判断本次探测不成功的标记),这个时候我们再按照原方法枚举一遍,找到某个非平凡因子返回即可。

    同时(gcd)的地方也是可以优化很多的,大致就是仿照高精度(gcd)那样用二进制处理,常数会小很多。

    Code:

    #pragma GCC optimize(3, "inline", "Ofast")
    #include <bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int pri[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
    int T;
    ll n, ans = 0;
    
    ll Qpow(ll x, ll y, ll Md) {
      ll ans = 1;
      for(; y; y >>= 1, x = x * x % Md) if(y & 1) ans = ans * x % Md;
      return ans;
    }
    
    bool Miller_Rabin(ll x) {
      if(x == 2) return 1;
      if(!(x & 1) || x == 1) return 0;
      ll s = 0, t = x - 1;
      while(!(t & 1)) {
        t >>= 1;
        s ++;
      }
      for(int i = 0; i < 10 && pri[i] < x; i++) {
        ll a = pri[i], b = Qpow(a, t, x);
        ll k;
        for(int j = 1; j <= s; j++) {
          k = (__int128)b * b % x;
          if(k == 1 && b != 1 && b != x - 1) return 0;
          b = k;
        }
        if(b != 1) return 0;
      }
      return 1;
    }
    
    ll Gcd(ll a, ll b) {
      if(!a || !b) return a | b;
      int t = __builtin_ctzll(a | b); //__builtin_ctzll是得到某个数二进制下末尾0的个数
      a >>= __builtin_ctzll(a);
      do {
        b >>= __builtin_ctzll(b);
        if(b < a) swap(a, b);
        b -= a;
      }while(b);
      return a << t;
    }
    
    ll Pollard_Rho(ll n, ll c) {
      ll x = 1ll * rand() * rand() % (n - 2) + 1, y = x;
      int i = 1, k = 2;
      ll Mu = 1, z = y;
      while(1) {
        i++;
        x = (__int128) x * x % n;
        x = (x + c) % n;
        Mu = (__int128) Mu * (y > x ? y - x : x - y) % n;
        if(y == x) return n;
        if(i == k) {
          ll G = Gcd(Mu, n);
          if(G == 1) { y = x; k <<= 1; continue; }
          if(G == n) {
            x = y;
            for(int t = 1; t <= k >> 1; t++) {
              x = (__int128) x * x % n;
              x = (x + c) % n;
              ll g = Gcd(y > x ? y - x : x - y, n);
              if(g > 1 && g < n) return g;
            }
          }
          return G;
        }
      }
    }
    
    
    void Find(ll x, ll c) {
      if(x == 1 || x < ans) return ;
      if(Miller_Rabin(x)) return (void) (ans = max(ans, x));
      ll p = x;
      while(p == x) p = Pollard_Rho(p, c--);
      while(x % p == 0) x /= p;
      Find(p, c); Find(x, c);
    }
    
    int main() {
      srand(2333);
      scanf("%d", &T);
      while(T--) {
        scanf("%lld", &n);
        ans = 0;
        bool f = Miller_Rabin(n);
        if(f) {
          puts("Prime");
          continue;
        }
        Find(n, 100000);
        printf("%lld
    ", ans);
      }
      return 0;
    }
    
    

    尝试再次交一发,愉快的发现它过了。爆OJ真舒服

    然后还有一个关于(127)这个数的优化,现在还没有看懂,咕咕咕~

  • 相关阅读:
    PAIP.MYSQL数据库比较
    paip.验证码识别----判断汉字还是英文
    SQLServer2008客户端软件
    paip.多个TOMCAT共存在一台主机上配置方法
    paip.银行卡号的发卡行归属地查询
    paip.获取当前实际北京时间API
    PAIP.HIBERNATE ORA02289 sequence does not exist的解决
    C51与汇编语言混合编程之一
    KEIL C51高级编程之二
    可重入函数与不可重入函数(转)
  • 原文地址:https://www.cnblogs.com/Apocrypha/p/10672354.html
Copyright © 2020-2023  润新知