洛谷P3768简单的数学题——莫比乌斯反演+杜教筛推导详解

最新推荐文章于 2023-08-14 06:30:00 发布

偶耶XJX

最新推荐文章于 2023-08-14 06:30:00 发布

阅读量229

点赞数

分类专栏：信息竞赛解题文章标签：数学

本文链接：https://blog.csdn.net/weixin_43960287/article/details/113483123

版权

信息竞赛解题专栏收录该内容

195 篇文章 5 订阅

订阅专栏

前置技能：数论基础、莫比乌斯反演、狄利克雷卷积、杜教筛

题目

传送门

题解

~~（由于笔者对富文本编辑器情有独钟，下面的推导公式可能都是以无法复制的图片展示的）~~

题意很简单，求 $ans=\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)$ ，模数p是≥5e8的素数，所以放心求逆元；

下面开始推式子：

发现有个gcd，没有好的方法直接快速求和，所以把gcd提出来：

$ans=\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)=\sum_{d=1}^nd^3\sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\sum_{j=1}^{\left\lfloor \frac{n}{d} \right\rfloor}ij[gcd(i,j)=1]$

后两个 $\sum_$ 的部分很熟悉，我们发现可以用莫比乌斯反演，设 $mf(x)=\sum_{i=1}^{n}\sum_{j=1}^{n}ij[gcd(i,j)=x]$ ， $mF(x)=\sum_{x|d}mf(x)$ （用m-大小f表示莫比乌斯反演推导的函数是为了与后面做区分），

感性分析可得， $mF(x)=\sum_{x|d}mf(x)=x^2(\sum_{i=1}^{\left\lfloor \frac{n}{x} \right\rfloor}i)^2=x^2sum(\left\lfloor \frac{n}{x} \right\rfloor)^2$ （定义 $sum(x)=\sum_{i=1}^xi=\frac{x(x+1)}{2}$ ），

所以有 $mf(x)=\sum_{x|d}\mu(\frac{d}{x})mF(d)=\sum_{x|d}\mu(\frac{d}{x})d^2sum(\left\lfloor \frac{n}{d} \right\rfloor)^2$ ，

将 mf(1) 带入原式子： $ans=\sum_{d=1}^nd^3\sum_{i=1}^{\left\lfloor \frac{n}{d} \right\rfloor}\mu(i)i^2sum(\left\lfloor \frac{n}{di} \right\rfloor)^2$ ，

然后我们发现循环里有个 $\mu$ ，已目前能力~~（是我太菜了）~~无法直接求和，除非把 i 变为某个数的因子形式，所以改变枚举顺序，不枚举 d、i，而枚举 d*i 和 d；

令 $T=di,i=\frac{T}{d}$ ，有 $ans=\sum_{T=1}^nsum(\left\lfloor \frac{n}{T} \right\rfloor)^2\sum_{d|T}d^3\frac{T^2}{d^2}\mu(\frac{T}{d})=\sum_{T=1}^nsum(\left\lfloor \frac{n}{T} \right\rfloor)^2T^2\sum_{d|T}d\mu(\frac{T}{d})$ ；

