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=650,mod=1e9+7;
int n,ans=1,c[N],c2[N],tot;
ll p[N],P=1;
bool vis[N],flag;
ll rdl()
{
  ll ret=0;bool fx=1;char ch=getchar();
  while(ch>'9'||ch<'0'){if(ch=='-')fx=0;ch=getchar();}
  while(ch>='0'&&ch<='9') ret=(ret<<3ll)+(ret<<1ll)+ch-'0',ch=getchar();
  return fx?ret:-ret;
}
void upd(ll &x,ll md){x-=(x>=md?md:0);}
ll mul(ll a,ll b,ll md)
{
  ll ret=0;//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=1;
  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 1;
  return mul(p-1,pw(p,k-1,mod),mod);
}
void add(ll a,bool flag)
{
  if(flag)
    {
      for(int i=1;i<=tot;i++)
    if(p[i]==a)
      { c[i]++;return;}
      p[++tot]=a; c[tot]=1;
    }
  else
    {
      for(int i=1;i<=tot;i++)
    if(p[i]==a&&c2[i]<c[i])
      c2[i]++;//Q=gcd()
    }
}
bool ml_rb(ll n)
{
  if(n<2)return false;if(n==2)return true;if((n&1)==0)return false;
  ll u=n-1,t=0;
  while((u&1)==0){u>>=1,t++;}
  int s=20;
  while(s--)
    {
      ll a=rand()%(n-2)+2;//2~n-1
      a=pw(a,u,n); ll pre=a;
      for(int i=1;i<=t;i++)
    {
      a=mul(a,a,n);
      if(a==1&&pre!=1&&pre!=n-1)return false;
      pre=a;
    }
      if(a!=1) 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=1,i=0;
  while(1)
    {
      x0=(mul(x0,x0,x)+c)%x;
      ll g=gcd(abs(x0-y0),x);
      if(g!=1&&g!=x) return g;
      if(x0==y0)return x;
      if((++i)==k){k<<=1;y0=x0;}
    }
}
void fd_fac(ll n,bool flag)
{
  if(n<2)return;
  if(ml_rb(n))
    {
      add(n,flag);return;
    }
  ll p=n;
  while(p==n)p=pl_rho(p,rand()%(n-1)+1);
  fd_fac(p,flag); fd_fac(n/p,flag);
}
int main()
{
  srand(time(0));  n=rdl();
  for(int i=1;i<=n;i++)
    {
      ll a=rdl(); P=mul(P,a,mod);
      fd_fac(a,1);
    }
  for(int i=1;i<=n;i++)
    {
      ll a=rdl(); if(!a){flag=1;continue;}
      if(flag)continue;
      fd_fac(a,0);
    }
  if(flag)for(int i=1;i<=tot;i++)c2[i]=c[i];
  for(int i=1,d,tp;i<=tot;i++)
    {
      tp=pw(p[i],c[i]-1,mod);
      d=mul(tp,mul(p[i]-1,c2[i]+1,mod)+(c[i]==c2[i]),mod);
      ans=mul(ans,d,mod);
    }
  printf("%d
",ans);
  return 0;
}