Lucas的数论

Description

• 给定整数 \(n\)，求

  \[\left[\sum_{i = 1}^n \sum_{j = 1}^n d(ij)\right] \bmod (10^9 + 7) \]

• 对于 \(100\%\) 的数据 \(n\le 10^9\)。

Solution

\[\begin{aligned} \sum_{i = 1}^n \sum_{j = 1}^n d(ij) & = \sum_{i = 1}^n \sum_{j = 1}^n \sum_{x\mid i} \sum_{y\mid j} [\gcd(x, y) = 1] \\ & = \sum_{x = 1}^n \sum_{y = 1}^n [\gcd(x, y) = 1] \sum_{i = 1}^n [x\mid i] \sum_{j = 1}^n [y\mid j] \\ & = \sum_{x = 1}^n \sum_{y = 1}^n [\gcd(x, y) = 1] \left\lfloor\dfrac{n}{x}\right\rfloor \left\lfloor\dfrac{n}{y}\right\rfloor \\ & = \sum_{d = 1}^n \mu(d) \sum_{x = 1}^n [d\mid x] \left\lfloor\dfrac{n}{x}\right\rfloor \sum_{y = 1}^n [d\mid y] \left\lfloor\dfrac{n}{y}\right\rfloor \\ & = \sum_{d = 1}^n \mu(d) \left(\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} \left\lfloor\dfrac{n}{di}\right\rfloor\right)^2 & (1)\\ & = \sum_{d = 1}^n \mu(d) \left(\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor} d(i)\right)^2 & (2) \end{aligned} \]

\((1)\) 中有两层整除分块。

\(\mu\) 用杜教筛，阈值为 \(10^6\)。

预处理 \(d\) 的前缀和。

以 \(k\) 为阈值，当 \(\left\lfloor\dfrac{n}{d}\right\rfloor \le k\) 时用 \((2)\) 的公式，以预处理过的 \(d\) 的前缀和来计算；当 \(\left\lfloor\dfrac{n}{d}\right\rfloor > k\) 时用 \((1)\) 中的整除分块计算。

那么预处理是 \(\Omicron(k)\) 的。

第一层整除分块中，若 \(\left\lfloor\dfrac{n}{d}\right\rfloor \le k\)，即 \(d \ge \left\lfloor\dfrac{n}{k}\right\rfloor\)，那么就会用 \(d\) 的前缀和 \(\Omicron(1)\) 计算，这一段是 \(\Omicron\left(\sqrt{n - \dfrac{n}{k}}\right) \le \Omicron(\sqrt{n})\) 的。

第一层整除分块中，若 \(\left\lfloor\dfrac{n}{d}\right\rfloor > k\)，即 \(d < \left\lfloor\dfrac{n}{k}\right\rfloor\)，那么就会用整除分块来计算，一共是 \(\sum_{d = 1}^{\left\lfloor\frac{n}{k}\right\rfloor} \sqrt{\left\lfloor\dfrac{n}{d}\right\rfloor}\)，即

\[\int_{0}^{\frac{n}{k}} \sqrt{\dfrac{n}{x}}\, dx = \dfrac{2n}{\sqrt{k}} \]

也就是 \(\Omicron\left(\dfrac{n}{\sqrt{k}}\right)\) 的。

综上，总时间复杂度为 \(\Omicron\left(k + \sqrt{n} + \dfrac{n}{\sqrt{k}} \right)\)，\(k\) 大约取 \(n^{\frac{2}{3}}\) 时最优，和杜教筛一样。

Code

// 18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#include <unordered_map>x
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

const int MAXN = 1e6 + 5;
const int N = 1e6;
const int MOD = 1e9 + 7;

int p[MAXN], mu[MAXN], sum_mu[MAXN], d[MAXN], sum_d[MAXN], num[MAXN];
bool vis[MAXN];

void pre()
{
	mu[1] = sum_mu[1] = d[1] = sum_d[1] = 1;
	for (int i = 2; i <= N; i++)
	{
		if (!vis[i])
		{
			p[++p[0]] = i;
			mu[i] = -1;
			d[i] = 2;
			num[i] = 1;
		}
		for (int j = 1; j <= p[0] && i * p[j] <= N; j++)
		{
			vis[i * p[j]] = true;
			if (i % p[j] == 0)
			{
				mu[i * p[j]] = 0;
				d[i * p[j]] = d[i] / (num[i] + 1) * (num[i] + 2);
				num[i * p[j]] = num[i] + 1;
				break;
			}
			mu[i * p[j]] = mu[i] * mu[p[j]];
			d[i * p[j]] = d[i] * d[p[j]];
			num[i * p[j]] = 1;
		}
		sum_mu[i] = sum_mu[i - 1] + mu[i];
		sum_d[i] = (sum_d[i - 1] + d[i]) % MOD;
	}
}

unordered_map<int, int> dp_mu;

int sublinear_mu(int n)
{
	if (n <= N)
	{
		return sum_mu[n];
	}
	if (dp_mu.find(n) != dp_mu.end())
	{
		return dp_mu[n];
	}
	int res = 1;
	for (int l = 2, r; l <= n; l = r + 1)
	{
		int k = n / l;
		r = n / k;
		res = (res - (ll)(r - l + 1) * sublinear_mu(k) % MOD + MOD) % MOD;
	}
	return dp_mu[n] = res;
}

int getsum_mu(int l, int r)
{
	return (sublinear_mu(r) - sublinear_mu(l - 1) + MOD) % MOD;
}

int getsum_d(int n)
{
	if (n <= N)
	{
		return sum_d[n];
	}
	int res = 0;
	for (int l = 1, r; l <= n; l = r + 1)
	{
		int k = n / l;
		r = n / k;
		res = (res + (ll)(r - l + 1) * k % MOD) % MOD;
	}
	return res;
}

int block(int n)
{
	int res = 0;
	for (int l = 1, r; l <= n; l = r + 1)
	{
		int k = n / l;
		r = n / k;
		int tmp = getsum_d(k);
		res = (res + (ll)getsum_mu(l, r) * tmp % MOD * tmp % MOD) % MOD;
	}
	return res;
}

int main()
{
	pre();
	int n;
	scanf("%d", &n);
	printf("%d\n", block(n));
	return 0;
}