原文链接https://www.cnblogs.com/zhouzhendong/p/CF1016G.html
题目传送门 - CF1016G
题意
给定 $n,x,y$ ，以及一个含有 $n$ 个元素的数组 $a$ 。
我们称一个数对 $(i,j)$ 是合法的，当且仅当存在一个 $v$ ，使得 $gcd(a_i,v)=x$ 且 ${ m lcm} (a_j,v)=y$ 。
请你统计有多少合法的 $(i,j)$ 数对。
$nleq 2 imes 10^5,1leq x,y,a_ileq 10^{18}$
题解
先用 Pollard_Rho 将 $y$ 分解一下。由于最小的 $16$ 个不同的质数的乘积大于 $10^{18}$ ，所以， $y$ 最多只有 $15$ 个不同的质因子。
然后，我们罗列以下简单性质：
性质1 ： $x|y$ 。 证明： 因为 $gcd(a_i,v)=x$ ，所以 $x|v$ ； 因为 ${ m lcm}(a_j,v)=y$ ，所以 $v|y$ 。因为 $x|v,v|y$ ，所以 $x|y$ 。
性质2 ： $x|a_i, a_j|y$ 。
性质3 ： 令集合 $P$ 为素数集。定义一个数 $alpha$ 分解质因数后，得到的结果中，含有质因子 $p(pin P)$ 的个数为 $c_{a_i,p}$ 。则因为 $x|y$ ，所以 $forall pin P , c_{x,p}leq c_{y,p}$ 。
性质4 ： 数对 $(i,j)$ 是合法的，对应的数是 $v$ 。那么 $forall pin P, c_{x,p}leq min(c_{a_i,p},c_{v,p}),max(c_{a_j,p},c_{v,p})leq c_{y,p}$ 。
于是，我们就可以得到快速判定两个 $(i,j)$ 是否合法了。
性质5 ： $X={a_i} (x|a_i),Y={a_i} (a_i|y)$ ，如果一个数对 $(i,j)$ 合法，则根据性质2和性质4，必然有： $a_iin X,a_jin Y; forall pin P, c_{a_i,p}=c_{x,p} 或 c_{a_j,p}=c_{y,p}$ 。
由于一开始说到的， $y$ 最多只有 $15$ 个不同的质因子，所以对于 $pin P$ 且 $c_{y,p}=0$ ，由于 $a_j|y$ ，所以 $c_{a_j,p}=c_{y,p}=0$ ，所以这个不需要考虑。
考虑枚举 $a_iin X$ ，对于一个数 $i$ ，对应的合法的 $j$ 有多少？我们考虑写出合法的 $a_j$ 的范围。
对于每一个 $pin P$ ，那么我们需要询问的是一下两种情况下满足条件的数的个数：如果 $c_{a_i,p}=c_{x,p}$ ， $c_{a_j,p}=c_{y,p} { m OR} c_{a_j,p}<c_{y,p}$ ；如果 $c_{a_i,p}>c_{x,p}$ ，那么 $c_{a_j,p}=c_{y,p}$ 。
由于我们最多只需要涉及 $15$ 个质因子，我们可以把集合 $y$ 中的数的信息压缩一下，得到一个二进制表示。令 $y$ 的第 $i$ 个质数为 $p_i$ ，则对于一个数 $a_j$ ，如果 $c_{a_j,p_i}=c_{y,p_i}$ 那么对应二进制位为 $1$ ，否则为 $0$ 。于是，每一次询问就是要求一下对应的数的个数。不难发现，我们只需要统计一下每一个二进制表示对应的在集合 $y$ 中有多少数字，然后做一个 and 卷积的FWT，就可以得到需要询问的东西了。
令 $d=15$ ，所以时间复杂度为 $O ((n+2^d)log d)$ 。
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
typedef long double LD;
namespace Pollard_Rho{
int prime[9]={2,3,5,7,11,13,17,19,23};
ULL RR;
int Pcnt;
LL p[70];
vector <LL> res;
LL R(LL mod){
return (RR+=4179340454199820289LL)%mod;
}
LL Mul(LL x,LL y,LL mod){
LL d=(LL)floor((LD)x*y/mod+0.5);
LL res=x*y-d*mod;
if (res<0)
res+=mod;
return res;
}
LL Pow(LL x,LL y,LL mod){
LL ans=1%mod;
for (;y;y>>=1,x=Mul(x,x,mod))
if (y&1)
ans=Mul(ans,x,mod);
return ans;
}
bool Miller_Rabin(LL n){
if (n<=1)
return 0;
for (int i=0;i<9;i++)
if (n==prime[i])
return 1;
LL d=n-1;
int tmp=0;
while (!(d&1))
d>>=1,tmp++;
for (int i=0;i<9;i++){
LL x=Pow(prime[i],d,n),p=x;
for (int j=1;j<=tmp;j++){
x=Mul(x,x,n);
if (x==1&&p!=1&&p!=n-1)
return 0;
p=x;
}
if (x!=1)
return 0;
}
return 1;
}
LL f(LL x,LL c,LL mod){
return (Mul(x,x,mod)+c)%mod;
}
LL gcd(LL x,LL y){
return y?gcd(y,x%y):x;
}
LL Get_Factor(LL c,LL n){
LL x=R(n),y=f(x,c,n),p=n;
while (x!=y&&(p==n||p==1)){
p=gcd(n,max(x-y,y-x));
x=f(x,c,n);
y=f(f(y,c,n),c,n);
}
return p;
}
void Pollard_Rho(LL n){
if (n<=1)
return;
if (Miller_Rabin(n)){
res.push_back(n);
return;
}
while (1){
LL v=Get_Factor(R(n-1)+1,n);
if (v!=n&&v!=1){
Pollard_Rho(v);
Pollard_Rho(n/v);
return;
}
}
}
void work(LL n){
res.clear();
Pollard_Rho(n);
}
}
LL read(){
LL x=0,f=1;
char ch=getchar();
while (!isdigit(ch)&&ch!='-')
ch=getchar();
if (ch=='-')
f=-1,ch=getchar();
while (isdigit(ch))
x=(x<<1)+(x<<3)+ch-48,ch=getchar();
return x;
}
const int N=200005;
int n;
LL x,y,a[N];
LL v[16],vc=0,xc[16],yc[16],c[16];
vector <LL> Mul_x;
LL A[1<<15];
void Get_V(){
memset(xc,0,sizeof xc);
memset(yc,0,sizeof yc);
Pollard_Rho :: work(y);
vector <LL> xx=Pollard_Rho :: res;
sort(xx.begin(),xx.end());
v[vc=1]=xx[0],yc[vc]=1;
for (int i=1;i<xx.size();i++)
if (xx[i]==xx[i-1])
yc[vc]++;
else
v[++vc]=xx[i],yc[vc]=1;
LL xv=x;
for (int i=1;i<=vc;i++)
while (xv%v[i]==0)
xv/=v[i],xc[i]++;
}
void FWT(LL a[],int n,int flag){
for (int d=1;d<n;d<<=1)
for (int i=0;i<n;i+=(d<<1))
for (int j=0;j<d;j++)
a[i+j]+=a[i+j+d]*flag;
}
int main(){
n=read(),x=read(),y=read();
for (int i=1;i<=n;i++)
a[i]=read();
if (y%x){
putchar('0');
return 0;
}
if (y==1){
int ans=0;
for (int i=1;i<=n;i++)
if (a[i]==1)
ans++;
printf("%I64d",1LL*ans*ans);
return 0;
}
Get_V();
Mul_x.clear();
memset(A,0,sizeof A);
for (int i=1;i<=n;i++){
if (a[i]%x==0)
Mul_x.push_back(a[i]);
if (y%a[i]==0){
memset(c,0,sizeof c);
LL yv=a[i];
for (int i=1;i<=vc;i++)
while (yv%v[i]==0)
yv/=v[i],c[i]++;
int s=0;
for (int i=1;i<=vc;i++)
if (c[i]==yc[i])
s|=1<<(i-1);
A[s]++;
}
}
FWT(A,1<<vc,1);
LL ans=0;
for (int i=0;i<Mul_x.size();i++){
LL now=Mul_x[i];
memset(c,0,sizeof c);
LL yv=now;
for (int i=1;i<=vc;i++)
while (yv%v[i]==0)
yv/=v[i],c[i]++;
int s=(1<<vc)-1;
for (int i=1;i<=vc;i++)
if (c[i]==xc[i]||xc[i]==yc[i])
s^=1<<(i-1);
ans+=A[s];
}
printf("%I64d",ans);
return 0;
}