NOI 2016 循环之美

#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\n",G(tot,last),n/i,F(m/i));
        ans+=(G(tot,last)-G(tot,i-1))*(n/i)*F(m/i);

    }
    printf("%lld\n",ans);
    return 0;
}