大數因式分解 Pollard_rho 算法詳解


給你一個大數n,將它分解它的質因子的乘積的形式。

首先需要了解Miller_rabin判斷一個數是否是素數

大數分解最簡單的思想也是試除法,這里就不再展示代碼了,就是從2到sqrt(n),一個一個的試驗,直到除到1或者循環完,最后判斷一下是否已經除到1了即可。

 

但是這樣的做的復雜度是相當高的。一種很妙的思路是找到一個因子(不一定是質因子),然后再一路分解下去。這就是基於Miller_rabin的大數分解法Pollard_rho大數分解。

 

Pollard_rho算法的大致流程是 先判斷當前數是否是素數(Miller_rabin)了,如果是則直接返回。如果不是素數的話,試圖找到當前數的一個因子(可以不是質因子)。然后遞歸對該因子和約去這個因子的另一個因子進行分解。

 

那么自然的疑問就是,怎么找到當前數n的一個因子?當然不是一個一個慢慢試驗,而是一種神奇的想法。其實這個找因子的過程我理解的不是非常透徹,感覺還是有一點兒試的意味,但不是盲目的枚舉,而是一種隨機化算法。我們假設要找的因子為p,他是隨機取一個x1,由x1構造x2,使得{p可以整除x1-x2 && x1-x2不能整除n}則p=gcd(x1-x2,n),結果可能是1也可能不是1。如果不是1就找尋成功了一個因子,返回因子;如果是1就尋找失敗,那么我們就要不斷調整x2,具體的辦法通常是x2=x2*x2+c(c是自己定的)直到出現x2出現了循環==x1了表示x1選取失敗重新選取x1重復上述過程。(似乎還存在一個每次找尋范圍*2的優化,但是不太懂。。。)

 

因為x1和x2再調整時最終一定會出現循環,形成一個類似希臘字母rho的形狀,故因此得名。

 

 

另外通過find函數來分解素數,如果找到了一個素數因子則加入到因子map中,否則如果用Pollard找到一個因子則遞歸去找素數因子。

  1 #include<iostream>
  2 #include<ctime>
  3 #include<algorithm>
  4 #include<map>
  5 using namespace std;
  6 typedef long long ll;
  7 map<ll, int>m;
  8 const int mod = 10000019;
  9 const int times = 50;//測試50次
 10 ll mul(ll a, ll b, ll m)
 11 //求a*b%m
 12 {
 13     ll ans = 0;
 14     a %= m;
 15     while(b)
 16     {
 17         if(b & 1)ans = (ans + a) % m;
 18         b /= 2;
 19         a = (a + a) % m;
 20     }
 21     return ans;
 22 }
 23 ll pow(ll a, ll b, ll m)
 24 //a^b % m
 25 {
 26     ll ans = 1;
 27     a %= m;
 28     while(b)
 29     {
 30         if(b & 1)ans = mul(a, ans, m);
 31         b /= 2;
 32         a = mul(a, a, m);
 33     }
 34     ans %= m;
 35     return ans;
 36 }
 37 bool Miller_Rabin(ll n, int repeat)//n是測試的大數,repeat是測試重復次數
 38 {
 39     if(n == 2 || n == 3)return true;//特判
 40     if(n % 2 == 0 || n == 1)return false;//偶數和1
 41 
 42     //將n-1分解成2^s*d
 43     ll d = n - 1;
 44     int s = 0;
 45     while(!(d & 1)) ++s, d >>= 1;
 46     //srand((unsigned)time(NULL));在最開始調用即可
 47     for(int i = 0; i < repeat; i++)//重復repeat次
 48     {
 49         ll a = rand() % (n - 3) + 2;//取一個隨機數,[2,n-1)
 50         ll x = pow(a, d, n);
 51         ll y = 0;
 52         for(int j = 0; j < s; j++)
 53         {
 54             y = mul(x, x, n);
 55             if(y == 1 && x != 1 && x != (n - 1))return false;
 56             x = y;
 57         }
 58         if(y != 1)return false;//費馬小定理
 59     }
 60     return true;
 61 }
 62 ll gcd(ll a, ll b)
 63 {
 64     return b == 0 ? a : gcd(b, a % b);
 65 }
 66 ll pollard_rho(ll n, ll c)//找到n的一個因子
 67 {
 68     ll x = rand() % (n - 2) + 1;
 69     ll y = x, i = 1, k = 2;
 70     while(1)
 71     {
 72         i++;
 73         x = (mul(x, x, n) + c) + n;//不斷調整x2
 74         ll d = gcd(y - x, n);
 75         if(1 < d && d < n)
 76             return d;//找到因子
 77         if(y == x)
 78             return n;//找到循環,返回n,重新來
 79         if(i == k)//一個優化
 80         {
 81             y = x;
 82             k <<= 1;
 83         }
 84     }
 85 }
 86 void Find(ll n, ll c)
 87 {
 88     if(n == 1)return;//遞歸出口
 89 
 90     if(Miller_Rabin(n, times))//如果是素數,就加入
 91     {
 92         m[n]++;
 93         return;
 94     }
 95 
 96     ll p = n;
 97     while(p >= n)
 98         p = pollard_rho(p, c--);//不斷找因子,知道找到為止,返回n說明沒找到
 99 
100     Find(p, c);
101     Find(n / p, c);
102 }
103 int main()
104 {
105     ll n;srand((unsigned)time(NULL));
106     while(cin >> n)
107     {
108         m.clear();
110 Find(n, rand() % (n - 1) + 1);//這是自己設置的一個數 111 cout<<n<<" = "; 112 for(map<ll ,int>::iterator it = m.begin(); it != m.end();) 113 { 114 cout<<it->first<<" ^ "<<it->second; 115 if((++it) != m.end()) 116 cout<<" * "; 117 } 118 cout<<endl; 119 } 120 return 0; 121 }

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM