[Luogu P3768] 简单的数学题

洛谷传送门

题目描述

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

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

刚才题面打错了，已修改

输入输出格式

输入格式：

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

输出格式：

一行一个整数$({\sum }_{i=1}^n{\sum }_{j=1}^{n}ijgcd(i,j))modp$。

输入输出样例

输入样例#1：

998244353 2000

输出样例#1：

883968974

说明

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

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

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

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

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

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

解题分析

大力推式子：
$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j) \\ =\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)=d] \\ =\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij{d}^2[gcd(i,j)=1] \\ =\sum _{d=1}^n{d}^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}^n{d}^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 \end{gathered}$$
设$s(n)={\sum }_{i=1}^{n}i$， 并设$T=dp$， 那么：
$$\begin{gathered} \sum _{d=1}^n{d}^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}{)}^2{T}^2\sum _{d|T}d\mu (\frac{T}{d}) \end{gathered}$$
而我们知道$\varphi (i)=\mu (i)\ast id(i)$， 那么实际上这玩意就是：
$$\sum _{T=1}^{n}s(\frac{n}{T}{)}^2{T}^2\varphi (T)$$
前面可以下底分块， 关键在于后面的$f(x)=x^2\varphi (x)$上面， 这玩意该怎么科学筛出前缀和。

考虑到$\varphi (i)\ast I=id(i)$， 我们可以设$g(x)=x^2$， 那么
$$\begin{gathered} h(x)=g(x)\ast f(x) \\ =\sum _{d|x}{d}^2\ast \varphi (\frac{x}{d})(\frac{x}{d}{)}^2 \\ =\sum _{d|x}{x}^2\varphi (x) \\ =x^3 \end{gathered}$$
那么根据这个套路的式子：
$$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$$
和${\sum }_{i=1}^{n}h(i)=\frac{(x+1{)}^2{x}^2}{4},{\sum }_{i=1}^{n}g(i)=\frac{n\ast (n+1)\ast (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);
}