素數判定(Miller_Rabin) & 大數分解(Pollard_Pho)

    · 博主語文體育老師教的.

    · 本文年齡限定 16+

    · 吐槽上面 2條的都是⑨

    事情要從昨天考試說起.

    Pro 2 : 給定n, 求n最少能分成多少個完全平方數.

    那啥做法就不扯淡了, 各種特判. 

    核心演算法的素數判定以及大數分解完全就被一大堆特判淹沒了.

    推薦:

        1、《64位以內Rabin-Miller+強偽素數測試和Pollard+rho+因數分解演算法的實現》 (裡面的虛擬碼非常奇葩.... 傳說中的PC語言 ?!)

        2、CSDN Fisher_jiang《Miller_Rabin素數測試》- http://blog.csdn.net/fisher_jiang/article/details/986654

        3、《演算法導論》自己翻去.

    素數判定 (Miller_Rabin) 

    

    個人就不吐槽了.

    總之就是用費馬小定理 a^(p - 1) ≡ 1 (mod p)  (其中 a, p 互質且 p 為質數)  確定 p 是否為質數的較高正確性.

    證明不管 (確切地說沒有精力去看), 只剽演算法=. =

    只要多選取幾個 a , 就可以使得演算法正確率非常高. 目前有以下二種方法:

        1、隨機取.

        2、選取前 k 個質數.      

    第一種出錯率 1 / 4^s ( s為隨機選取 a 的個數),  第二種在 [1] 中有詳細介紹, 取前 10 個質數就可以滿足在int64範圍內100%正解.

    那啥怎麼求a^(p - 1) ?  快速冪...... ( 不可能有人都看到這種東西了還不會快速冪吧......)

    大數分解 (Pollard_Pho)

    極度吐槽. 讓我悲劇了一下午的東西.

    總之就是隨機選 2 個數a, b (b < a < n),  判斷 gcd(a - b, n) 是否大於 1 就可以了.

    也許會說這種隨機怎麼可能那麼快出解啊! 什麼的.

    但是若構造一個迴圈節 (a1, a2, a3, a4, a5 ... ak) (在模n意義下迴圈),  並且使 b 也有一定規律的話.

    根據各種論文,  能在O(√p) 內找出 n 的一個因子 p (注意是因子, 不是質因子, 所以找出 p 後要遞迴處理)

    一般用 f(x) = x² + c 來構造迴圈節. 其餘證明見 [1]  (個人認為應該是算導以外最好的了)

    

    以下是虛擬碼 :

    Function Pho(n)

        if N IS A PRIME  then                                 // n 是個質數就直接退出

          return error;

        x = y = x0;   c = random[1, n);                   //  初始化

        k = 0;   i = 1;   d = 1;

        while (true)   do

            ++k;                                                          

            x = f(x), d = gcd(x - y, n);                       //  構造 x

            if  d ∈ (1, n)  then  return d;                //   找到一個因子

            if  d = n  then c = random[1, n);          //   x - y = 0, 表示 c 這個常數因子有點萎... 換一個.  （悲劇一下午, Orz 屏屏哥)

            if  k = i  then y = x, i <<= 1;                  //   更新 y

    End Function

    以上是Pollard 的 Brent最佳化代碼, 同樣, 具體見 [1].

    下面是完全平方數分解的 Code ( 裡面內建素數判定以及大數的質因數分解 )

#include <cmath>#include <cstdio>#include <cstdlib>#include <iostream>#include <algorithm>#define error -1#define ll long longconst int po[11] = {0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29};using namespace std;int test;ll n, p[10010];ll gcd(ll a, ll b){    ll x;    while (x = a, b) a = b, b = x % b;    return a; }ll mul(ll a, ll b, ll c){    ll plas = 0;    for (ll i = 1; i <= b; i <<= 1, a = (a + a) % c)      if (b & i)  plas = (plas + a) % c;    return plas;  }ll Rollard_Brent(ll n){    ll x = 1, y = 1, d = 1;  y = x;    int k = 0, i = 1, z = rand() + 1;    while (true)      {  ++k;         x = (mul(x, x, n) + z) % n;         d = gcd((x - y + n) % n, n);         if (d > 1  &&  d < n)  return d;         if (k == i)  y = x, i <<= 1;         if (d == n) z = rand() + 1;      }}    bool Miller_Rabin(ll n){    ll a, w, sec;    for (int i = 1; i <= 10; ++i)      {   sec = 1;  if (n == po[i])  continue;          for (a = po[i], w = 1; w < n; w <<= 1, a = a * a % n)            if ((n - 1) & w)  sec = sec * a % n;           if (sec != 1)  return false;      }        return true;}int tap(ll n){    ll plas = n;    int head = 1, tail = 1, top = 0;    p[tail] = n;     while (head <= tail  &&  n != 1)      if (Miller_Rabin(p[head]))  {  bool wis = 0;         while (n % p[head] == 0)  wis = !wis, n /= p[head];        if (wis  &&  (p[head] & 3) == 3)  goto Compare;  ++head;      }            else  {        ll vec = Rollard_Brent(p[head]);        p[++tail] = p[head] / vec;  p[++tail] = vec;          ++head;  }      return 2;    Compare:    while (!(plas & 3))  plas >>= 2;    if (((plas - 7) & 7) == 0)  return 4;    else  return 3;}int main(){    freopen("p2.in", "r", stdin);    freopen("p2.out", "w", stdout);    srand((unsigned)time(NULL));    scanf("%d", &test);    for (; test--; )      {                scanf("%I64d", &n);            int w = (int) sqrt((double) n);  if ((ll) w * w == n)  {  printf("1\n");  continue;  }          printf("%d\n", tap(n));      }    }