bzoj 3481 DZY Loves Math III——反演+rho分解质因数

时间:2023-03-09 14:47:43
bzoj 3481 DZY Loves Math III——反演+rho分解质因数

题目:https://www.lydsy.com/JudgeOnline/problem.php?id=3481

推推式子发现:令Q=gcd(P,Q),ans=Σ(d|Q) d*phi(P/d)。把 d 质因数分解,设 t 为 Q 的指数, w 为 P 的指数,ans变成每个质数的 Σ(i=0~t) p^i * phi( p^(w-i) ) 连乘。

分解质因数用 Pollar Rho 。

注意 Q=0 就是 Q=P,要特判!而且不要以为答案变成  (!x || !y) 了!

d从0到P-1 就是 d从1到P!不要特判 P==Q时给答案减P !因为算的时候就没算d=0的!

每个质数的那个式子可以化简,把 phi 写开,和 p^i 合并,就不用枚举 i 。但要注意 w-i ==0 时 phi 的式子有些不同了。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#define ll long long
using namespace std;
const int N=,mod=1e9+;
int n,ans=,c[N],c2[N],tot;
ll p[N],P=;
bool vis[N],flag;
ll rdl()
{
ll ret=;bool fx=;char ch=getchar();
while(ch>''||ch<''){if(ch=='-')fx=;ch=getchar();}
while(ch>=''&&ch<='') ret=(ret<<3ll)+(ret<<1ll)+ch-'',ch=getchar();
return fx?ret:-ret;
}
void upd(ll &x,ll md){x-=(x>=md?md:);}
ll mul(ll a,ll b,ll md)
{
ll ret=;//0!
while(b)
{
if(b&1ll)ret=ret+a,upd(ret,md);
a=a+a,upd(a,md);b>>=1ll;
}
return ret;
}
ll pw(ll x,ll k,ll md)
{
ll ret=;
while(k)
{
if(k&1ll)ret=mul(ret,x,md);
x=mul(x,x,md);k>>=1ll;
}
return ret;
}
int phi(ll p,int k)
{
if(!k) return ;
return mul(p-,pw(p,k-,mod),mod);
}
void add(ll a,bool flag)
{
if(flag)
{
for(int i=;i<=tot;i++)
if(p[i]==a)
{ c[i]++;return;}
p[++tot]=a; c[tot]=;
}
else
{
for(int i=;i<=tot;i++)
if(p[i]==a&&c2[i]<c[i])
c2[i]++;//Q=gcd()
}
}
bool ml_rb(ll n)
{
if(n<)return false;if(n==)return true;if((n&)==)return false;
ll u=n-,t=;
while((u&)==){u>>=,t++;}
int s=;
while(s--)
{
ll a=rand()%(n-)+;//2~n-1
a=pw(a,u,n); ll pre=a;
for(int i=;i<=t;i++)
{
a=mul(a,a,n);
if(a==&&pre!=&&pre!=n-)return false;
pre=a;
}
if(a!=) return false;
}
return true;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll pl_rho(ll x,ll c)
{
ll x0=rand()%x,y0=x, k=,i=;
while()
{
x0=(mul(x0,x0,x)+c)%x;
ll g=gcd(abs(x0-y0),x);
if(g!=&&g!=x) return g;
if(x0==y0)return x;
if((++i)==k){k<<=;y0=x0;}
}
}
void fd_fac(ll n,bool flag)
{
if(n<)return;
if(ml_rb(n))
{
add(n,flag);return;
}
ll p=n;
while(p==n)p=pl_rho(p,rand()%(n-)+);
fd_fac(p,flag); fd_fac(n/p,flag);
}
int main()
{
srand(time()); n=rdl();
for(int i=;i<=n;i++)
{
ll a=rdl(); P=mul(P,a,mod);
fd_fac(a,);
}
for(int i=;i<=n;i++)
{
ll a=rdl(); if(!a){flag=;continue;}
if(flag)continue;
fd_fac(a,);
}
if(flag)for(int i=;i<=tot;i++)c2[i]=c[i];
for(int i=,d,tp;i<=tot;i++)
{
tp=pw(p[i],c[i]-,mod);
d=mul(tp,mul(p[i]-,c2[i]+,mod)+(c[i]==c2[i]),mod);
ans=mul(ans,d,mod);
}
printf("%d\n",ans);
return ;
}