bzoj 3481 DZY loves math —— 反演+Pollard_rho分解质因数

题目：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_m ^ t_m，Q = p₁ ^ k₁ * p₂ ^ k₂ * ... * p_m ^ k_m，则原式

= ∏( 1 <= i <= m ) ∑( 0 <= x <= k[i] ) p_i * φ(p_i^t_i-x)

而 φ(p_i^t_i) = p_i^t_i * ( p_i - 1 ) / p_i

所以原式最终变成 ∏( 1 <= i <= m ) p_i^(t_i-1) * ( p_i - 1 ) * ( k_i + 1 )

注意 t[i] = k[i] 时略有不同，求和的那一项 = p_i ， 所以如果有 t[i] = k[i]，则再加一个 p_i^(t_i-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=1005;
int m,S=20,t[xn],k[xn];
ll pri[xn],ans;
ll rd()
{
  ll ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=(ret<<3ll)+(ret<<1ll)+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll Rand(ll n)
{
    return (((ll)rand()<<31)|rand())%(n-1)+1;
//    return rand()%(n-1)+1;//也可，但无ll 
}
ll upt(ll x,ll mod){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
ll mul(ll a,ll b,ll mod)
{
  ll ret=0;  a=a%mod; b=b%mod;
  for(;b;b>>=1ll,a=upt(a+a,mod))
    if(b&1)ret=upt(ret+a,mod);
  return ret;
}
ll pw(ll a,ll b,ll mod)
{
  ll ret=1;  a=a%mod;
  for(;b;b>>=1ll,a=mul(a,a,mod))
    if(b&1)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=1;i<=t;i++)
    {
      x=mul(x,x,p);
      if(x==1&&lst!=1&&lst!=p-1)return 0; // a^2 % p = 1 -> a=1,a=p-1
      lst=x;
    }
  if(x!=1)return 0; // a^p-1 % p = 1
  return 1;
}
bool MR(ll x)
{
  if(x<=1||x&1==0)return 0;
  if(x==2)return 1;
  ll t=0,u=x-1;
  while(u&1==0)u>>=1ll,t++;
  for(int i=1;i<=S;i++)
    {
      ll a=Rand(x)%(x-1)+1;
      if(!ck(a,u,t,x))return 0;
    }
  return 1;
}
ll rho(ll n,ll c)
{
  ll x=Rand(n)%n,y=x;
  ll i=1,k=2;
  while(1)
    {
      i++;
      x=(mul(x,x,n)+c)%n;
      ll d=gcd((x-y+n)%n,n);
      if(d!=1&&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=1;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==1)return;//!
  if(MR(x)){add(x,tp); return;}
  ll p=x;
  while(p>=x)p=rho(p,Rand(x)%(x-1)+1);
  find(p,tp); find(x/p,tp);
}
int main()
{
  srand(23333);//或 time(0); 
  int mod=1e9+7;
  int n=rd(); ll x; bool fl=0;
  for(int i=1;i<=n;i++)x=rd(),find(x,1);//1:P-t[]
  for(int i=1;i<=n;i++)
  {
      x=rd();
      if(!x){fl=1; continue;}
    else find(x,0);
  }
  if(fl)for(int i=1;i<=m;i++)k[i]=t[i];
  else for(int i=1;i<=m;i++)k[i]=min(k[i],t[i]);//(P,Q)
  ans=1;
  for(int i=1;i<=m;i++)
    {
      ll pp=pw(pri[i],t[i]-1,mod);
      ll mm=mul(pri[i]-1,k[i]+1,mod);
      if(k[i]==t[i])mm++;
      ans=mul(ans,mul(pp,mm,mod),mod);
    }
  printf("%lld\n",ans);
  return 0;
}