Pollard-rho的质因数分解

  1 #include <iostream>
  2 #include <time.h>
  3 #include <map>
  4 using namespace std;
  5 long long an[] = {2,3,5,7,11,13,17,61};
  6 map<long long,int> mp;//存因数和对应出现次数
  7 int num = 0;
  8 long long Random(long long n)//生成0到n之间的整数
  9 {
 10     return (double) rand()/RAND_MAX*n+0.5;//(doubel)rand()/RAND_MAX生成0-1之间的浮点数
 11 }
 12 
 13 
 14 long long q_mod(long long a,long long n,long long p)//快速幂
 15 {
 16     a = a%p;
 17     //首先降a的规模
 18     long long sum = 1;//记录结果
 19     while(n)
 20     {
 21         if(n&1)
 22         {
 23             sum = (sum*a)%p;//n为奇数时单独拿出来乘
 24         }
 25         a = (a*a)%p;//合并a降n的规模
 26         n /= 2;
 27     }
 28     return sum;
 29 }
 30 
 31 
 32 long long q_mul(long long a,long long b,long long p)//大数模
 33 {
 34     long long sum = 0;
 35     while(b)
 36     {
 37         if(b&1)//如果b的二进制末尾是零
 38         {
 39             (sum += a)%=p;//a要加上取余
 40         }
 41         (a <<= 1)%=p;//不断把a乘2相当于提高位数
 42         b >>= 1;//把b右移
 43     }
 44     return sum;
 45 }
 46 
 47 
 48 //Miller-Rabin
 49 bool witness(long long a,long long n)
 50 {
 51     long long d = n-1;
 52     long long r = 0;
 53     while(d%2==0)
 54     {
 55         d/=2;
 56         r++;
 57     }//n-1分解成d*2^r，d为奇数
 58     long long x = q_mod(a,d,n);
 59     //cout << "d " << d << " r " << r << " x " << x << endl;
 60     if(x==1||x==n-1)//最终的余数是1或n-1则可能是素数
 61     {
 62         return true;
 63     }
 64     while(r--)
 65     {
 66         x = q_mul(x,x,n);
 67         if(x==n-1)//考虑开始在不断地往下余的过程
 68         {
 69             return true;//中间如果有一个余数是n-1说明中断了此过程，则可能是素数
 70         }
 71     }
 72     return false;//否则如果中间没有中断但最后是余数又不是n-1和1说明一定不是素数
 73 }
 74 bool miller_rabin(long long n)
 75 {
 76     const int times = 50;//试验次数
 77     if(n==2)
 78     {
 79         return true;
 80     }
 81     if(n<2||n%2==0)
 82     {
 83         return false;
 84     }
 85     for(int i = 0;i<times;i++)
 86     {
 87         long long a = Random(n-2)+1;//1到(n-1)
 88         //cout << a << endl;
 89         if(!witness(a,n))
 90         {
 91             return false;
 92         }
 93     }
 94     return true;
 95 }
 96 
 97 
 98 //求gcd
 99 long long gcd(long long a,long long b)
100 {
101     return b==0?a:gcd(b,a%b);
102 }
103 
104 //Pollard-rho
105 long long Pollard_rho(long long n,long long c)
106 {
107     //cout << "n "  << n << " c " << c << endl;
108     long long i = 1,k = 2;
109     long long x = Random(n-1)+1;
110     long long y = x;
111     while(true)
112     {
113         i++;
114         x = (q_mul(x,x,n)+c)%n;
115         long long d = gcd(y-x,n);
116         if(1<d&&d<n)
117         {
118             return d;
119         }
120         if(x==y)
121         {
122             return n;
123         }
124         if(i==k)
125         {
126             y = x;
127             k<<=1;
128         }
129     }
130 }
131 void find(long long n,long long c)
132 {
133     if(n==1)//找完了
134     {
135         return;
136     }
137     if(miller_rabin(n))//找到了质数
138     {
139         num++;
140         mp[n]++;
141         return;
142     }
143     long long p = n;
144     while(p>=n)//找p的因数
145     {
146         p = Pollard_rho(p,c--);//返回p的因数或1或本身
147     }
148     find(p,c);//递归地找p的因子
149     find(n/p,c);
150 }
151 int main()
152 {
153     long long n;
154     while(cin >> n)
155     {
156         num = 0;
157         mp.clear();
158         find(n,2137342);//随机选取的c
159         cout << n << " = ";
160         if(mp.empty())
161         {
162             cout << n << endl;
163         }
164         for(auto ite = mp.begin();ite!=mp.end();ite++)
165         {
166             cout << ite->first << "^" << ite->second;
167             auto i = ite;
168             if(++i!=mp.end()) //如果不是最后一个
169             {
170                 cout << "*";//输出乘号
171             }
172         }
173     }
174     return 0;
175 
176 }

 1 map<int,int> mp;
 2 void decom1(int n)
 3 {
 4     for(int i = 2;i<=n;i++)
 5     {
 6         while(n%i==0)
 7         {
 8             mp[i]++;
 9             n /= i;
10         }
11     }
12 }