[BZOJ4176]Lucas的数论

[BZOJ4176]Lucas的数论

试题描述

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。

在整理以前的试题时，发现了这样一道题目“求 (sum_{i=1}^n f(i))”，其中 (f(i)) 表示 (i) 的约数个数。他现在长大了，题目也变难了。

求如下表达式的值：

[sum_{i=1}^n sum_{j=1}^n f(ij) ]

其中 (f(ij)) 表示 (ij) 的约数个数。

他发现答案有点大，只需要输出模 (1000000007) 的值。

输入

第一行一个整数 (n)。

输出

一行一个整数 (ans)，表示答案模 (1000000007) 的值。

输入示例

2

输出示例

8

数据规模及约定

(T le 10000)

(1 le a,b le 10^7)

题解

首先用到个结论

[f(ij) = sum_{x|i} sum_{y|j} [gcd(x, y) = 1] ]

这个东西用归纳法证明，先考虑 (i) 和 (j) 都只有一种质因子的情况，若该质因子的次数分别为 (a) 和 (b)，那么等号两边显然都是 (a+b+1)；然后考虑给两个数都加进去一个质因子 (p')，它的次数分别为 (a') 和 (b')，那么等号左边会乘上 (a'+b'+1)，等号右边对于所有之前存在的互质数对可以选择在 (x) 上乘上 (p^1, p^2, cdots , p^{a'}) 或者在 (y) 上乘上 (p^1, p^2, cdots, p^{b'})，还有一种情况就是两边都不乘，所以总共是 (a'+b'+1) 种。

那么这题就可以开始大力反演了

[sum_{i=1}^n sum_{j=1}^n f(ij) \ = sum_{i=1}^n sum_{j=1}^n sum_{x|i} sum_{y|j} [gcd(x, y) = 1] \ = sum_{x=1}^n sum_{y=1}^n { [gcd(x, y) = 1] lfloor frac{n}{x} floor lfloor frac{n}{y} floor } \ = sum_{x=1}^n { lfloor frac{n}{x} floor sum_{y=1}^n { lfloor frac{n}{y} floor sum_{d|x, d|y} mu(d) } } \ = sum_{d=1}^n { mu(d) sum_{x=1}^{lfloor frac{n}{d} floor} lfloor frac{n}{xd} floor sum_{y=1}^{lfloor frac{n}{d} floor} lfloor frac{n}{yd} floor } \ = sum_{d=1}^n { mu(d) left( sum_{x=1}^{lfloor frac{n}{d} floor} lfloor frac{n}{xd} floor ight) ^2 } ]

(mu(d)) 的前缀和可以很容易杜教筛了，现在就是要快速求 (sum_{i=1}^n lfloor frac{n}{i} floor)，直接做 (O(sqrt n))，总复杂度是 (O(n^{frac{3}{4}})) 的，考虑用杜教筛优化这个。

[sum_{i=1}^n lfloor frac{n}{i} floor \ = sum_{i=1}^n sum_{j=1}^{lfloor frac{n}{i} floor} 1 \ = sum_{i=1}^n sum_{ij le n} 1 ]

现在我们用 (j) 替换 (ij)

[= sum_{i=1}^n sum_{i|j, j le n} 1 \ = sum_{j=1}^n sum_{i|j} 1 \ = sum_{j=1}^n f(j) ]

(f(j)) 显然可以线性筛，这样就可以预处理前 (n^{frac{2}{3}}) 项了！

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
#include <map>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 1000010
#define MOD 1000000007
#define LL long long

bool vis[maxn];
int prime[maxn], cp, mu[maxn], f[maxn], smu[maxn], sf[maxn], largest[maxn], nos[maxn];
void init() {
	int n = maxn - 1;
	mu[1] = smu[1] = f[1] = sf[1] = 1; largest[1] = 0;
	rep(i, 2, n) {
		if(!vis[i]) prime[++cp] = i, mu[i] = -1, f[i] = 2, largest[i] = 1, nos[i] = 1;
		for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) {
				mu[i*prime[j]] = 0;
				f[i*prime[j]] = f[nos[i]] * (largest[i/nos[i]] + 2);
				largest[i*prime[j]] = max(largest[nos[i]], largest[i/nos[i]] + 1);
				nos[i*prime[j]] = nos[i];
				break;
			}
			mu[i*prime[j]] = -mu[i];
			f[i*prime[j]] = f[i] << 1;
			largest[i*prime[j]] = max(largest[i], 1);
			nos[i*prime[j]] = i;
		}
		smu[i] = smu[i-1] + mu[i];
		if(smu[i] >= MOD) smu[i] -= MOD;
		if(smu[i] < 0) smu[i] += MOD;
		sf[i] = sf[i-1] + f[i];
		if(sf[i] >= MOD) sf[i] -= MOD;
	}
	return ;
}

map <int, int> hSmu, hSf;

int Smu(int n) {
	if(n < maxn) return smu[n];
	if(hSmu.count(n)) return hSmu[n];
	int ans = 1;
	for(int i = 2; i <= n; ) {
		int r = min(n / (n / i), n);
		ans -= (LL)(r - i + 1) * Smu(n / i) % MOD;
		if(ans < 0) ans += MOD;
		i = r + 1;
	}
	return hSmu[n] = ans;
}

int Sf(int n) {
	if(n < maxn) return sf[n];
	if(hSf.count(n)) return hSf[n];
	int ans = 0;
	for(int i = 1; i <= n; ) {
		int r = min(n / (n / i), n);
		ans += (LL)(r - i + 1) * (n / i) % MOD;
		if(ans >= MOD) ans -= MOD;
		i = r + 1;
	}
	return hSf[n] = ans;
}

LL sqr(int x) { return (LL)x * x % MOD; }

int main() {
	init();
	
	int n = read(), ans = 0;
	for(int d = 1; d <= n; ) {
		int r = min(n / (n / d), n);
		ans += (LL)(Smu(r) - Smu(d - 1) + MOD) * sqr(Sf(n / d)) % MOD;
		if(ans >= MOD) ans -= MOD;
		d = r + 1;
	}
	
	printf("%d
", ans);
	
	return 0;
}

稍微改一改你还可以用它 A 掉 BZOJ3994（注意这题有 (n, m)，不再取模，且必须开 long long，哦对还有多组数据）。