补充几个前置知识

莫比乌斯反演

φ∗1=id\*1=idφ∗1=id

杜教筛

杜教筛能够在Θ(n23)\Theta(n^\frac{2}{3})Θ(n32​)求出一个积性函数的前缀和。

给出一个积性函数f(i)f(i)f(i)，求∑i=1nf(i)\sum_{i=1}^{n}f(i)i=1∑n​f(i)

构造一个前缀和容易求的积性函数ggg，令h=f∗gh=f*gh=f∗g（其中∗*∗是狄利克雷卷积），使得hhh的前缀和也容易求，那么

∑i=1nh(i)=∑i=1n∑d∣ig(d)f(id)=∑d=1ng(d)∑i=1⌊nd⌋f(i)=∑d=1ng(d)S(⌊nd⌋)=g(1)S(n)+∑d=2ng(d)S(⌊nd⌋)\sum_{i=1}^{n}h(i) \\ =\sum_{i=1}^{n}\sum_{d|i}g(d)f(\frac{i}{d}) \\ =\sum_{d=1}^{n}g(d)\sum_{i=1}^{\\frac{n}{d}\}f(i) \\ =\sum_{d=1}^{n}g(d)S(\\frac{n}{d}\) \\ =g(1)S(n)+\sum_{d=2}^{n}g(d)S(\\frac{n}{d}\)i=1∑n​h(i)=i=1∑n​d∣i∑​g(d)f(di​)=d=1∑n​g(d)i=1∑⌊dn​⌋​f(i)=d=1∑n​g(d)S(⌊dn​⌋)=g(1)S(n)+d=2∑n​g(d)S(⌊dn​⌋)

所以

g(1)S(n)=∑i=1nh(i)−∑d=2ng(d)S(⌊nd⌋)g(1)S(n)=\sum_{i=1}^{n}h(i)-\sum_{d=2}^{n}g(d)S(\\frac{n}{d}\)g(1)S(n)=i=1∑n​h(i)−d=2∑n​g(d)S(⌊dn​⌋)

不难发现左边部分就是h(i)h(i)h(i)的前缀和，右边可以通过整除分块做。

实现的时候需要先预处理前Θ(n23)\Theta(n^\frac{2}{3})Θ(n32​)项，然后用map存下求出的已知项。时间复杂度Θ(n23)\Theta(n^\frac{2}{3})Θ(n32​)（不会证qwq）

int S(int n)
{
	if(n<=5000000) return sum[n];
	if(mp[n]!=0) return mp[n];
	int ans=g(n);
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans-=S(n/l)*(t(r)-t(l-1))%mod)%=mod;
	}
	return mp[n]=ans;
}

如何构造g？

我们的目的是让h=f∗gh=f*gh=f∗g的前缀和也容易求，所以可以先把(f∗g)(n)(f*g)(n)(f∗g)(n)写出来，看看h(n)h(n)h(n)长啥样，努力凑出一个可以求前缀和的形式。

例题 【模板】杜教筛 题目描述

给定一个正整数，求

ans1=∑i=1nφ(i)ans_1=\sum_{i=1}^n\(i)ans1​=i=1∑n​φ(i)

ans2=∑i=1nμ(i)ans_2=\sum_{i=1}^n \mu(i)ans2​=i=1∑n​μ(i)

输入格式

本题单测试点内有多组数据。

输入的第一行为一个整数，表示数据组数TTT。

接下来TTT 行，每行一个整数nnn，表示一组询问。

输出格式

对于每组询问，输出一行两个整数，分别代表​ 和​。

样例 #1 样例输入 #1

6
1
2
8
13
30
2333

样例输出 #1

1 1
2 0
22 -2
58 -3
278 -3
1655470 2

提示 数据规模与约定

对于全部的测试点，保证1≤T≤101 \leq T \leq 101≤T≤10，1≤n 题解 莫比乌斯函数前缀和

μ∗1=εμ∗1=εμ∗1=ε

不难想到令g=1g=1g=1，那么h=εh=εh=ε

