题解链接 (不会markdown只能甩链接了)
#include<bits/stdc++.h> #define LL long long #define SIZ 4011007 #define deg using namespace std; int f[3007],smu[SIZ],k,n,m,mu[SIZ],pim[SIZ],tig,s[71],tot,l; LL ans; bool u[SIZ]; map<pair<int,int>,LL> Ma1; map<int,int> Ma; LL F(int x) {return 1ll*x/k*f[k]+f[x%k];} int gcd(int x,int y){return y?gcd(y,x%y):x;} void pre() { for (int i=1;i<=k;i++) f[i]=f[i-1]+(gcd(i,k)==1); smu[1]=1; for (int i=2;i<SIZ;i++) { if (!u[i]) pim[++tig]=i,mu[i]=-1; for (int j=1;j<=tig&&pim[j]*i<SIZ;j++) { u[i*pim[j]]=1; if (i%pim[j]) mu[i*pim[j]]=-mu[i];else break; } smu[i]=smu[i-1]+mu[i]; } } int Su(int x){ if (x<SIZ) return smu[x]; if (Ma.count(x)) return Ma[x]; int ans=1; for (int i=2,last;i<=x;i=last+1) { last=x/(x/i); ans-=(last-i+1)*Su(x/i); } Ma[x]=ans; return ans; } LL G(int x,int y){ if (!x) return Su(y); if (y<=1) return y; pair<int,int> T=make_pair(x,y); int ans=0; if (Ma1.count(T)) return Ma1[T]; Ma1[T]=ans=G(x-1,y)+G(x,y/s[x]); return ans; } int main () { freopen("a.in","r",stdin); scanf("%d%d%d",&n,&m,&k); pre();int KK=k; for (int i=1;pim[i]<=KK;i++) { if (KK%pim[i]==0) s[++tot]=pim[i]; while (KK%pim[i]==0) KK/=pim[i]; } l=min(n,m); for (int i=1,last;i<=l;i=last+1) { last=min(n/(n/i),m/(m/i)); deg("%lld %d %lld ",G(tot,last),n/i,F(m/i)); ans+=(G(tot,last)-G(tot,i-1))*(n/i)*F(m/i); } printf("%lld ",ans); return 0; }