【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} floor)
(sum_{d = 1}^{n} mu(d) lfloor frac{n}{d} floor sum_{y = 1}^{lfloor frac{m}{d} floor}[(yd,k) == 1])
(sum_{d = 1}^{n} mu(d) lfloor frac{n}{d} floor sum_{y = 1}^{lfloor frac{m}{d} floor}[(d,k) == 1][(y,k) == 1])
(sum_{d = 1}^{n} [(d,k) == 1] mu(d) lfloor frac{n}{d} floor sum_{y = 1}^{lfloor frac{m}{d} floor}[(y,k) == 1])

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

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

我们可以用数论分块处理(lfloor frac{n}{d} floor) 和 (lfloor frac{m}{d} floor)
我们尝试算
(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} floor} [(y,k) == 1]mu(y))
(g(n,k) = g(n,q) + g(lfloor frac{n}{p} floor,k))

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

代码

#include <bits/stdc++.h>
//#define ivorysi
#define enter putchar('
')
#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;
}