所以式子就变成了

1×S(n)=1−∑d=2nS(⌊nd⌋)1\times S(n)=1-\sum_{d=2}^{n}S(\\frac{n}{d}\)1×S(n)=1−d=2∑n​S(⌊dn​⌋)

欧拉函数前缀和

方法一：根据定义用莫比乌斯反演做。

∑i=1nφ(i)=∑i=1n∑j=i+1n[gcd(i,j)=1]\sum_{i=1}^{n}\(i)=\sum_{i=1}^{n}\sum_{j=i+1}^{n}[gcd(i,j)=1]i=1∑n​φ(i)=i=1∑n​j=i+1∑n​[gcd(i,j)=1]

这不就是莫比乌斯反演板子题吗

方法二：杜教筛

φ∗1=id\*1=idφ∗1=id

令g=1g=1g=1得

S(n)=∑i=1ni−∑i=2nS(⌊nd⌋)S(n)=\sum_{i=1}^{n}i-\sum_{i=2}^nS(\\frac{n}{d}\)S(n)=i=1∑n​i−i=2∑n​S(⌊dn​⌋)

代码

#include
#define int long long
using namespace std;
const int N=1e6+7;
int mu[N],phi[N],s1[N],s2[N];
bool bz[N];
vector p;
map mp1,mp2;
void init(int n)
{
	mu[1]=1; phi[1]=1;
	for(int i=2; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i);
			mu[i]=-1; phi[i]=i-1;
		}
		for(int j:p)
		{
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0)
			{
				phi[i*j]=phi[i]*j;
				break;
			}
			mu[i*j]=-mu[i];
			phi[i*j]=phi[i]*(j-1);
		}
	}
	for(int i=1; i<=n; i++)
	{
		s1[i]=mu[i]+s1[i-1];
		s2[i]=phi[i]+s2[i-1];
	}
}
int S1(int n)
{
	if(n<2) return n;
	if(n<=N-7) return s1[n];
	if(mp1[n]) return mp1[n];
	int ans=1;
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		ans-=S1(n/l)*(r-l+1);
	}
	return mp1[n]=ans;
}
int S2(int n)
{
	if(n<2) return n;
	if(n<=N-7) return s2[n];
	if(mp2[n]) return mp2[n];
	int ans=n*(n+1)/2;
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		ans-=S2(n/l)*(r-l+1);
	}
	return mp2[n]=ans;
}
void O_o()
{
	int n;
	cin>>n;
	cout<>T;
	while(T--)
	{
		O_o();
	}
}

简单的数学题 题目描述

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

输入一个整数nnn 和一个整数ppp，你需要求出：

(∑i=1n∑j=⁡(i,j)) mod p\left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p(i=1∑n​j=1∑n​ijgcd(i,j))modp

其中gcd⁡(a,b)\gcd(a,b)gcd(a,b) 表示aaa 与bbb 的最大公约数。

输入格式

一行两个整数p,np,np,n。

输出格式

一行一个整数表示答案。

样例 #1 样例输入 #1

998244353 2000

样例输出 #1

883968974

提示

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

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

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

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

对于最后20%的数据，n≤1010n \leq 10^{10}n≤1010，时限4s。

对于100%的数据，5×108≤p≤1.1×1095 \times 10^8 \leq p \leq 1.1 \times 10^95×108≤p≤1.1×109 且ppp 为质数。

题解

引入欧拉反演（用莫比乌斯反演能做，推出来是一样的，但是很麻烦）

