POJ——2429(數論之大整數分解)
阿新 • • 發佈:2019-01-24
解析:題目的關鍵就是把lcm/gcd的大數分解成質因子之積,用到了pollard—rho尋找因子和Miller_rabin素數檢測。(可以用容斥原理求出所有的情況,然後選出來a+b最小的情況。)
#include <iostream> #include <cmath> #include <string> #include <cstring> #include <cstdlib> #include <ctime> #include <algorithm> #include <cstdio> using namespace std; typedef long long ll; #define MAX(a,b) a>b?a:b #define MIN(a,b) a>b?b:a #define INF ((ll)1<<59) #define N 501 #define C 201 //pollard_rho裡面加的常數C #define Times 11 ll mini; int ct;//記錄因子的個數 ll r[N]; //儲存N的質因子 ll gcd(ll a,ll b){return b?gcd(b,a%b):a;} //歐幾里得 ll random(ll n) //生成隨機數 { return (ll)((double)rand()/RAND_MAX*n+0.5); } ll mul(ll a,ll b,ll MOD) //a*b%MOD { ll ans=0; while(b>0) { if(b&1) { ans=(ans+a)%MOD; b--; } b=b>>1; a=(a<<1)%MOD; } return ans; } ll pow_mod_1(ll x,ll y,ll MOD) //快速冪mod(非遞迴實現) { ll ans=1; x%=MOD; while(y>0) //x,y的值可能會很大,因此選擇這個可能會較好。 { if(y&1) { ans=mul(ans,x,MOD); y--; } y=y>>1; x=mul(x,x,MOD); } return ans; } bool Witness(ll a,ll n) //Miller_rabin,素數測試 { ll m=n-1; int j=0; while(!(m&1)) { j++; m=m>>1; } ll x=pow_mod_1(a,m,n); if(x==1||x==n-1) return false; // while(j--) { x=x*x%n; if(x==n-1) return false; //若x==n-1,則x*x也為n-1 } return true; } bool Miller_rabin(ll n) { if(n<2) return false; if(n==2) return true; if(!(n&1)) return false; for(int i=1;i<=Times;i++) //測試了Times次,錯誤概率降到很低,幾乎為0 { ll a=random(n-2)+1; if(Witness(a,n)) return false; } return true; } ll pollard_rho(ll n,int c) //尋找N的因子。 { ll x,y,d,i=1,k=2; x=random(n-1)+1; y=x; while(1) { i++; x=(mul(x,x,n)+c)%n; d=gcd(y-x,n); if(d>1&&d<n) return d; if(y==x) return n; if(i==k) { y=x; k=k<<1; } } } void find(ll n,int k) //找到因子,為素數則儲存。 { if(n==1) return ; if(Miller_rabin(n)) { r[ct++] =n; if(n<mini) mini=n; return ; } ll p=n; while(p>=n) p=pollard_rho(p,k--); //p為n的因子 find(p,k); find(n/p,k); } int temp[N][2]; int q; //記錄含有不同質因子的個數,為容斥原理做鋪墊。 void settle() //整理得到的質因子到temp中,temp[N][0]記錄的是質因子,temp[N][1]記錄的是該質因子的個數 { sort(r,r+ct); q=0; temp[q][0]=r[0]; temp[q][1]=1; for(ll i=1;i<ct;i++) { if(r[i]==r[i-1]) temp[q][1]++; else { q++; temp[q][0]=r[i]; temp[q][1]=1; } } q++; //符合習慣 } ll Pow(ll x,int y) //x^y { if(y==1) return x; ll ans=Pow(x,y/2); return y%2?(ans*ans)*x:(ans*ans); } void swap(ll &a,ll &b) //交換 { ll t=a; a=b; b=t; } void resolve(ll lcm_gcd,ll gcdd) { ll mini=INF; ll ed=1<<q; ll temp_p; ll ans_a,ans_b; ll mid_l; ll mid_r; for(ll i=1;i<ed-1;i++) //容斥,見有人是dfs實現的,好像複雜度差不多。 { temp_p=i; mid_l=1; int k=0; while(temp_p) { if(temp_p&1) mid_l*=Pow(temp[k][0],temp[k][1]); k++; temp_p=temp_p>>1; } ll mid_r=lcm_gcd/mid_l+mid_l; //(a+b)兩者之和 if(mid_r<mini) { mini=mid_r; ans_a=mid_l; ans_b=mid_r-mid_l; } } if(ans_b<ans_a) swap(ans_b,ans_a);//調整下。 cout<<ans_a*gcdd<<" "<<ans_b*gcdd<<endl; } int main() { int i,j; ll lcmm,gcdd; ll lcm_gcd; while(cin>>gcdd>>lcmm) { if(gcdd==lcmm) { cout<<gcdd<<" "<<lcmm<<endl; continue; } ct=0; lcm_gcd=lcmm/gcdd; find(lcm_gcd,C); settle(); resolve(lcm_gcd,gcdd); } return 0; }