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

傳送門

考慮如果 x x x确定

那麼相當于解一個模意義下的方程

則隻有當 g c d ( x , P ) ∣ Q gcd(x,P)|Q gcd(x,P)∣Q時有解

且解得個數為 g c d ( x , p ) gcd(x,p) gcd(x,p)

枚舉 g c d gcd gcd可得

a n s = ∑ d ∣ Q 且 d ∣ P d ∗ ∑ x [ g c d ( x , p d ) = 1 ] ans=\sum_{d|Q且d|P}d*\sum_x[gcd(x,\frac pd)=1] ans=d∣Q且d∣P∑​d∗x∑​[gcd(x,dp​)=1]

= ∑ d ∣ g c d ( P , Q ) d ∗ ϕ ( p d ) =\sum_{d|gcd(P,Q)}d*\phi(\frac pd) =d∣gcd(P,Q)∑​d∗ϕ(dp​)

可以發現這個函數是一個積性函數

于是考慮每個質數 x x x的貢獻乘起來

若 p = x a p ′ , q = x b q ′ p=x^ap',q=x^bq' p=xap′,q=xbq′

那麼此時貢獻為

∑ i = 0 m i n ( a , b ) x i ϕ ( x a − i ) \sum_{i=0}^{min(a,b)}x^i\phi(x^{a-i}) i=0∑min(a,b)​xiϕ(xa−i)

= ∑ i x i ∗ ( x − 1 ) ∗ x a − i − 1 =\sum_{i}x^i*(x-1)*x^{a-i-1} =i∑​xi∗(x−1)∗xa−i−1

= ( x − 1 ) m i n ( a , b ) x a − 1 =(x-1)min(a,b)x^{a-1} =(x−1)min(a,b)xa−1

但如果有 i = a i=a i=a時

ϕ ( 1 ) = 1 \phi(1)=1 ϕ(1)=1，是以這時候貢獻是 x a x^a xa而不是 ( x − 1 ) x a − 1 (x-1)x^{a-1} (x−1)xa−1

用 P o l l a r d − R h o \mathrm{Pollard-Rho} 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<<'\n';
}