[Luogu P3768] 简单的数学题

最新推荐文章于 2022-12-22 11:30:41 发布

LPA20020220

最新推荐文章于 2022-12-22 11:30:41 发布

阅读量250

点赞数

CC 4.0 BY-SA版权

分类专栏：数学莫比乌斯反演文章标签：数学莫比乌斯反演

本文链接：https://blog.youkuaiyun.com/LPA20020220/article/details/84868985

数学同时被 2 个专栏收录

122 篇文章

订阅专栏

莫比乌斯反演

19 篇文章

订阅专栏

本文详细解析了一道洛谷平台上的数学题目，该题目要求计算特定条件下的模运算结果。通过对题目的深入分析，文章提供了从数学角度理解并解决此类问题的方法，包括使用最大公约数、分块和前缀和等技巧。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

洛谷传送门

题目描述

由于出题人懒得写背景了，题目还是简单一点好。

输入一个整数 $n$ 和一个整数 $p$ ，你需要求出 $(\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p$ ，其中 $g c d (a, b)$ 表示 $a$ 与 $b$ 的最大公约数。

刚才题面打错了，已修改

输入输出格式

输入格式：

一行两个整数 $p$ 、 $n$ 。

输出格式：

一行一个整数 $(\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))~mod~p$ 。

输入输出样例

输入样例#1：

998244353 2000

输出样例#1：

883968974

说明

对于20%的数据， $\leq 1000$ 。

对于30%的数据， $\leq 5000$ 。

对于60%的数据， $\leq 10^6$ ，时限1s。

对于另外20%的数据， $\leq 10^9$ ，时限3s。

对于最后20%的数据， $\leq 10^{10}$ ，时限6s。

对于100%的数据， $\times 10^8 \leq p \leq 1.1 \times 10^9$ 且 $p$ 为质数。

解题分析

大力推式子：
$\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j) \\ =\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d] \\ =\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ijd^2[gcd(i,j)=1] \\ =\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}ij\sum_{p|gcd(i,j)}\mu(p) \\ =\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2\sum_{i=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}ij$
设 $s(n)=\sum_{i=1}^ni$ ，并设 $T = d p$ ，那么：
$\sum_{d=1}^nd^3\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)p^2\sum_{i=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}ij \\ =\sum_{T=1}^{n}s(\frac{n}{T})^2T^2\sum_{d|T}d\mu(\frac{T}{d})$
而我们知道 $\phi(i)=\mu(i)*id(i)$ ，那么实际上这玩意就是：
$\sum_{T=1}^{n}s(\frac{n}{T})^2T^2\phi(T)$
前面可以下底分块，关键在于后面的 $f(x)=x^2\phi(x)$ 上面，这玩意该怎么科学筛出前缀和。

考虑到 $\phi(i)*I=id(i)$ ，我们可以设 $g(x)=x^2$ ，那么
$\\=\sum_{d|x}d^2*\phi(\frac{x}{d})(\frac{x}{d})^2 \\=\sum_{d|x}x^2\phi(x) \\ =x^3$
那么根据这个套路的式子：
$g(1)S(n)=\sum_{i=1}^nh(i)-\sum_{d=2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)$
和 $\sum_{i=1}^{n}h(i)=\frac{(x+1)^2x^2}{4},\sum_{i=1}^{n}g(i)=\frac{n*(n+1)*(2n+1)}{6}$ ，就可以科学筛出来了。

代码如下：

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <tr1/unordered_map>
#include <algorithm>
#define R register
#define IN inline
#define W while
#define gc getchar()
#define MX 5005000
#define ll long long
int pcnt;
std::tr1::unordered_map <ll, ll> mp;
int phi[MX], pri[MX];
ll sum[MX], val[MX], MOD, inv2, inv4, inv6;
bool npr[MX];
void get()
{
	phi[1] = 1, val[1] = 1;
	R int i, j, tar;
	for (i = 2; i <= 5e6; ++i)
	{
		if (!npr[i]) pri[++pcnt] = i, phi[i] = i - 1;
		for (j = 1; j <= pcnt; ++j)
		{
			tar = i * pri[j];
			if (tar > 5e6) break;
			npr[tar] = true;
			if (!(i % pri[j])) {phi[tar] = phi[i] * pri[j]; break;}
			phi[tar] = phi[i] * phi[pri[j]];
		}
	}
	for (i = 1; i <= 5e6; ++i) sum[i] = (sum[i - 1] + 1ll * i * i % MOD * phi[i] % MOD) % MOD;
}
IN ll f1(ll val)
{return val % MOD * val % MOD * (val + 1) % MOD * (val + 1) % MOD * inv4 % MOD;}
IN ll f2(ll val)
{return val % MOD * (val + 1) % MOD * (2 * val % MOD + 1) % MOD * inv6 % MOD;}
IN ll fpow(ll base, ll tim)
{
	ll ret = 1;
	W (tim)
	{
		if (tim & 1)ret = ret * base % MOD;
		base = base * base % MOD, tim >>= 1;
	}
	return ret;
}
ll solve(R ll bd)
{
	if (bd <= 5e6) return sum[bd];
	if (mp.count(bd)) return mp[bd];
	ll ret = f1(bd);
	for (R ll lef = 2, rig; lef <= bd; lef = rig + 1)
	{
		rig = bd / (bd / lef);
		ret -= ((f2(rig) - f2(lef - 1) + MOD) % MOD) * solve(bd / lef) % MOD;
	}
	return mp[bd] = (ret + MOD) % MOD;
}
int main(void)
{
	ll n, ans = 0, bd, mul;
	scanf("%lld%lld", &MOD, &n);
	get();
	inv2 = fpow(2, MOD - 2);
	inv4 = fpow(4, MOD - 2);
	inv6 = fpow(6, MOD - 2);
	for (R ll lef = 1, rig; lef <= n; lef = rig + 1)
	{
		rig = n / (n / lef); bd = n / lef;
		mul = (1 + bd) % MOD * bd % MOD * inv2 % MOD;
		mul = mul * mul % MOD;
		(ans += mul * ((solve(rig) - solve(lef - 1) + MOD) % MOD) % MOD) %= MOD;
	}
	printf("%lld", ans);
}