代码源[R11F] 题解

前言

一道有点难的数学题，主要考欧拉。

（欧拉筛和欧拉公式）

题面

给定\(N (1 \leq N \leq 5 \times 10^7)\),求

\[(\sum_{i=1}^{n} \sum_{j=1}^{n} \gcd(i,j) ) \mod 998244353 \]

思路

做这道题，我们要先会欧拉函数。

欧拉函数\(\phi\),读作fi，\(\phi(n)\)表示\(1 \sim n\) 中与\(n\)互质的数的个数。
对于所有的质数\(x\)，\(\phi(x) = x-1\)
对于所有数字的求法：
若$n = {p_1}^{r_1} \times {p_2}^{r_2} \times ... \times {p_k}^{r_k} $(分解质因数)，那么 \(\phi(n) = n(1 - \frac{1}{p_1})(1 - \frac{1}{p_2})...(1 - \frac{1}{p_k})\)

接下来，我们要会用欧拉筛\(O(n)\)求\(\phi(1), \phi(2),..., \phi(n)\)

这部分看代码(\(f_i\)表示\(\phi(i)\)）：

const int N = 5e7 + 5;
int p[N], f[N];
long long s[N];
bool b[N];
int main()
{
    for (int i = 2; i <= n; i++)
    {
    	if (!b[i]) // 判断质数（欧拉筛）
    	{
    		p[++c] = i; // 放入质数表（欧拉筛）
    		f[i] = i - 1; // 质数的欧拉函数
    	}
    	for (int j = 1; j <= c && i * p[j] <= n; j++) // 遍历质数（欧拉筛）
    	{
    	    b[i * p[j]] = true; // 筛掉
    	    if (i % p[j] == 0) f[i * p[j]] = f[i] * p[j]; // 如果这个质数已经在i中存在，那么从f[i]变为f[i * p[j]]只需要把n(1-1/p1)(1-1/p2)...(1-1/pk)中的第一个乘数n从i改成i*p[j]，后面已经有了，就不用再乘了。f[i * p[j]] = f[i] * p[j]
    	    else f[i * p[j]] = f[i] * (p[j] - 1); // 如果质数不在i中，那么在把f[i]变为f[i * p[j]]时，不仅把n从i改成i*p[j]，还要再把(1-1/p[j])也乘上去。f[i * p[j]] = f[i] * p[j] * (1 - 1/p[j]) = f[i] * (p[j] - 1)
    	    if (i % p[j] == 0) break;  // 筛过了，停止。
    	}
    }
    return 0;
}

再接下来，我们要对于所有\(n\)，求出每个\(\displaystyle{\sum_{i=1}^{n} \phi(n)}\)

直接对上面的\(f\)进行前缀和即可

for (int i = 1; i <= n; i++)
{
  	s[i] = s[i - 1] + f[i]; //  对f求出前缀和s
}

预处理：

此时，我们便有了所有的预处理代码

const int N = 5e7 + 5;
int p[N], f[N];
long long s[N];
bool b[N];
int main()
{
    int n;
    cin >> n;
    f[1] = 1, f[2] = 1;
    int c = 0;
    for (int i = 2; i <= n; i++) // 求出1~n每个数字的φ
    {
    	if (!b[i])
    	{
    		p[++c] = i;
    		f[i] = i - 1;
    	}
    	for (int j = 1; j <= c && i * p[j] <= n; j++)
    	{
    	    b[i * p[j]] = true;
    	    if (i % p[j] == 0) f[i * p[j]] = f[i] * p[j];
    	    else f[i * p[j]] = f[i] * (p[j] - 1);
    	    if (i % p[j] == 0) break;
    	}
    }
    for (int i = 1; i <= n; i++) // 前缀和
    {
    	s[i] = s[i - 1] + f[i];
    }
    return 0;
}

正式求答案

有了上面这些后，我们枚举每个\(x(1\leq x\leq n)\)，求出有多少对数字\((i,j)\)满足\(\gcd(i,j) = x\)。

这部分对于某个\(x\)的求法是：先求出有多少对\((p,q)\)数字满足\(p,q \leq \frac{n}{x} 且 \gcd(p,q) = 1\)，这个的对数和\((i,j)\)满足\(\gcd(i,j) = x\)的对数相同，原因是对于每一对\(p,q\)，都能存在一对\((i,j)\)满足要求，且\(i=xp,j=xq\)

那么如果考虑\((p,q)\)数字满足\(p,q \leq \frac{n}{x} 且 \gcd(p,q) = 1 且 p \leq q\),那么对数是\(f_1 + f_2 + f_3 + ... + \frac{n}{x}\)，也就是\(q\)的所有方案中\(p\)的方案和。
而\(f_1 + f_2 + f_3 + ... + \frac{n}{x} = s_{\frac{n}{x}}\),所以如果\(p \leq q（或者i \leq j)\)，对数就是s[n / x]

\(p \geq q\)的方案同理（具有对称性），也有s[n / x]对。

但是如果直接求s[n / x] * 2，那么\((p,q)=(1,1)\)的方案会被算两次，所以方案数为s[n / x] * 2 - 1，

所以有s[n / x] * 2 - 1对数字\((i,j)\)满足\(\gcd(i,j) = x\)。所以累加时\(x\)会被算s[n / x] * 2 - 1次，和为(s[n / x] * 2 - 1) * x

这一部分代码：

long long ans = 0; // 记录所有x的答案和
for (int x = 1; x <= n; x++)
{
	ans += (2 * s[n / x] - 1) * x % mod; // 按照公式累加，记得取模
	ans %= mod; // 上面只是把x的和取模了，接下来要把已统计下来的ans整个取模。
}

代码

最终的代码：

#include<bits/stdc++.h>
using namespace std;
const int mod = 998244353;
const int N = 5e7 + 5;
int p[N], f[N];
long long s[N];
bool b[N];
int main()
{
    int n;
    cin >> n;
    f[1] = 1, f[2] = 1;
    int c = 0;
    for (int i = 2; i <= n; i++)
    {
    	if (!b[i])
    	{
    		p[++c] = i;
    		f[i] = i - 1;
    	}
    	for (int j = 1; j <= c && i * p[j] <= n; j++)
    	{
    	    b[i * p[j]] = true;
    	    if (i % p[j] == 0) f[i * p[j]] = f[i] * p[j];
    	    else f[i * p[j]] = f[i] * (p[j] - 1);
    	    if (i % p[j] == 0) break;
    	}
    }
    for (int i = 1; i <= n; i++)
    {
    	s[i] = s[i - 1] + f[i];
    }
    long long ans = 0;
    for (int x = 1; x <= n; x++)
    {
    	ans += (2 * s[n / x] - 1) * x % mod;
    	ans %= mod;
    }
    cout << ans << endl;
	return 0;
}

完成~~~~