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

6

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=100010;
LL read(){
    LL x=0,f=1;char ch=getchar();
    while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while(ch>='0' && ch<='9'){x=x*10-'0'+ch;ch=getchar();}
    return x*f;
}
LL Rand(LL n){return (((LL)rand()<<31)|rand())%(n-1)+1;}
//
map<LL,int>primap;
LL pri[mxn*10];int cnt=0;
bool vis[mxn];
void init(){
    int mxn=200;
    for(int i=2;i<mxn;i++){
        if(!vis[i]){
            pri[++cnt]=i;
            primap[i]=cnt;
        }
        for(int j=1;j<=cnt && pri[j]*i<mxn;j++){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0)break;
        }
    }
    return;
}
//
LL mul(LL x,LL k,LL P)
{
    LL res=0;
    while(k){
        if(k&1)res=(res+x)%P;
        x=(x+x)%P;
        k>>=1;
    }
    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=1;
    while(k){
        if(k&1)res=mul(res,a,P);
        a=mul(a,a,P);
        k>>=1;
    }
    return res;
}
///
bool check(LL a,LL n,LL r,LL s){
    LL res=ksm(a,r,n);LL b=res;
    for(int i=1;i<=s;i++){
        res=mul(res,res,n);
        if(res==1 && b!=1 && b!=n-1)return 1;
        b=res;
    }
    return (res!=1);
}
bool MR(LL n){//素数测试 
    if(n<=1)return 0;
    if(n==2)return 1;
    if(~n&1)return 0;
    LL r=n-1,s=0;
    while(!(r&1))r>>=1,s++;
    for(int i=1;i<=10;i++)
        if(check(Rand(n),n,r,s))return 0;
    return 1;
}
///
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]=1;
    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<0)while(1);
    LL k=2,x=Rand(n),y=x,p=1;
    for(LL i=1;p==1;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<<=1;
    }
    return p;
}
void solve(LL x,bool flag){//分解质因数 
    if(x==1)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+7;
int n;
LL P[50],Q[50];
void Calc(){
    LL ans=1;
    for(int i=1;i<=cnt;i++){
        if(!p[i])continue;
        LL R=mul((q[i]+1),(pri[i]-1),mod);
        if(p[i]==q[i])R++;
        R=mul(R,ksm(pri[i],p[i]-1,mod),mod);
        ans=mul(ans,R,mod);
    }
    printf("%lld\n",ans);
    return;
}
int main(){
    srand(19260817);
    init();
    int i,j;
    n=read();
    for(i=1;i<=n;i++){
        P[i]=read();
        solve(P[i],0);
    }
    for(i=1;i<=n;i++){
        Q[i]=read();
        if(!Q[i]){//特判0 
            for(j=1;j<=cnt;j++)q[j]=p[j];
            break;
        }
        else solve(Q[i],1);
    }
    Calc();
    return 0;
}