BZOJ4652 [Noi2016]循环之美 【数论 + 莫比乌斯反演 + 杜教筛】

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<‘ ‘; puts("");
#define mp(a,b) make_pair<LL,LL>(a,b)
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000;
map<LL,LL> _mu;
map<LL,LL>::iterator it;
map<pair<LL,LL>,LL> G;
map<pair<LL,LL>,LL>::iterator IT;
int N;
int p[maxn],pi,isn[maxn];
LL Smu[maxn],f[maxn],mu[maxn];
LL n,m,K;
int gcd(int a,int b){return b ? gcd(b,a % b) : a;}
void init(){
    N = 1000000;
    mu[1] = 1;
    for (int i = 2; i <= N; i++){
        if (!isn[i]) p[++pi] = i,mu[i] = -1;
        for (int j = 1; j <= pi && i * p[j] <= N; j++){
            isn[i * p[j]] = true;
            if (i % p[j] == 0){
                mu[i * p[j]] = 0;
                break;
            }
            mu[i * p[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= N; i++) Smu[i] = Smu[i - 1] + mu[i];
    for (int i = 1; i <= K; i++) f[i] = f[i - 1] + (gcd(i,K) == 1);
}
LL S(LL n){
    if (n <= N) return Smu[n];
    if ((it = _mu.find(n)) != _mu.end()) return it->second;
    LL ans = 1;
    for (LL i = 2,nxt; i <= n; i = nxt + 1){
        nxt = n / (n / i);
        ans -= (nxt - i + 1) * S(n / i);
    }
    return _mu[n] = ans;
}
LL F(LL n){
    return (n / K) * f[K] + f[n % K];
}
LL g(LL n,LL k){
    if ((IT = G.find(mp(n,k))) != G.end())
        return IT->second;
    if (n == 0) return 0;
    if (k == 1) return S(n);
    LL ans = 0;
    for (LL i = 1; i * i <= k; i++){
        if (k % i == 0){
            if (mu[i]) ans += g(n / i,i);
            if (i * i != k && mu[k / i])
                ans += g(n / (k / i),k / i);
        }
    }
    return G[mp(n,k)] = ans;
}
int main(){
    cin >> n >> m >> K;
    init();
    LL ans = 0,now,last = 0;
    for (LL i = 1,nxt; i <= min(n,m); i = nxt + 1){
        nxt = min(n / (n / i),m / (m / i));
        now = g(nxt,K);
        ans += 1ll * (now - last) * (n / i) * F(m / i);
        last = now;
    }
    cout << ans << endl;
    return 0;
}