首先题目是在求:
ans=∑i=1n∑j=1m[i⊥j][j⊥k]ans=∑i=1n∑j=1m[i⊥j][j⊥k]
注意到j和k互质来得简单,因为k只有2000并且是常数.
ans=∑j=1m[j⊥k]∑i=1n[i⊥j]ans=∑j=1m[j⊥k]∑i=1n[i⊥j]
ans=∑j=1m[j⊥k]∑i=1n∑d|i,d|jμ(d)ans=∑j=1m[j⊥k]∑i=1n∑d|i,d|jμ(d)
ans=∑d=1min(n,m)μ(d)⌊nd⌋∑d|j,j≤m[j⊥k]ans=∑d=1min(n,m)μ(d)⌊nd⌋∑d|j,j≤m[j⊥k]
注意到若j=tdj=td那么[j⊥k]=[td⊥k]=[d⊥k][t⊥k][j⊥k]=[td⊥k]=[d⊥k][t⊥k]
ans=∑d=1nμ(d)[d⊥k]⌊nd⌋∑j=1⌊md⌋[j⊥k]ans=∑d=1nμ(d)[d⊥k]⌊nd⌋∑j=1⌊md⌋[j⊥k]
=∑d=1nμ(d)[d⊥k]⌊nd⌋f(⌊md⌋)=∑d=1nμ(d)[d⊥k]⌊nd⌋f(⌊md⌋)
其中f(n)=∑ni=1[i⊥k]=⌊nk⌋ϕ(k)+f(n mod k)f(n)=∑i=1n[i⊥k]=⌊nk⌋ϕ(k)+f(n mod k)
因此f(n)f(n)可以在O(klgk)−O(1)O(klgk)−O(1)内解决
因此对最后的式子做数论分块,那么我们需要求:
g(n,k)=∑nd=1μ(d)[d⊥k]g(n,k)=∑d=1nμ(d)[d⊥k]
这个怎么做呢?设n=pcqn=pcq,其中p是n的某个质因子,并且p⊥qp⊥q
那么考虑g(n,q)g(n,q)中的所有数字,显然除了答案的部分还有和q互质但是不和p互质的部分要减去;
记这些中的数字为y=tp,t≤⌊np⌋,t⊥qy=tp,t≤⌊np⌋,t⊥q(显然因为p是质数所以不和p互质意味着是p的倍数)
考虑我们求的是这些yy的和,那么如果t不和p互质的话μ(y)=0μ(y)=0,因此可以认为t也和p互质,因此t就和k互质。
g(n,k)=g(n,q)−∑t=1⌊np⌋[t⊥k]μ(tp)g(n,k)=g(n,q)−∑t=1⌊np⌋[t⊥k]μ(tp)
注意到莫比乌斯函数是积性函数,而t⊥pt⊥p,因此μ(tp)=μ(t)μ(p)μ(tp)=μ(t)μ(p)
g(n,k)=g(n,q)−μ(p)∑t=1⌊np⌋[t⊥k]μ(t)g(n,k)=g(n,q)−μ(p)∑t=1⌊np⌋[t⊥k]μ(t)
注意到p是质数,所以μ(p)=−1μ(p)=−1
g(n,k)=g(n,q)+∑t=1⌊np⌋[t⊥k]μ(t)=g(n,q)+g(⌊np⌋,k)g(n,k)=g(n,q)+∑t=1⌊np⌋[t⊥k]μ(t)=g(n,q)+g(⌊np⌋,k)
按照固定的顺序除以p,那么有效的状态只有d(k)2d(k)2,其中d(k)d(k)表示k的因子个数,实现可以记忆化,因此只传k此时的质因子个数,会好些一些。
边界条件,当n=0n=0的时候g(n,k)=0g(n,k)=0,当k=1的时候需要杜教筛。综上总复杂度O(d(k)2+n23)O(d(k)2+n23).
#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<map>
#define lint long long
#define N 6000000
#define hv(n,c) (n*10000ll+c)
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
int ms[N],m[N],pri[N];bool np[N];
map<int,int> sav;
map<lint,int> gs;
int v[50],f[2010];
inline int prelude(int n)
{
ms[1]=m[1]=1;
for(int i=2,c=0;i<=n;i++)
{
if(!np[i]) pri[++c]=i,m[i]=-1;ms[i]=ms[i-1]+m[i];
for(int j=1;j<=c&&(lint)i*pri[j]<=n;j++)
{
int x=i*pri[j];np[x]=true;
if(i%pri[j]) m[x]=-m[i];
else { m[x]=0;break; }
}
}
return 0;
}
inline int S(int n)
{
if(n<N) return ms[n];int ans=0ll;
if(sav.count(n)) return sav[n];
for(int s=2,t;s<=n;s=++t)
t=n/(n/s),ans+=(t-s+1)*S(n/t);
return sav[n]=1-ans;
}
inline int g(int n,int c)
{
lint nc=hv(n,c);if(gs.count(nc)) return gs[nc];
if(!n) return 0;if(!c) return gs[nc]=S(n);
return gs[nc]=g(n,c-1)+g(n/v[c],c);
}
inline int g(int a,int b,int c) { return g(b,c)-g(a-1,c); }
inline int F(int n,int k) { return n<=k?f[n]:(n/k)*f[k]+f[n%k]; }
inline int gcd(int a,int b) { return a?gcd(b%a,a):b; }
int main()
{
int n,m,k,c=0;scanf("%d%d%d",&n,&m,&k);
prelude(N-1);lint ans=0ll;int ks=k;
for(int i=2;i<=ks;i++) if(ks%i==0)
{ v[++c]=i;while(ks%i==0) ks/=i; }
for(int i=1;i<=k;i++) f[i]=f[i-1]+(gcd(i,k)==1);
for(int s=1,t;s<=min(n,m);s=++t)
t=min(n/(n/s),m/(m/s)),ans+=(lint)F(m/s,k)*(n/t)*g(s,t,c);
return !printf("%lld\n",ans);
}