求:$sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1][(j,k)=1]$
这个时候可以拆前面的,也可以拆后面的.
由于后面的 $k$ 是一个定值,考虑拆解后面的部分.
得:$sum_{d|k} mu(d) sum_{i=1}^{n} sum_{j=1}^{frac{m}{d}} [(i,jd)=1]$
$Rightarrow sum_{d|k} mu(d) sum_{i=1}^{n} sum_{j=1}^{frac{m}{d}} [(i,j)=1][(i,d)=1]$
令 $S(n,m,k)=sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1][(j,k)=1]$
则有 $S(n,m,k)=sum_{d|k} mu(d)S(frac{m}{d},n,d)$
其中边界条件是 $n,m leqslant 0$ 和 $k=1$.
$k=1$ 时求的是 $sum_{i=1}^{n} sum_{j=1}^{m} [(i,j)=1]$
也就是 $sum_{d=1}^{min(n,m)} mu(d)frac{n}{d}frac{m}{d}$
这个需要用到杜教筛来计算 $mu$ 的前缀和.
$Rightarrow S(n)=frac{sum_{i=1}^{n}(f*g)(i)-sum_{i=2}^{n}S(frac{n}{i})g(i)}{g(1)}$
我们选 $g=I(x)$,则 $(f*g)(i)=[i=1]$.
则 $S(n)=1-sum_{i=2}^{n} S(frac{n}{i})$
时间复杂度不会计算,但是跑的挺快.
代码:
#include <map> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #define N 20000008 #define ll long long #define pb push_back #define setIO(s) freopen(s".in","r",stdin) using namespace std; map<int,ll>ansmu; map<int,int>vised; vector<int>G[2003]; int vis[N],mu[N],prime[N],sum[N],tot; void init() { mu[1]=1; for(int i=2;i<N;++i) { if(!vis[i]) prime[++tot]=i,mu[i]=-1; for(int j=1;j<=tot&&i*prime[j]<N;++j) { vis[i*prime[j]]=1; if(i%prime[j]) { mu[i*prime[j]]=-mu[i]; } else { mu[i*prime[j]]=0; break; } } } for(int i=1;i<N;++i) { sum[i]=sum[i-1]+mu[i]; } for(int i=1;i<2003;++i) { for(int j=1;j<=i;++j) if((i%j)==0&&mu[j]) G[i].pb(j); } } ll gcd(ll a,ll b) { return b?gcd(b,a%b):a; } ll getmu(int n) { if(n<N) return sum[n]; if(vised[n]) return ansmu[n]; ll ans=0; for(int i=2,j;i<=n;i=j+1) { j=n/(n/i),ans+=(j-i+1)*getmu(n/i); } vised[n]=1,ansmu[n]=1ll-ans; return 1ll-ans; } ll calc(int n,int m) { if(n>m) swap(n,m); ll ans=0; for(int i=1,j;i<=n;i=j+1) { j=min(n/(n/i),m/(m/i)); ans+=1ll*(n/i)*(m/i)*(getmu(j)-getmu(i-1)); } return ans; } ll solve(int n,int m,int k) { if(n<=0||m<=0) return 0; if(k==1) { return calc(n,m); } ll ans=0; for(int i=0;i<G[k].size();++i) { int cur=G[k][i]; ans+=1ll*mu[cur]*solve(m/cur,n,cur); } return ans; } int main() { // setIO("input"); int n,m,k; scanf("%d%d%d",&n,&m,&k); ll ans=0; init(); printf("%lld ",solve(n,m,k)); return 0; }