【BZOJ3481】 DZY Loves Math III（裴蜀定理、欧拉函数）

传送门

考虑如果$x$确定
那么相当于解一个模意义下的方程
则只有当$gcd(x,P)|Q$时有解
且解得个数为$gcd(x,p)$

枚举$gcd$可得
$ans=sum_{d|Q且d|P}d*sum_x[gcd(x,frac pd)=1]$
$=sum_{d|gcd(P,Q)}d*phi(frac pd)$

可以发现这个函数是一个积性函数
于是考虑每个质数$x$的贡献乘起来
若$p=x^ap',q=x^bq'$
那么此时贡献为
$sum_{i=0}^{min(a,b)}x^iphi(x^{a-i})$
$=sum_{i}x^i*(x-1)*x^{a-i-1}$
$=(x-1)min(a,b)x^{a-1}$
但如果有$i=a$时
$phi(1)=1$，所以这时候贡献是$x^a$而不是$(x-1)x^{a-1}$

用$mathrm{Pollard-Rho}$求质因子

#include<bits/stdc++.h>
using namespace std;
#define cs const
#define re register
#define pb push_back
#define pii pair<int,int>
#define ll long long
#define fi first
#define se second
#define bg begin
cs int RLEN=1<<20|1;
inline char gc(){
    static char ibuf[RLEN],*ib,*ob;
    (ib==ob)&&(ob=(ib=ibuf)+fread(ibuf,1,RLEN,stdin));
    return (ib==ob)?EOF:*ib++;
}
inline int read(){
    char ch=gc();
    int res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
inline ll readl(){
    char ch=gc();
    ll res=0;bool f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
template<class tp>inline void chemx(tp &a,tp b){a<b?a=b:0;}
template<class tp>inline void chemn(tp &a,tp b){a>b?a=b:0;}
map<ll,int> vt;
int cnt[2][150],num;
ll pr[150];
namespace Pollard_Rho{
    inline ll mul(ll a,ll b,cs ll &mod){
        return (a*b-(ll)((long double)a/mod*b)*mod+mod)%mod;
    }
    inline ll ksm(ll a,ll b,cs ll &mod){
        ll res=1;
        for(;b;b>>=1,a=mul(a,a,mod))(b&1)&&(res=mul(res,a,mod));
        return res;
    }
    inline bool check_p(ll x){
        static int pr[9]={2,3,5,7,11,13,17,19,23},l=9;
        for(int i=0;i<l;i++)if(x%pr[i]==0)return x==pr[i];
        if(x<23)return false;
        ll t=x-1,s=0;
        while(!(t&1))s++,t>>=1;
        for(int i=0;i<l;i++){
            ll b=ksm(pr[i],t,x);
            for(int j=1;j<=s;j++){
                ll k=mul(b,b,x);
                if(k==1&&b!=1&&b!=x-1)return false;
                b=k;
            }
            if(b!=1)return false;
        }
        return true;
    }
    inline ll F(ll x,ll c,ll mod){return (mul(x,x,mod)+c)%mod;}
    inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
    inline ll find_p(ll x){
        ll s=0,t=0,c=1ll*rand()*rand()%(x-1)+1;
        for(int len=1;;len<<=1,s=t){
            ll val=1;
            for(int step=1;step<=len;step++){
                t=F(t,c,x);
                val=mul(val,abs(t-s),x);
                if(!(step&127)){
                    ll d=gcd(val,x);
                    if(d>1)return d;
                }
            }
            ll d=gcd(val,x);
            if(d>1)return d;
        }
    }
    inline void ins(ll x,int kd){
        if(vt.count(x))++cnt[kd][vt[x]];
        else pr[++num]=x,vt[x]=num,cnt[kd][num]=1;
    }
    inline void rho(ll x,int kd){
        if(x<2)return;
        if(check_p(x)){ins(x,kd);return;}
        ll p=x;
        while(p>=x)p=find_p(x);
        rho(x/p,kd),rho(p,kd);
    }
}
cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?(a-mod):a;}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){static ll r;r=1ll*a*b;return (r>=mod)?(r%mod):r;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline void Dec(int &a,int b){a-=b,a+=a>>31&mod;}
inline void Mul(int &a,int b){static ll r;r=1ll*a*b;a=(r>=mod)?(r%mod):r;}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,Mul(a,a))(b&1)&&(Mul(res,a),1);return res;}
inline int Inv(int x){return ksm(x,mod-2);}
int n;ll q[12];
int main(){
    #ifdef Stargazer
    freopen("lx.in","r",stdin);
    #endif
    n=read();
    for(int i=1;i<=n;i++){
        ll x=readl();Pollard_Rho::rho(x,0);
    }bool fg=0;
    for(int i=1;i<=n;i++){
        q[i]=readl();if(q[i]==0)fg=1;
    }
    if(!fg)for(int i=1;i<=n;i++)Pollard_Rho::rho(q[i],1);
    else for(int i=1;i<=num;i++)cnt[1][i]=cnt[0][i];
    int ret=1;
    for(int i=1;i<=num;i++)if(cnt[0][i]){
        int a=cnt[0][i],b=cnt[1][i],x=pr[i]%mod;
        if(a>b)Mul(ret,mul(mul(dec(x,1),b+1),ksm(x,a-1)));
        else Mul(ret,add(mul(mul(dec(x,1),a),ksm(x,a-1)),ksm(x,a)));
    }
    cout<<ret<<'
';
}