[bzoj3481]dzy loves math III 解题报告

显然，(Q mod P,P)=(P,Q)，那么我们考虑x,(x,P)=d,d|(P,Q)，显然它的加法运算循环节是$P \over d$，所以对于y它有${P \over {P \over d}}=d$种选择使得$xy \equiv Q(mod \ P)$,而这样的x有$φ({P \over d})$个。而如果(x,P)不能整除(P,Q)，那么显然不存在y使得$xy \equiv Q(mod P)$，因为如果不这样的话，就以为着(Q mod P,P)>(P,Q).
所以答案就是求

$$\sum_{d|(P,Q)}dφ({P \over d})$$
设P=$\prod p_i^{m_i},(P,Q)=\prod p_i^{k_i}$,$p_i$是质数。则答案就是求 $$\prod p_i^{m_i-1}((k_i+1)(p_i-1)+[k_i==m_i])$$
所以假如我们可以对P分解质因子，我们就可以在$O(logP)$的时间复杂度内统计答案了.

所以现在我们的难点就是分解因式了！
这需要用到Miller-Rabin算法和rho算法。

rho算法基于生日悖论，用来在$O(n^{1 \over 4}logn)$时间内获得n的一个约数。
考虑基于一个简单的递推式$x_i=x_i^2+1 \ mod \ n$和任意一个首项的数列$\{ x_n \}$,我们可以假定它的每一项是随机分布的，而事实上它确实几乎如此，且它每一项都是由前一项决定的。这就导致根据生日悖论，它在$O(\sqrt n)$的时候，将出现循环节。而我们在mod n的时候，我们也保存了mod n的约数的信息。假如说n是一个合数，那么它最小的约数最大是$O(\sqrt n)$,所以最小的循环节会在$O(n^{1\over 4})$出现。
但是如果我们暴力判断是否出现循环了的话，需要$O({n^{1\over 4}}^2logn)$的时间复杂度，所以我们使用类似维护vector一样的方法，只求$（x_i,x_{2^{\lfloor log_2n\rfloor}})$，这样便可以在$O(n^{1\over 4}logn)$的时间复杂度内找到一个因子。
而如果我们把n较少的质因子先除去，那么根据上面的分析，实际上rho第一个找到的因子很有可能就是n最小的质因子。
对于时间复杂度$O(n^{1\over 4}logn)$中的$logn$，这实际上不仅仅是对于求gcd，而且还有快速乘的复杂度。