∑i=1n∑j=1ngcd(i,j)ij=∑i=1n∑j=1nid(gcd(i,j))ij=∑i=1n∑j=1n∑d∣gcd(i,j)φ(d)ij=∑d=1nφ(d)∗d2∑i=1⌊nd⌋∑j=1⌊nd⌋ij\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)ij \\ =\sum_{i=1}^{n}\sum_{j=1}^{n}id(gcd(i,j))ij \\ =\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{d|gcd(i,j)}\(d)ij \\ =\sum_{d=1}^{n}\(d)*d^2\sum_{i=1}^{\\frac{n}{d}\}\sum_{j=1}^{\\frac{n}{d}\}iji=1∑n​j=1∑n​gcd(i,j)ij=i=1∑n​j=1∑n​id(gcd(i,j))ij=i=1∑n​j=1∑n​d∣gcd(i,j)∑​φ(d)ij=d=1∑n​φ(d)∗d2i=1∑⌊dn​⌋​j=1∑⌊dn​⌋​ij

令sum(n)=∑i=1n∑j=(n)=\sum_{i=1}^{n}\sum_{j=1}^{n}ijsum(n)=∑i=1n​∑j=1n​ij

显然sum(n)=n2(n+1)24sum(n)=\frac{n^2(n+1)^2}{4}sum(n)=4n2(n+1)2​

令f(i)=φ(i)∗i2f(i)=\(i)*i^2f(i)=φ(i)∗i2，现在我们如果能求出它的前缀和，我们就可以整除分块解决问题。这个时候就要杜教筛了：

(f∗g)(n)=∑d∣nφ(d)∗d2∗g(nd)(f*g)(n)\\=\sum_{d|n}\(d)*d^2*g(\frac{n}{d})(f∗g)(n)=d∣n∑​φ(d)∗d2∗g(dn​)

d2d^2d2特别烦人怎么办？消掉！

令g(i)=i2g(i)=i^2g(i)=i2

(f∗g)(n)=∑d∣nφ(d)∗d2∗g(nd)=n2∑d∣nφ(d)=n3(f*g)(n)\\=\sum_{d|n}\(d)*d^2*g(\frac{n}{d})\\=n^2\sum_{d|n}\(d)\\=n^3(f∗g)(n)=d∣n∑​φ(d)∗d2∗g(dn​)=n2d∣n∑​φ(d)=n3

用我们小学二年级就学过的立方和公式得

∑i=1ni3=n2(n+1)24=sum(n)\sum_{i=1}^{n}i^3=\frac{n^2(n+1)^2}{4}=sum(n)i=1∑n​i3=4n2(n+1)2​=sum(n)

于是乎，整除分块即可

代码

#include
#define int long long
using namespace std;
const int N=5e6+7;
int mod,inv6,phi[N],sum[N];
bool bz[N];
vector p;
map mp;
int power(int x,int t)
{
	int b=1;
	while(t)
	{
		if(t&1) b=b*x%mod;
		x=x*x%mod; t>>=1;
	}
	return b;
}
void init(int n)
{
	phi[1]=1;
	for(int i=2; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i); phi[i]=i-1;
		}
		for(int j:p)
		{
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0)
			{
				phi[i*j]=phi[i]*j%mod;
				break;
			}
			phi[i*j]=phi[i]*phi[j]%mod;
		}
	}
	for(int i=1; i<=n; i++) sum[i]=(sum[i-1]+phi[i]*i%mod*i%mod)%mod;
}
int t(int n)
{
	n%=mod;
	return n*(n+1)%mod*(2*n%mod+1)%mod*inv6%mod;
}
int g(int n)
{
	n%=mod;
	return (n*(n+1)/2%mod)*(n*(n+1)/2%mod)%mod;
}
int S(int n)
{
	if(n<=5000000) return sum[n];
	if(mp[n]!=0) return mp[n];
	int ans=g(n);
	for(int l=2,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans-=S(n/l)*(t(r)-t(l-1))%mod)%=mod;
	}
	return mp[n]=ans;
}
void O_o()
{
	int n,ans=0;
	cin>>mod>>n;
	inv6=power(6,mod-2);
	init(5000000);
	for(int l=1,r; l<=n; l=r+1)
	{
		r=n/(n/l);
		(ans+=(S(r)-S(l-1))*g(n/l)%mod)%=mod;
	}
	cout<<(ans+mod)%mod;
}
signed main()
{
	ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
	cout<>T;
	while(T--)
	{
		O_o();
	}
}