NJUST1722(大数分解的应用)

最新推荐文章于 2018-04-11 21:13:35 发布

转载最新推荐文章于 2018-04-11 21:13:35 发布 · 916 阅读

算法专栏收录该内容

538 篇文章

订阅专栏

本文介绍了一道关于大数分解与平方差的应用题，通过解析算法流程，包括快速幂运算、Miller-Rabin素性测试、Pollard's rho质因数分解等，解决了求特定形式平方差的问题。

NJUST1722(大数分解的应用)

分类：数论 2013-08-01 20:40 125人阅读评论(0) 收藏举报

题目：所有的平方差

题意：n = a1*a1 - b1*b1, n = a2*a2 - b2*b2, ... , n = am*am - bm*bm

其中对于任意的1<=k1,k2<=m，ak1 != ak2 或者 bk1 != bk2，并且ak1,ak2,bk1,bk2>=0

那么我们需要求的是n^(a1*a1+b1*b1+a2*a2+b2*b2+...+am*am+bm*bm)%大质数10000000019的值。

本题注意一点，由于10000000019很大，就是最后求出来的ak与ak在相乘的时候注意一下，比如ak大于10^10时，直接相乘

就会爆LL，我们这里有技巧，就是把乘法变为加法，二分加法，也就是下面的multi函数。

通过本题学到了一招，就是对于求x*x%MOD怎么处理，x与MOD大于10^10。

[cpp]view plaincopy 
    
 #include <stdio.h>  
 #include <stdlib.h>  
 #include <string.h>  
 #include <algorithm>  
 #include <iostream>  
 #include <math.h>  
   
 using namespace std;  
 typedef unsigned long long LL;  
   
 const LL Times=10;  
 const LL N=555;  
 const LL MOD=10000000018;  
   
 LL ct,cnt,c,n;  
 LL fac[N],num[N];  
 LL arr[N];  
   
 LL gcd(LL a,LL b)  
 {  
     return b? gcd(b,a%b):a;  
 }  
   
 LL multi(LL a,LL b,LL m)  
 {  
     LL ans=0;  
     while(b)  
     {  
         if(b&1)  
         {  
             ans=(ans+a)%m;  
             b--;  
         }  
         b>>=1;  
         a=(a+a)%m;  
     }  
     return ans;  
 }  
   
 LL quick_mod(LL a,LL b,LL m)  
 {  
     LL ans=1;  
     a%=m;  
     while(b)  
     {  
         if(b&1)  
         {  
             ans=multi(ans,a,m);  
             b--;  
         }  
         b>>=1;  
         a=multi(a,a,m);  
     }  
     return ans;  
 }  
   
 bool Miller_Rabin(LL n)  
 {  
     if(n==2) return true;  
     if(n<2||!(n&1)) return false;  
     LL a,m=n-1,x,y;  
     int k=0;  
     while((m&1)==0)  
     {  
         k++;  
         m>>=1;  
     }  
     for(LL i=0; i<Times; i++)  
     {  
         a=rand()%(n-1)+1;  
         x=quick_mod(a,m,n);  
         for(LL j=0; j<k; j++)  
         {  
             y=multi(x,x,n);  
             if(y==1&&x!=1&&x!=n-1) return false;  
             x=y;  
         }  
         if(y!=1) return false;  
     }  
     return true;  
 }  
   
 LL Pollard_rho(LL n,LL c)  
 {  
     LL x,y,d,i=1,k=2;  
     y=x=rand()%(n-1)+1;  
     while(true)  
     {  
         i++;  
         x=(multi(x,x,n)+c)%n;  
         d=gcd((y-x+n)%n,n);  
         if(1<d&&d<n) return d;  
         if(y==x) return n;  
         if(i==k)  
         {  
             y=x;  
             k<<=1;  
         }  
     }  
 }  
   
 void find(LL n,LL c)  
 {  
     if(n==1) return;  
     if(Miller_Rabin(n))  
     {  
         fac[ct++]=n;  
         return ;  
     }  
     LL p=n;  
     LL k=c;  
     while(p>=n) p=Pollard_rho(p,c--);  
     find(p,k);  
     find(n/p,k);  
 }  
   
 void dfs(LL dept,LL product=1)  
 {  
     if(dept==cnt)  
     {  
         if(product<=(LL)sqrt(1.0*n)&&(product+n/product)%2==0&&(n/product-product)%2==0)  
             arr[c++]=product;  
         return;  
     }  
     for(LL i=0; i<=num[dept]; i++)  
     {  
         dfs(dept+1,product);  
         product*=fac[dept];  
     }  
 }  
   
 int main()  
 {  
     LL t,tt=1;  
     LL ans,x,y;  
     cin>>t;  
     while(t--)  
     {  
         cin>>n;  
   
         printf("Case %d: ",tt++);  
         if(n==1)  
         {  
             puts("1");  
             continue;  
         }  
         ct=0;  
         c=0;  
         find(n,120);  
         sort(fac,fac+ct);  
         num[0]=1;  
         LL k=1;  
         for(LL i=1; i<ct; i++)  
         {  
             if(fac[i]==fac[i-1])  
                 ++num[k-1];  
             else  
             {  
                 num[k]=1;  
                 fac[k++]=fac[i];  
             }  
         }  
         cnt=k;  
         dfs(0,1);  
         sort(arr,arr+c);  
         ans=0;  
         for(LL i=0; i<c; i++)  
         {  
             x=(n/arr[i]+arr[i])/2;  
             y=(n/arr[i]-arr[i])/2;  
             ans+=multi(x,x,MOD)+multi(y,y,MOD);  
             ans%=MOD;  
         }  
         if(c>0) cout<<quick_mod(n,ans,MOD+1)<<endl;  
         else    puts("-1");  
     }  
     return 0;  
 }  