Bzoj3481 DZY Loves Math III

Time Limit: 5 Sec Memory Limit: 64 MB
Submit: 310 Solved: 65

Description

Input

Output

Sample Input

3
1
2
3
2
4
2

Sample Output

HINT

1<=N<=10,0<=Qi<=10^18,1<=Pi<=10^18,P>=2

本题仅四组数据。

Source

By Jc

数学问题欧拉函数 Miller-Rabin Pollard-rho

花了一整晚来怼这道题……在long long的领域遇到了许多问题。

假设我们有充足的时间枚举每一个x，那么在x确定的情况下，原式变成了一个模方程。
根据裴蜀定理，我们知道只有当$ gcd(x,P)|Q $ 的情况下方程有 $ gcd(x,P) $ 个解。
现在我们可以愉快地枚举每一个gcd，那么答案就是：
$$ans=\sum_{d|P,d|Q} d * \sum_{x}[gcd(x,\frac{P}{d})==1]$$
后面那个sum显然是欧拉函数
那么$$ ans=\sum_{d|P,d|Q} d · \varphi(\frac{P}{d}) $$
分(an)析(zhao)一(tao)波(lu)，发现这是一个积性函数，所以我们可以分别考虑每个质因子的贡献，再把它们累乘起来。
我们现在考虑一个质因子p，它在P中有$q$个，在Q中有$q'$个：
它对答案的贡献是
$$ \sum_{i=0}^{q'} p^i * \varphi(p^{q-i})$$
我们知道
$$\varphi(p^{q})=p^{q}·\frac{p-1}{p}$$
所以上面的sum可以化成：
$$p^{q-1}·(p-1)·(q'+1)$$
但是有一个特殊情况，当i=q的时候,$\varphi(1)=1$，不能表示成$\frac{p-1}{p}$的形式，而我们却把它用这种形式算进去了。
也就是说我们把一个$p^q$ 的贡献算成了 $(p-1)p^{q-1}$，特判一下消除即可。

 #include<iostream>

 #include<algorithm>

 #include<cstring>

 #include<cstdio>

 #include<cmath>

 #include<map>

 #define LL long long

 using namespace std;

 const int mxn=;

 LL read(){

     LL x=,f=;char ch=getchar();

     while(ch<'' || ch>''){if(ch=='-')f=-;ch=getchar();}

     while(ch>='' && ch<=''){x=x*-''+ch;ch=getchar();}

     return x*f;

 }

 LL Rand(LL n){return (((LL)rand()<<)|rand())%(n-)+;}

 //

 map<LL,int>primap;

 LL pri[mxn*];int cnt=;

 bool vis[mxn];

 void init(){

     int mxn=;

     for(int i=;i<mxn;i++){

         if(!vis[i]){

             pri[++cnt]=i;

             primap[i]=cnt;

         }

         for(int j=;j<=cnt && pri[j]*i<mxn;j++){

             vis[pri[j]*i]=;

             if(i%pri[j]==)break;

         }

     }

     return;

 }

 //

 LL mul(LL x,LL k,LL P)

 {

     LL res=;

     while(k){

         if(k&)res=(res+x)%P;

         x=(x+x)%P;

         k>>=;

     }

     return res;

 }

 /*

 LL mul(LL a,LL b,LL P){

     LL d=(long double)a/P*b+1e-8;

     a=a*b-d*P;

     a=(a<0)?a+P:a;

     printf("%lld \n",a);

     return a;

 }*/

 LL ksm(LL a,LL k,LL P){

     LL res=;

     while(k){

         if(k&)res=mul(res,a,P);

         a=mul(a,a,P);

         k>>=;

     }

     return res;

 }

 ///

 bool check(LL a,LL n,LL r,LL s){

     LL res=ksm(a,r,n);LL b=res;

     for(int i=;i<=s;i++){

         res=mul(res,res,n);

         if(res== && b!= && b!=n-)return ;

         b=res;

     }

     return (res!=);

 }

 bool MR(LL n){//素数测试

     if(n<=)return ;

     if(n==)return ;

     if(~n&)return ;

     LL r=n-,s=;

     while(!(r&))r>>=,s++;

     for(int i=;i<=;i++)

         if(check(Rand(n),n,r,s))return ;

     return ;

 }

 ///

 LL p[mxn],q[mxn];

 void addpri_P(LL x){

     if(primap.count(x)){

         ++p[primap[x]];return;

     }

     pri[++cnt]=x;

     primap[x]=cnt;

     p[cnt]=;

     return;

 }

 void addpri_Q(LL x){

     if(primap.count(x)){

         int t=primap[x];

         if(q[t]<p[t])++q[t];

     }

     return;

 }

 LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}

 LL Rho(LL n,LL c){

     if(n<)while();

     LL k=,x=Rand(n),y=x,p=;

     for(LL i=;p==;i++){

         x=(mul(x,x,n)+c)%n;

         p=(y>x)?(y-x):(x-y);

         p=gcd(n,p);

         if(i==k)y=x,k<<=;

     }

     return p;

 }

 void solve(LL x,bool flag){//分解质因数

     if(x==)return;

     if(MR(x)){

         if(!flag)addpri_P(x);//P

         else addpri_Q(x);//Q

         return;

     }

     LL t=x;

     while(t==x)t=Rho(x,Rand(x));

     solve(t,flag);

     solve(x/t,flag);

     return;

 }

 //

 const int mod=1e9+;

 int n;

 LL P[],Q[];

 void Calc(){

     LL ans=;

     for(int i=;i<=cnt;i++){

         if(!p[i])continue;

         LL R=mul((q[i]+),(pri[i]-),mod);

         if(p[i]==q[i])R++;

         R=mul(R,ksm(pri[i],p[i]-,mod),mod);

         ans=mul(ans,R,mod);

     }

     printf("%lld\n",ans);

     return;

 }

 int main(){

     srand();

     init();

     int i,j;

     n=read();

     for(i=;i<=n;i++){

         P[i]=read();

         solve(P[i],);

     }

     for(i=;i<=n;i++){

         Q[i]=read();

         if(!Q[i]){//特判0

             for(j=;j<=cnt;j++)q[j]=p[j];

             break;

         }

         else solve(Q[i],);

     }

     Calc();

     return ;

 }

巴特西