题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3481
推式子:xy % P = Q 的个数
由于 0 <= x,y < P,所以对于一个固定的 x,如果 (x,P) | Q,则有 (x,P) 个解;
所以个数为 ∑(0 <= x < P ) (x,P) * [ (x,P) | Q ] ( [...] 表示 ... 为真则为1,否则为0)
= ∑( d|P, d|Q ) d * φ( P/d )
令 Q = (P,Q),则
= ∑( d|Q ) d * φ( P/d )
单独考虑每个质因子的贡献,若P = p ^ t * p ^ t2 p ^ t,Q = p ^ k * p ^ k * ... * p ^ k,则原式
= ∏( 1 <= i <= m ) ∑( 0 <= x <= k[i] ) p * φ(p^t-x)
而 φ(p^t) = p^t * ( p - 1 ) / p
所以原式最终变成 ∏( 1 <= i <= m ) p^(t-1) * ( p - 1 ) * ( k + 1 )
注意 t[i] = k[i] 时略有不同,求和的那一项 = p , 所以如果有 t[i] = k[i],则再加一个 p^(t-1);
然后...质因数分解要用 Pollar_rho,Pollard_rho 内又用到 Miller-Rabin,参考了一下博客:
https://blog.csdn.net/sunshine_cfbsl/article/details/52512706
https://www.cnblogs.com/kuangbin/archive/2012/08/19/2646396.html
https://blog.csdn.net/forever_shi/article/details/80547041
套上板子就可以了;
但要注意,如果把 P 和 Q 的所有质因子都提取出来会 TLE,因为令 Q = (P,Q),所以只有 Q 有的质因子就不用专门记录下来了;
直接 rand() 似乎只能得到 int,所以可以专门写一下 ll 的,虽然只有 int 好像也可以;
注意各处都是 ll 。
代码如下:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<ctime>
using namespace std;
typedef long long ll;
int const xn=;
int m,S=,t[xn],k[xn];
ll pri[xn],ans;
ll rd()
{
ll ret=,f=; char ch=getchar();
while(ch<''||ch>''){if(ch=='-')f=; ch=getchar();}
while(ch>=''&&ch<='')ret=(ret<<3ll)+(ret<<1ll)+ch-'',ch=getchar();
return f?ret:-ret;
}
ll Rand(ll n)
{
return (((ll)rand()<<)|rand())%(n-)+;
// return rand()%(n-1)+1;//也可,但无ll
}
ll upt(ll x,ll mod){while(x>=mod)x-=mod; while(x<)x+=mod; return x;}
ll mul(ll a,ll b,ll mod)
{
ll ret=; a=a%mod; b=b%mod;
for(;b;b>>=1ll,a=upt(a+a,mod))
if(b&)ret=upt(ret+a,mod);
return ret;
}
ll pw(ll a,ll b,ll mod)
{
ll ret=; a=a%mod;
for(;b;b>>=1ll,a=mul(a,a,mod))
if(b&)ret=mul(ret,a,mod);
return ret;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
bool ck(ll a,ll u,ll t,ll p)
{
ll x=pw(a,u,p),lst=x;
for(int i=;i<=t;i++)
{
x=mul(x,x,p);
if(x==&&lst!=&&lst!=p-)return ; // a^2 % p = 1 -> a=1,a=p-1
lst=x;
}
if(x!=)return ; // a^p-1 % p = 1
return ;
}
bool MR(ll x)
{
if(x<=||x&==)return ;
if(x==)return ;
ll t=,u=x-;
while(u&==)u>>=1ll,t++;
for(int i=;i<=S;i++)
{
ll a=Rand(x)%(x-)+;
if(!ck(a,u,t,x))return ;
}
return ;
}
ll rho(ll n,ll c)
{
ll x=Rand(n)%n,y=x;
ll i=,k=;
while()
{
i++;
x=(mul(x,x,n)+c)%n;
ll d=gcd((x-y+n)%n,n);
if(d!=&&d!=n)return d;
if(x==y)return n;
if(i==k)y=x,k+=k;
}
}
void add(ll x,bool tp)//1:P-t[]
{
for(int i=;i<=m;i++)
if(pri[i]==x)
{
if(tp)t[i]++; else k[i]++;
return;
}
if(!tp)return;//否则TLE!!!
pri[++m]=x;
if(tp)t[m]++; else k[m]++;
}
void find(ll x,bool tp)
{
if(x==)return;//!
if(MR(x)){add(x,tp); return;}
ll p=x;
while(p>=x)p=rho(p,Rand(x)%(x-)+);
find(p,tp); find(x/p,tp);
}
int main()
{
srand();//或 time(0);
int mod=1e9+;
int n=rd(); ll x; bool fl=;
for(int i=;i<=n;i++)x=rd(),find(x,);//1:P-t[]
for(int i=;i<=n;i++)
{
x=rd();
if(!x){fl=; continue;}
else find(x,);
}
if(fl)for(int i=;i<=m;i++)k[i]=t[i];
else for(int i=;i<=m;i++)k[i]=min(k[i],t[i]);//(P,Q)
ans=;
for(int i=;i<=m;i++)
{
ll pp=pw(pri[i],t[i]-,mod);
ll mm=mul(pri[i]-,k[i]+,mod);
if(k[i]==t[i])mm++;
ans=mul(ans,mul(pp,mm,mod),mod);
}
printf("%lld\n",ans);
return ;
}