众所周知,反演题只有套路题和神仙题
概括题意:
求$\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]$
蛤?你问我为啥是这个式子?
道可道,非常道。只可意会,不可言传。请出门右转题解区。(一言以蔽之,我不会推)
大跃进反演:
$\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]$
$=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{a=1}^n\mu(a)[a|i][a|j]\sum\limits_{b=1}^m\mu(b)[b|j][b|k]$
$=\sum\limits_{a=1}^n\mu(a)\left\lfloor\dfrac{n}{a}\right\rfloor\sum\limits_{b|k}\mu(b)\left\lfloor\dfrac{m}{lcm(a,b)}\right\rfloor$
这什么沙雕?!
事实证明,大跃进不可取
五五开反演:
设$f(n,m,k)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=1][gcd(j,k)=1]$
$f(n,m,k)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=1]\sum\limits_{d|k}\mu(d)[d|j]$
$=\sum\limits_{d|k}\mu(d)\sum\limits_{i=1}^n\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}[gcd(i,jd)=1]$
$=\sum\limits_{d|k}\mu(d)\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum\limits_{i=1}^n[gcd(i,j)=1][gcd(i,d)=1]$
$=\sum\limits_{d|k}\mu(d)f(\left\lfloor\dfrac{m}{d}\right\rfloor,n,d)$
暴力递归算就好了。
边界:$f(n,0,k)=f(0,m,k)=0$,$f(n,m,1)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)=1]=\sum\limits_{d=1}^n\mu(d)\left\lfloor\dfrac{n}{d}\right\rfloor\left\lfloor\dfrac{m}{d}\right\rfloor$,还要套个杜教筛。
蛤?你问我复杂度?
当然是O(能过)查阅一下题解是$O(\log n\log m\sqrt{n}+n^{\frac{2}{3}})$的。(我TM真不会推啊)
代码:
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <vector>
#define maxn 10000005
#define inf 0x3f3f3f3f
using namespace std;
inline int read(){
int x=0,y=0;
char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')y=1;ch=getchar();}
while(ch>='0'&&ch<='9')x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
return y?-x:x;
}
int mu[maxn],prime[maxn>>2],cnt,sq;
bool is[maxn];
vector<int>fac[2005];
const int mod =3e6 + 7;
struct hash{
int cmp[mod],dat[mod],h[mod],nex[mod],cnt;
void insert(int x,int d){nex[++cnt]=h[x%mod],h[x%mod]=cnt,cmp[cnt]=x,dat[cnt]=d;}
int find(int x){
for(register int i=h[x%mod];i;i=nex[i])
if(cmp[i]==x)return dat[i];
return -inf;
}
}mmu;
int smu(int n){
if(n<maxn)return mu[n];
int ans=mmu.find(n);
if(ans!=-inf)return ans;
ans=1;
for(register int l=2,r;l<=n;l=r+1){
r=n/(n/l);
ans-=(r-l+1)*smu(n/l);
}
mmu.insert(n,ans);
return ans;
}
long long g(int n,int m){
if(n>m)swap(n,m);
long long ans=0;
for(register int l=1,r;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
ans+=1ll*(smu(r)-smu(l-1))*(n/l)*(m/l);
}
return ans;
}
long long f(int n,int m,int k){
if(!n||!m)return 0;
if(k==1)return g(n,m);
long long ans=0;
for(register int i=0,x;i<fac[k].size();++i){
x=fac[k][i];
ans+=f(m/x,n,x)*(mu[x]-mu[x-1]);
}
return ans;
}
int main(){
mu[1]=1;
for(register int i=2;i<maxn;++i){
if(!is[i])prime[++cnt]=i,mu[i]=-1;
for(register int j=1;j<=cnt&&i*prime[j]<maxn;++j){
is[i*prime[j]]=1;
if(i%prime[j]==0){mu[i*prime[j]]=0;break;}
mu[i*prime[j]]=-mu[i];
}
}
for(register int i=1;i<=2000;++i){
if(!mu[i])continue;
for(register int j=i;j<=2000;j+=i)
fac[j].push_back(i);
}
for(register int i=2;i<maxn;++i)mu[i]+=mu[i-1];
int n=read(),m=read(),k=read(),ans=0;
printf("%lld\n",f(n,m,k));
}