[科技]$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\)函数的一个原因，因为系统的随机并不一定会出现循环。而这个循环也构成了一个像\(\rho\)一样的形状，这也是其名称后缀\(Rho\)的来源。

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

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

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

然后发现最后一个点我们死活跑不过去md真毒瘤。然后我们继续考虑优化算法。发现整个算法的复杂度的瓶颈在于\(Pollard\_Rho\)中\(gcd\)上，虽然我们用了倍增算法减少了\(gcd\)的次数，但是这个复杂度仍然是无法承受的。但是需要注意到一件事情，我们只需要知道\(x\)的一个因子即可，并不需要做出规定一定是某个。然后一个显然的结论就是，根据欧几里得算法，\(gcd(a, b) > 1 \to 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\n", ans);
  }
  return 0;
}

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

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

原文地址：https://www.cnblogs.com/Apocrypha/p/10672354.html

时间： 2024-08-30 14:11:23

[科技]$Miller\_Rabin$ 和 $Pollard\_Rho$ 及各种玄学优化的相关文章

51_1037最长循环节 （miller rabin算法 pollard rho算法 原根）

1037 最长的循环节 V2 基准时间限制:1 秒 空间限制:131072 KB 分值: 320 难度:7级算法题 收藏 关注 正整数k的倒数1/k,写为10进制的小数如果为无限循环小数,则存在一个循环节,求<=n的数中,倒数循环节长度最长的那个数. 1/6= 0.1(6) 循环节长度为1 1/7= 0.(142857) 循环节长度为6 1/9= 0.(1)  循环节长度为1 Input 输入n(10 <= n <= 10^18) Output 输出<=n的数中倒数循环节长度最长的

bzoj3667: Rabin-Miller算法

Description Input 第一行:CAS,代表数据组数(不大于350),以下CAS行,每行一个数字,保证在64位长整形范围内,并且没有负数.你需要对于每个数字:第一,检验是否是质数,是质数就输出Prime 第二,如果不是质数,输出它最大的质因子是哪个. Output 第一行CAS(CAS<=350,代表测试数据的组数) 以下CAS行:每行一个数字,保证是在64位长整形范围内的正数. 对于每组测试数据:输出Prime,代表它是质数,或者输出它最大的质因子,代表它是和数 Sample In

poj1881：素因子分解+素数测试

很好的入门题 先测试是否为素数,若不是则进行素因子分解,算法详见总结贴 miller robin 和pollard rho算法 AC代码 #include <iostream> #include<stdio.h> #include<algorithm> #include<math.h> using namespace std; long long ans; long long gcd(long long a,long long b) { return b?g

ACM数论模板

w2w原创 1.欧几里得定理 int gcd(int a, int b){     return b==0?a:gcd(b, a%b); } 2.扩展欧几里得(求ax+by = gcd(a,b)的特解) void e_gcd(LL a, LL b, LL &d, LL &x, LL &y){ if(b==0){ x = 1; y = 0; d = a; } else{ e_gcd(b, a%b, d, y, x); y-= x*(a/b); } } 3.中国剩余定理 同余方程组 x

信息学竞赛知识点一览

C++语言 基础算法 位运算 快速幂 模拟 枚举 递推 递归 分治 二分 三分 排序 归并排序 离散化 倍增 贪心 高精度 数据结构 前缀和 差分 栈 对顶栈 单调栈 队列 双端队列 循环队列 单调队列 ST表 链表 链式前向星 Hash表 二叉堆 Huffman树 并查集 路径压缩 按秩合并 扩展域 边带权 树状数组 线段树 延迟标记 扫描线 动态开点线段树 分块 莫队 点分治 BST 平衡树 Treap Splay 红黑树 AVL SBT 替罪羊树 LCT CDQ分治 三维偏序 整体二分 可

HDU 3864 D_num Miller Rabin 质数判断+Pollard Rho大整数分解

链接:http://acm.hdu.edu.cn/showproblem.php?pid=3864 题意:给出一个数N(1<=N<10^18),如果N只有四个约数,就输出除1外的三个约数. 思路:大数的质因数分解只能用随机算法Miller Rabin和Pollard_rho,在测试多的情况下正确率是由保证的. 代码: #include <iostream> #include <cstdio> #include <cstring> #include <c

Pollard rho算法+Miller Rabin算法 BZOJ 3668 Rabin-Miller算法

BZOJ 3667: Rabin-Miller算法 Time Limit: 60 Sec  Memory Limit: 512 MBSubmit: 1044  Solved: 322[Submit][Status][Discuss] Description Input 第一行:CAS,代表数据组数(不大于350),以下CAS行,每行一个数字,保证在64位长整形范围内,并且没有负数.你需要对于每个数字:第一,检验是否是质数,是质数就输出Prime 第二,如果不是质数,输出它最大的质因子是哪个. O

Miller&amp;&amp;Pollard HDOJ 4344 Mark the Rope

题目传送门 题意:一个长为n(n<2^63)的管子,在管子上做标记,每隔L个长度单位做一个标记,从管子头端开始,保证最后一次标记恰好在管子的尾端.让你找出有多少个这样的L(L<n),且他们之间两两互素,然后求出这些L的和最大值. 分析:转换一下就是求n有多少个质因子用pollard_rho大整数分解分解n,因为素数之间两两互质,所以每段L都由每个质因子的k次幂组成,如果n是素数,由于L<n,所以只能L==1 收获:接触到随机算法 代码: /************************

Miller&amp;&amp;Pollard POJ 1811 Prime Test

题目传送门 题意:素性测试和大整数分解, N (2 <= N < 254). 分析:没啥好讲的,套个模板,POJ上C++提交 收获:写完这题得到模板 代码: /************************************************ * Author :Running_Time * Created Time :2015-8-28 13:02:38 * File Name :POJ_1811.cpp ************************************