POJ 2429 GCD & LCM Inverse(Miller-Rabbin素性测试,Pollard rho质因子分解)
2015-10-20 22:09
369 查看
x = lcm/gcd,假设答案为a,b,那么a*b = x且gcd(a,b) = 1,因为均值不等式所以当a越接近sqrt(x),a+b越小。
x的范围是int64的,所以要用Pollard_rho算法去分解因子。因为a,b互质,所以我们把相同因子一起处理。
最多16个不同的因子:2,3,5,7,11,13,17,19,23,29,31,37,41,43,47, 乘积为 614889782588491410, 乘上下一个质数53会爆int64范围。
所以剩下暴力枚举一下就好。
x的范围是int64的,所以要用Pollard_rho算法去分解因子。因为a,b互质,所以我们把相同因子一起处理。
最多16个不同的因子:2,3,5,7,11,13,17,19,23,29,31,37,41,43,47, 乘积为 614889782588491410, 乘上下一个质数53会爆int64范围。
所以剩下暴力枚举一下就好。
#include<cstdio> #include<iostream> #include<string> #include<cstring> #include<queue> #include<vector> #include<stack> #include<vector> #include<map> #include<set> #include<algorithm> #include<cmath> #include<ctime> //#include<bits/stdc++.h> using namespace std; typedef long long ll; ll mulMod(ll a, ll b, ll m) // a %= m { ll re = 0; while(b){ if(b&1){ re = (re+a); if(re >= m) re -= m; } a <<= 1; if(a >= m) a -= m; b >>= 1; } return re; } ll powMod(ll a,ll q,ll m) { ll re = 1; while(q){ if(q&1){ re = mulMod(re,a,m); } a = mulMod(a,a,m); q >>= 1; } return re; } bool witness(ll a,ll d,int t,ll n) { ll x = powMod(a,d,n), y; while(t--){ y = mulMod(x,x,n); if(y == 1 && x != 1 && x != n-1) return true; x = y; } return x != 1; } bool RabinMiller(ll n,int k = 10) { if(n == 2) return true; if(n < 2 || ((n&1)^1) ) return false; int t = 0; ll d = n-1; while((d&1)^1) { d>>=1; t++; } while(k--){ if(witness(rand()%(n-1)+1,d,t,n)) return false; } return true; } ll gcd(ll a,ll b) { while(b){ ll t = a%b; a = b; b = t; } return a < 0? -a: a; } ll rho(ll n) { ll x = rand()%(n-1) + 1; ll y = x, d; int i = 1, k = 2; ll c = 1+rand()%(n-1); while(true){ i++; x = (mulMod(x,x,n) + c); if(x >= n) x -= n; d = gcd(y-x,n); if(d != 1 && d != n) return d; if(y == x) return n; if(i == k){ y = x; k <<= 1; } } } ll Prms[150], stk[150]; int tot; void pollard(ll n) { tot = 0; int top = 0; stk[++top] = n; while(top){ n = stk[top--]; if(RabinMiller(n)){ Prms[tot++] = n; }else{ ll p = n; while(p >= n) p = rho(p); stk[++top] = p; stk[++top] = n/p; } } } ll mpow(ll a,int q) { ll re = 1; while(q){ if(q&1) re *= a; a *= a; q>>=1; } return re; } //#define LOCAL int main() { #ifdef LOCAL freopen("in.txt","r",stdin); #endif ll g,l; while(~scanf("%I64d%I64d",&g,&l)){ if(l == g){ printf("%I64d %I64d\n",g,l); continue; } ll x = l/g; pollard(x); sort(Prms,Prms+tot); Prms[tot] = -1; int k = 0; for(int j = 0,i = 1, cnt = 1; i <= tot; i++) { if(Prms[i] == Prms[j]) cnt++; else { Prms[k++] = mpow(Prms[j],cnt); cnt = 1; j = i; } } ll sqr = floor(sqrt(x)), best = 1; for(int S = 1<<(tot = k); --S; ){ ll cur = 1; for(int i = 0; i < tot; i++){ if(S>>i&1) cur *= Prms[i]; if(cur > sqr) break; } if(cur <= sqr && cur > best) best = cur; } printf("%I64d %I64d\n",g*best,l/best); } return 0; }
相关文章推荐
- C++中常用的三种参数传递方式
- 类似花生壳的代理本机服务至公网的ngrok以及pagekite
- 学习笔记06-多线程
- 【转】windows 签发 pfx 证书
- 深度优先算法
- 如何生成带注释的DLL文件
- iOS线程的简单学习<1>
- iOS初探+load和+initialize
- Linux下LVM配置过程
- 单维与多维线性回归代码( machine-learning-ex1 ) Stanford machine learning
- 音乐播放器项目技术之一音乐播放进度及音量的控制
- JAVA基础学习day25--Socket基础二-多线程
- CentOS 6.X学习心得
- 提升网站性能开发的10个技巧
- java按行读取txt文件并按顺序放到map对象里面实例
- sign starfieldtech
- iOS AppDelegate浅析
- java定时器使用
- python文件操作之文件写入
- iOS限制UITextField的输入字符数