[bzoj3481]dzy loves math III 解题报告

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

所以答案就是求∑d|(P,Q)dφ(Pd)

设P=∏pmii,(P,Q)=∏pkii,pi是质数。则答案就是求∏pmi−1i((ki+1)(pi−1)+[ki==mi])

所以假如我们可以对P分解质因子，我们就可以在O(logP)的时间复杂度内统计答案了.

所以现在我们的难点就是分解因式了！

这需要用到Miller-Rabin算法和rho算法。

rho算法基于生日悖论，用来在O(n14logn)时间内获得n的一个约数。

考虑基于一个简单的递推式xi=x2i+1 mod n和任意一个首项的数列{xn},我们可以假定它的每一项是随机分布的，而事实上它确实几乎如此，且它每一项都是由前一项决定的。这就导致根据生日悖论，它在O(n√)的时候，将出现循环节。而我们在mod n的时候，我们也保存了mod n的约数的信息。假如说n是一个合数，那么它最小的约数最大是O(n√),所以最小的循环节会在O(n14)出现。

但是如果我们暴力判断是否出现循环了的话，需要O(n142logn)的时间复杂度，所以我们使用类似维护vector一样的方法，只求（xi,x2⌊log2n⌋)，这样便可以在O(n14logn)的时间复杂度内找到一个因子。

而如果我们把n较少的质因子先除去，那么根据上面的分析，实际上rho第一个找到的因子很有可能就是n最小的质因子。

对于时间复杂度O(n14logn)中的logn，这实际上不仅仅是对于求gcd，而且还有快速乘的复杂度。

那么我们为什么需要Miller-Rabin呢？因为如果你用rho来判断n是否是质数的话，你需要O(n√)的时间，那么一个数量级为1018的质数就足以卡掉了。同样对于rho出来的因子，虽然它很有可能是一个质数，但是我们最好还是check一下，这就需要用到Miller-Rabin算法了。

Miller-Rabin算法是对费马小定理的改进。费马小定理给出了n是质数的一个充分不必要条件：an−1=a,∃a∈[1,n)，或者更强的an−1=a,∀a∈[1,n)。但是实际上，有的合数也具有这些性质，它们被称为伪素数。检验出它们就是Miller-Rabin算法的核心部分。

作为引理，我们先来看一下关于群论的部分。

考虑(S′,⊕)是群(S,⊕)的一个子群，设a∈S,A={a⊕x|x∈S′}，则称A为S’关于a的左陪集。且∀b∈A,{b⊕x|x∈S′}=A，即对于S’的任意两个左陪集A，B，要么A=B，要么A∩B=∅，即要么|A∪B|=|A|=|B|，要么|A∪B|=|A|+|B|.设S’所有左陪集的并集必为C，那么显然有|S′|∣∣|C|，且S’关于S中的任意一个元素a的左陪集中必有a，因为S’中有幺元e使得S’关于a的左陪集中有a⊕e=a，所以有C=S，所以|S′|∣∣|C|，这就是群论上的拉格朗日定理。

那么，若(S′,⊕)是群(S,⊕)的一个真子集，就显然有|S′|≤|S|2，这种放缩便可以用来分析时间复杂度了。

我们知道若x2≡1 mod n,x∈[1,n),则n|(x−1)(x+1),所以当n质数时，必然有x∈{1,n−1}。注意到这与费马小定理的关系。我们便不妨设n−1=2tu,Xa=(au mod n,a2u mod n,...,a2tu mod n),那么Xa有4

种情况(...,d),d≠1(...,d,...,1),d≠±1(...,−1,...,1)(1,...,1)我们将属于后两种情况的称为是非证据，即它们不能证明n是否为合数，而显然，如果有前两种情况出现，那么就意味着n是一个合数。

那么下面我们要说明，如果n是一个合数，则非证据数目小于等于n2。

如果n是偶数，那么显然所有偶数都是情况1，所以非证据数目≤奇数数目=n2.

如果n是奇数，再分情况讨论：

一、n不是Carmichael数，即∃x∈Z∗n，使得xn−1≠1。

那么考虑B={x∈Z∗n|xn−1≡1 mod n}。显然1∈B，所以B有幺元；考虑B关于mod n的乘法，显然有结合律；而若a,b∈B，则(ab)n−1≡1 mod n，所以ab∈B，所以有封闭律；而因为a∈Z∗n，所以a有逆元，对于an−1≡1 mod n，两边同时乘以(a−1)n，则an−1(a−1)n−1≡(a∗a−1)n−1≡1≡(a−1)n−1，所以a−1∈B.所以（B，⋅n）是Z∗n的一个子群。且x∈Z∗n，x∉B，所以B⊂Z∗n.

二、n是Carmichael数。

我们首先试着来构建B。设k为最大的j使得∃a∈Z∗n,au2j≡−1(modn),这样的(a,j)是一定存在的，因为(n−1)u≡−1(modn).令B={x|xu2k≡1 mod n}，显然所有非证据均在B内，且与情况一类似我们可以证明(B,⋅n)是Z∗n的子群，那么我们只需构造一个元素x，x∈Z∗n−B。

我们考虑中国剩余定理，那么我们需要n至少有两个质因子，而事实上正是如此。因为若n=pe(e>1)，且xpe−1≡1 mod n,xϕ(pe)≡xpe−1(p−1)≡1 mod n，就有pe−1(p−1)|pe−1,但p|pe−1(p−1),p/|pe−1，所以pe−1(p−1)/|pe−1.（事实上，更强的，n不含平方因子且至少有3个质因子）

设n=n1n2,(n1,n2)=1，a2ku≡−1 mod n.考虑b使得b≡a mod n1b≡1 mod n2，那么就有bu2k≡−1 mod n1bu2k≡1 mod n2，所以bu2k≠±1 mod n,b∈Z∗n−B，即B⊂Z∗n

综上，B⊂Z∗n,|B|≤|Z∗n|2≤n−12<n2

所以非证据数小于等于n2，证据数大于等于n2。对一个合数check x个数，出错的概率是12x。然而事实上，更强地（我并不会证明），对于一个奇合数n而言，非证据数小于等于n−14，而且确实有这样的数。

最后我们考虑一下时间复杂度，设X=1018，P的质因子有约100个，预处理10^5以内的质数，Miller-Rabin时每次随机选取6个数。

那么总时间复杂度是O(Y+n1/4logX+NZlogSXlog2X+NS)≈107…至于这是不是最优的，我感觉这起码是比较优的。。如果真要求最优时间复杂度实在太困难了。。

#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;
}