那么我们为什么需要Miller-Rabin呢？因为如果你用rho来判断n是否是质数的话，你需要$O(\sqrt n)$的时间，那么一个数量级为$10^{18}$的质数就足以卡掉了。同样对于rho出来的因子，虽然它很有可能是一个质数，但是我们最好还是check一下，这就需要用到Miller-Rabin算法了。
Miller-Rabin算法是对费马小定理的改进。费马小定理给出了n是质数的一个充分不必要条件：$a^{n-1}=a,\exists a \in [1,n)$，或者更强的$a^{n-1}=a,\forall a \in [1,n)$。但是实际上，有的合数也具有这些性质，它们被称为伪素数。检验出它们就是Miller-Rabin算法的核心部分。
作为引理，我们先来看一下关于群论的部分。
考虑$(S',\oplus )$是群$(S,\oplus )$的一个子群，设$a \in S,A=\{ a \oplus x|x \in S' \}$，则称A为S’关于a的左陪集。且$\forall b \in A,\{ b \oplus x|x \in S'\} =A$，即对于S’的任意两个左陪集A，B，要么A=B，要么$A \cap B=\varnothing$，即要么$|A \cup B|=|A|=|B|$，要么$|A \cup B|=|A|+|B|$.设S’所有左陪集的并集必为C，那么显然有$|S'| \big | |C|$，且S’关于S中的任意一个元素a的左陪集中必有a，因为S’中有幺元e使得S’关于a的左陪集中有$a \oplus e=a$，所以有C=S，所以$|S'| \big | |C|$，这就是群论上的拉格朗日定理。
那么，若$(S',\oplus)$是群$(S,\oplus )$的一个真子集，就显然有$|S'| \le {|S| \over 2}$，这种放缩便可以用来分析时间复杂度了。

我们知道若$x^2 \equiv 1\ mod \ n,x \in [1,n)$,则$n|(x-1)(x+1)$,所以当n质数时，必然有$x \in \{ 1,n-1 \}$。注意到这与费马小定理的关系。我们便不妨设$n-1=2^tu,X_a=(a^u\ mod \ n,a^2u\ mod \ n,...,a^{2^tu}\ mod \ n)$,那么$X_a$有4
种情况

$$(...,d),d \not= 1$$ $$(...,d,...,1),d \not= \pm 1$$ $$(...,-1,...,1)$$ $$(1,...,1)$$ 我们将属于后两种情况的称为是非证据，即它们不能证明n是否为合数，而显然，如果有前两种情况出现，那么就意味着n是一个合数。
那么下面我们要说明，如果n是一个合数，则非证据数目小于等于$n \over 2$。
如果n是偶数，那么显然所有偶数都是情况1，所以非证据数目$\le$奇数数目=$n \over 2$.
如果n是奇数，再分情况讨论：
一、n不是Carmichael数，即$\exists x \in Z_n^*$，使得$x^{n-1} \not= 1$。
那么考虑$B=\{ x \in Z_n^*|x^{n-1} \equiv 1 \ mod \ n \}$。显然1∈B，所以B有幺元；考虑B关于mod n的乘法，显然有结合律；而若a,b∈B，则$(ab)^{n-1} \equiv 1 \ mod \ n$，所以ab∈B，所以有封闭律；而因为$a\in Z_n^*$，所以a有逆元，对于$a^{n-1} \equiv 1 \ mod \ n$，两边同时乘以$(a^{-1})^n$，则$a^{n-1} (a^{-1})^{n-1} \equiv (a*a^{-1})^{n-1} \equiv 1 \equiv (a^{-1})^{n-1}$，所以$a^{-1} \in B$.所以$（B，·_n）$是$Z_n^*$的一个子群。且$x \in Z_n^*$，$x \notin B$，所以$B \subset Z_n^*$.
二、n是Carmichael数。
我们首先试着来构建B。设k为最大的j使得$\exists a \in Z_n^*,a^{u2^j} \equiv -1(mod n)$,这样的(a,j)是一定存在的，因为$(n-1)^u \equiv -1 (mod n)$.令$B=\{ x|x^{u2^k} \equiv 1 \ mod \ n \}$，显然所有非证据均在B内，且与情况一类似我们可以证明$(B,·_n)$是$Z_n^*$的子群，那么我们只需构造一个元素x，$x \in Z_n^*-B$。
我们考虑中国剩余定理，那么我们需要n至少有两个质因子，而事实上正是如此。因为若$n=p^e(e>1)$，且$x^{p^e-1} \equiv 1 \ mod \ n,x^{\phi (p^e)} \equiv x^{p^{e-1}(p-1)} \equiv 1 \ mod \ n$，就有$p^{e-1}(p-1)|p^e-1$,但$p|p^{e-1}(p-1),p\not| p^e-1$，所以$p^{e-1}(p-1) \not| p^e-1$.（事实上，更强的，n不含平方因子且至少有3个质因子）
设$n=n_1n_2,(n_1,n_2)=1$，$a^{2^ku}\equiv -1 \ mod \ n$.考虑b使得 $$b \equiv a \ mod \ n_1$$ $$b \equiv 1 \ mod \ n_2$$ ，那么就有 $$b^{u2^k} \equiv -1 \ mod \ n_1$$ $$b^{u2^k} \equiv 1 \ mod \ n_2$$ ，所以$b^{u2^k} \not= \pm 1\ mod \ n,b \in Z_n^*-B$，即$B \subset Z_n^*$
综上，$B \subset Z_n^*,|B| \le {|Z_n^*|\over 2}\le {n-1 \over 2}<{n \over 2}$
所以非证据数小于等于$n \over 2$，证据数大于等于$n \over 2$。对一个合数check x个数，出错的概率是$1 \over 2^x$。然而事实上，更强地（我并不会证明），对于一个奇合数n而言，非证据数小于等于$n-1 \over 4$，而且确实有这样的数。

最后我们考虑一下时间复杂度，设$X=10^{18}$，P的质因子有约100个，预处理10^5以内的质数，Miller-Rabin时每次随机选取6个数。
那么总时间复杂度是$O(Y+n^{1/4}logX+NZlog_SXlog^2X+NS)≈10^7$…至于这是不是最优的，我感觉这起码是比较优的。。如果真要求最优时间复杂度实在太困难了。。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
using namespace std;
#include<cstdlib>
int prm[100005];
bool is_prm[100005];
typedef long long LL;
unsigned rand_32(){
    return rand()<<17^rand()<<5^rand();
}
LL rand_63(){
    return (LL)rand_32()<<31^rand_32();
}
LL multi(LL a,LL b,LL Mod){
    LL ans=0,sum=a%Mod;
    for(;b;b>>=1,sum=(sum<<1)%Mod)
        if(b&1)
            ans=(ans+sum)%Mod;
    return ans;
}
LL power(LL a,LL b,LL Mod){
    LL ans=1,pro=a%Mod;
    for(;b;b>>=1,pro=multi(pro,pro,Mod))
        if(b&1)
            ans=multi(ans,pro,Mod);
    return ans;
}
bool mr(LL x){
    LL pre,now,tmp,cnt;
    tmp=x-1;
    while(~tmp&1)tmp>>=1;
    for(int i=6;i--;){
        pre=now=power(rand_63()%x+1,tmp,x);
        if(pre==x-1||pre==1)continue;
        for(cnt=tmp<<1;cnt<=x-1;cnt<<=1){
            now=multi(pre,pre,x);
            if(now==1)
                if(pre!=1&&pre!=x-1)return 0;
                else break;
            pre=now;
        }
        if(now!=1)return 0;
    }
    return 1;
}
LL gcd(LL a,LL b){
    return b?gcd(b,a%b):a;
}
LL dvs[205];
int tot,Pcnt[205],Gcnt[205];
void add(LL g,LL &n){
    //cout<<g<<" "<<n<<endl;
    for(dvs[++tot]=g;n%g==0;n/=g)++Pcnt[tot];
    //cout<<g<<" "<<n<<endl;
}
void rho(LL &n){
    //cout<<"->";
    if(mr(n)){
        add(n,n);
        return;
    }
    LL rcd=rand_63()%n,x=rcd,g;
    for(int i=1,k=1;;++i){
        if((g=gcd((((x=(multi(x,x,n)+1)%n)-rcd)%n+n)%n,n))!=1){
            if(mr(g))add(g,n);
            //cout<<"<-\n";
            return;
        }
        if(i==k){
            k<<=1;
            rcd=x;
        }
    }
    //cout<<"<-\n";
}
#define Mod 1000000007
int main(){
    freopen("bzoj_3481.in","r",stdin);
    freopen("bzoj_3481_std.out","w",stdout);
    int i,j;
    for(i=2;i<=100;++i){
        if(!is_prm[i])prm[++prm[0]]=i;
        for(j=1;j<=prm[0]&&i*prm[j]<=100;++j){
            is_prm[i*prm[j]]=1;
            if(i%prm[j]==0)break;
        }
    }
    LL p,q;
    int n;
    LL root;
    cin>>n;
    for(j=n;j--;){
        cin>>p;
        for(i=tot;i;--i)
            for(;p%dvs[i]==0;p/=dvs[i])
                ++Pcnt[i];
        for(i=1,root=sqrt(p);i<=prm[0]&&prm[i]<=root;++i)
            if(p%prm[i]==0)
                add(prm[i],p);
        if(prm[i]>root&&p!=1)add(p,p);
        while(p!=1)rho(p);
    }
    for(j=n;j--;){
        cin>>q;
        if(q==0){
            for(i=tot;i;--i)Gcnt[i]=Pcnt[i];
            break;
        }
        //cout<<"---"<<q<<"--\n";
        for(i=tot;i;--i)
            for(;q%dvs[i]==0;q/=dvs[i]){
                ++Gcnt[i];
                //cout<<"Get "<<q<<" "<<dvs[i]<<endl;
            }
    }
    LL ans=1,pwr;
    for(i=tot;i;--i){
        Gcnt[i]=min(Gcnt[i],Pcnt[i]);
        pwr=1,dvs[i]%=Mod;
        for(j=Pcnt[i]-1;j>0;--j)pwr=pwr*dvs[i]%Mod;
        ans=ans*((Gcnt[i]+1)*(dvs[i]-1)%Mod+(Gcnt[i]==Pcnt[i]))*pwr%Mod;
    }
    cout<<(ans+Mod)%Mod;
}