观看尾部的 $\sum_{d|T}d\mu(\frac{T}{d})$ ，我们又双叒发现有点熟悉，因为根据刚学过的狄利克雷卷积的知识， $\left\{\begin{matrix} \phi*I=id\\ \mu*I=\epsilon \end{matrix}\right. \Rightarrow \phi=\mu*id$ ，也就是说 $\phi(x)=\sum_{d|x}d\mu(\frac{x}{d})$ ，

这就很棒了，因为我们已经可以把带 $\mu$ 的一堆换成欧拉函数： $ans=\sum_{T=1}^nsum(\left\lfloor \frac{n}{T} \right\rfloor)^2T^2\sum_{d|T}d\mu(\frac{T}{d})=\sum_{T=1}^nsum(\left\lfloor \frac{n}{T} \right\rfloor)^2T^2\phi(T)$ ，

设 $f(x)=x^2\phi(x)$ ，则 $ans=\sum_{T=1}^nsum(\left\lfloor \frac{n}{T} \right\rfloor)^2f(T)$ ；

显然， $\left\lfloor \frac{n}{T} \right\rfloor$ 最多只有 $2\sqrt{n}$ 种取值，如果我们能较快求得 f 函数的前缀和，就可以用数论分块做；

所以我们用一用杜教筛：

设 $S(n)=\sum_{i=1}^nf(i)$ ， $h(i)=(f*g)(i)=\sum_{d|i}f(d)g(\frac{i}{d})$ ，

$S(n)=(\sum_{i=1}^nh(i)-\sum_{i=2}^ng(i)S(\left\lfloor \frac{n}{i} \right\rfloor))/g(1)$ ；

根据套路，g 函数要根据 f 函数设置为最适合的一个完全积性函数，此处将 g 设置为 g(x)=id^2(x)=x^2 ，则 $h(i)=\sum_{d|i}\phi(d)d^2\frac{i^2}{d^2}=i^2\sum_{d|i}\phi(d)$ ，

由 $\phi*I=id$ ，也就是 $\sum_{d|i}\phi(d)=id(i)=i$ ，可以得到 h(i)=i^3 ；

由此我们得到了快速求 h 函数前缀和的方法，也就是立方前缀和公式： $\sum_{i=1}^ni^3=(\sum_{i=1}^ni)^2=sum(n)^2$ ，

然后套杜教筛： $S(n)=\sum_{i=1}^ni^3-\sum_{i=1}^ni^2S(\left\lfloor \frac{n}{i} \right\rfloor)=sum(n)^2-\sum_{i=1}^ni^2S(\left\lfloor \frac{n}{i} \right\rfloor)$ ，（平方前缀和公式： $\sum_{i=1}^ni^2=\frac{n(n+1)(2n+1)}{6}$ ）

此时再用一个数论分块，就可以较快求得 f 函数的前缀和，

最后把两个过程套起来，根据杜教筛的时间复杂度分析，全部平摊下来复杂度并不会达到 O(n) ，是 $O(n^{\frac{2}{3}})$ 的。

代码

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<tr1/unordered_map>
#define ll long long
#define MAXN 10000005
using namespace std;
using namespace tr1;
inline ll read(){
	ll x=0;bool f=1;char s=getchar();
	while((s<'0'||s>'9')&&s>0){if(s=='-')f^=1;s=getchar();}
	while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+s-'0',s=getchar();
	return f?x:-x;
}
ll n,p[MAXN],sp[MAXN],MOD=read();
unordered_map<ll,ll>ph;
bool nop[MAXN];
vector<ll>h;
inline ll ksm(ll a,ll b){
	ll res=1;
	for(;b;b>>=1){
		if(b&1)res=res*a%MOD;
		a=a*a%MOD;
	}
	return res;
}
ll NI=ksm(6,MOD-2);//提前求逆元
inline void build(){
	nop[0]=nop[1]=1;
	for(int i=1;i<MAXN-4;i++)p[i]=i;
	for(int a=2;a<MAXN-4;a++){
		if(!nop[a])h.push_back(a),p[a]=a-1;
		for(int i=0,u;i<h.size()&&h[i]*a<MAXN-4;i++){
			u=a*h[i],nop[u]=1;
			if(a%h[i]==0)p[u]=p[a]*h[i],i=MAXN;
			else p[u]=p[a]*p[h[i]];
		}
	}
	sp[1]=p[1];
	for(int i=2;i<MAXN-4;i++)sp[i]=(sp[i-1]+p[i]*i%MOD*i%MOD)%MOD;
}
inline ll sumx(ll n){n%=MOD;//n的范围比模数大就离谱
	return (n*(n+1)>>1)%MOD;   //一次前缀和（sum）
}
inline ll sumy(ll n){n%=MOD;
	return n*(n+1)%MOD*(n*2+1)%MOD*NI%MOD;  //平方前缀和
}
inline ll sumphi(ll n){
	if(n<MAXN-4)return sp[n];
	if(ph.find(n)!=ph.end())return ph[n];
	ll res=sumx(n)*sumx(n)%MOD;
	for(ll i=2,ls;i<=n;i=ls+1)
		ls=n/(n/i),res=(res-(sumy(ls)-sumy(i-1)+MOD)%MOD*sumphi(n/i)%MOD+MOD)%MOD;
	return ph[n]=res;
}
inline ll getans(ll n){
	ll res=0;
	for(ll i=1,ls;i<=n;i=ls+1)
		ls=n/(n/i),res=(res+(sumphi(ls)-sumphi(i-1)+MOD)%MOD*sumx(n/i)%MOD*sumx(n/i)%MOD)%MOD;
	return res;
}
int main()
{
	build();n=read();
	printf("%lld\n",getans(n));
	return 0;
}