【LOJ】#2085. 「NOI2016」循环之美

题解

我们要求的其实是这个东西= =
\(\sum_{i = 1}^{n}\sum_{j = 1}^{n}[(i,j) == 1][(j,k) == 1]\)
然后变一下形
\(\sum_{j = 1}^{n}[(j,k) == 1]\sum_{i = 1}^{n}[(i,j) == 1]\)
\(\sum_{j = 1}^{n}[(j,k) == 1]\sum_{i = 1}^{n}\sum_{d|i,j}\mu(d)\)
\(\sum_{j = 1}^{n}[(j,k) == 1]\sum_{d | j} \mu(d) \lfloor \frac{n}{d} \rfloor\)
\(\sum_{d = 1}^{n} \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{y = 1}^{\lfloor \frac{m}{d} \rfloor}[(yd,k) == 1]\)
\(\sum_{d = 1}^{n} \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{y = 1}^{\lfloor \frac{m}{d} \rfloor}[(d,k) == 1][(y,k) == 1]\)
\(\sum_{d = 1}^{n} [(d,k) == 1] \mu(d) \lfloor \frac{n}{d} \rfloor \sum_{y = 1}^{\lfloor \frac{m}{d} \rfloor}[(y,k) == 1]\)

设\(f(n) = \sum_{i = 1}^{n}[(i,k) == 1]\)
我们可以预处理\(f(1)\)到\(f(k)\)
那么就有
\(f(n) = \lfloor\frac{n}{k}\rfloor f(k) + f(n mod k)\)
因为\((a,b) == 1\) 等价于\((a mod b,b) == 1\)

我们现在可以\(O(n)\)解决这个式子了，但是还不够

我们可以用数论分块处理\(\lfloor \frac{n}{d} \rfloor\) 和 \(\lfloor \frac{m}{d} \rfloor\)
我们尝试算
\(\sum_{d = 1}^{n} [(d,k) == 1] \mu(d)\)
设\(g(n,k) = \sum_{d = 1}^{n} [(d,k) == 1] \mu(d)\)
我们选择一个\(k\)的质因子\(p\)，然后把\(k\)表示成\(p^{c}q\)的形式
我们从和\(q\)互质的数里除掉和\(p\)互质的数

\(g(n,k) = \sum_{d = 1}^{n} [(d,q) == 1] \mu(d) - \sum_{y = 1}^{n} [(yp,q) == 1] \mu(yp)\)
由于\(p\)和\(q\)互质，所以我们只需要保证\([(y,p) == 1][(y,q) == 1]\)
\(g(n,k) = \sum_{d = 1}^{n} [(d,q) == 1] \mu(d) - \mu(p) \sum_{y = 1}^{n} [(y,p) == 1][(y,q) == 1] \mu(y)\)
\(g(n,k) = \sum_{d = 1}^{n} [(d,q) == 1] \mu(d) + \sum_{y = 1}^{\lfloor \frac{n}{d} \rfloor} [(y,k) == 1]\mu(y)\)
\(g(n,k) = g(n,q) + g(\lfloor \frac{n}{p} \rfloor,k)\)

边界是\(n <= 1\)返回\(n\)，\(k = 1\)返回一个莫比乌斯函数的前缀和，可以杜教筛

代码

#include <bits/stdc++.h>
//#define ivorysi
#define enter putchar('\n')
#define space putchar(' ')
#define fi first
#define se second
#define pb push_back
#define mp make_pair
#define eps 1e-8
#define mo 974711
#define MAXN 1000000
#define pii pair<int,int>
using namespace std;
typedef long long int64;
typedef double db;
template<class T>
void read(T &res) {
    res = 0;char c = getchar();T f = 1;
    while(c < '0' || c > '9') {
	if(c == '-') f = -1;
	c = getchar();
    }
    while(c >= '0' && c <= '9') {
	res = res * 10 + c - '0';
	c = getchar();
    }
    res *= f;
}
template<class T>
void out(T x) {
    if(x < 0) {putchar('-');x = -x;}
    if(x >= 10) {
	out(x / 10);
    }
    putchar('0' + x % 10);
}
int N,M,K;
int mu[MAXN + 5],prime[MAXN + 5],tot,Mu[MAXN + 5];
bool nonprime[MAXN + 5];
int f[2005];
struct HASH {
    struct node {
	int64 x,v;
	int next;
    }E[MAXN * 2];
    int head[mo + 5],sumE;
    HASH() {
	memset(head,0,sizeof(head));sumE = 0;
    }
    void add(int u,int64 x,int64 v) {
	E[++sumE].x = x;E[sumE].v = v;E[sumE].next = head[u];
	head[u] = sumE;
    }
    void Insert(int64 x,int64 v) {
	add(x % mo,x,v);
    }
    int64 Query(int64 x){
	for(int i = head[x % mo] ; i ; i = E[i].next) {
	    if(E[i].x == x) return E[i].v;
	}
	return -1;
    }
}H[2];
int gcd(int a,int b) {
    return b == 0 ? a : gcd(b,a % b);
}
int64 calcF(int x) {
    return 1LL * (x / K) * f[K] + f[x % K];
}
int64 S(int x) {
    if(x <= MAXN) return Mu[x];
    int64 c = H[0].Query(x);
    if(c != -1) return c;
    int64 res = 0;
    for(int i = 2 ; i <= x ; ++i) {
	int r = x / (x / i);
	res = res + 1LL * (r - i + 1) * S(x / i);
	i = r;
    }
    res = 1 - res;
    H[0].Insert(x,res);
    return res;
}
int64 G(int n,int k) {
    if(k == 1) return S(n);
    else if(n <= 1) return n;
    int64 c = H[1].Query(1LL * (n - 1) * K + k);
    if(c != -1) return c;
    for(int i = 1 ; i <= tot ; ++i) {
	if(k % prime[i] == 0) {
	    int x = k;
	    while(x % prime[i] == 0) x /= prime[i];
	    int64 res = G(n,x) + G(n / prime[i],x * prime[i]);
	    H[1].Insert(1LL * (n - 1) * K + k,res);
	    return res;
	}
    }
    
}
void Solve() {
    int t = min(N,M);
    int64 res = 0;
    for(int i = 1 ; i <= t ; ++i) {
	int r = min(N / (N / i),M / (M / i));
	int64 s = calcF(M / i) * (N / i);
	res = res + s * (G(r,K) - G(i - 1,K));
	i = r;
    }
    out(res);enter;
}
int main() {
#ifdef ivorysi
    freopen("f1.in","r",stdin);
#endif
    read(N);read(M);read(K);
    mu[1] = 1;Mu[1] = 1;
    for(int i = 2 ; i <= MAXN ; ++i) {
	if(!nonprime[i]) {
	    prime[++tot] = i;
	    mu[i] = -1;
	}
	for(int j = 1 ; j <= tot ; ++j) {
	    if(prime[j] > MAXN / i) break;
	    nonprime[i * prime[j]] = 1;
	    if(i % prime[j] == 0) break;
	    else mu[i * prime[j]] = -mu[i];
	}
	Mu[i] = Mu[i - 1] + mu[i];
    }
    for(int i = 1 ; i <= K ; ++i) {
	f[i] = f[i - 1] + (gcd(i,K) == 1);
    }
    Solve();
    return 0;
}