[洛谷P3768]简单的数学题(杜教筛)

题面

https://www.luogu.com.cn/problem/P3768

题解

前置知识

• 杜教筛 https://www.cnblogs.com/xh092113/p/12289949.html
• 求和中消gcd技巧 https://www.cnblogs.com/xh092113/p/12268356.html

看到gcd，开始反演

\[{\sum\limits_{i=1}^{n}}{\sum\limits_{j=1}^{n}}ij(i,j) \]

\[={\sum\limits_{i=1}^{n}}{\sum\limits_{j=1}^{n}}ij{\sum\limits_{p|(i,j)}}{\mu(p)} \]

\[={\sum\limits_{i=1}^{n}}{\sum\limits_{j=1}^{n}}ij{\sum\limits_{p|i,p|j}}{\mu(p)} \]

\[={\sum\limits_{p}}{\mu(p)}{\sum\limits_{i=1}^{n}}{\sum\limits_{j=1}^{n}}ij[p|i][p|j] \]

\[={\sum\limits_{p}}{\mu(p)}*p^2*({\frac{{\lfloor}{\frac{n}{p}}{\rfloor}({\lfloor}{\frac{n}{p}}{\rfloor}+1)}{2}})^2① \]

然后使用数论分块，剩下的问题就只是求解\(f(p)={\mu(p)}*p^2\)的前缀和sf了。

令\(g(p)=p^2，h(p)=p^3\)（这里g,h的作用详见前置知识->杜教筛）就可以使用杜教筛求解了。值得注意的是，杜教筛可以一次性求出所有的\(sf({\lfloor}{\frac{n}{i}}{\rfloor})(1{\leq}i{\leq}n)\)，而数论分块求①式的时候，也恰好只用到了这些sf值（分块时每一块的右端点R就是通过n/(n/L)得到的，所以每一个右端点都是\({\lfloor}{\frac{n}{i}}{\rfloor}(1{\leq}i{\leq}n)\)中的一个），所以总体相当于只需要做一次杜教筛，就\(O(n^{\frac{2}{3}})\)解决本题。

代码

• P.S.一定要注意爆long long的问题。虽然\(p^2\)并不会爆，但是\(n^2\)会爆。
• 如果此题的p再大一些，使得\(p^2\)也会爆的话，就必须要使用O(1)快速乘了。

#include<bits/stdc++.h>

using namespace std;

#define D 4641600
#define M 2160
#define ll long long
#define rg register

ll mod;

namespace ModCalc{
	inline void Inc(ll &x,ll y){
		x += y;if(x >= mod)x -= mod;
	}

	inline void Dec(ll &x,ll y){
		x -= y;if(x < 0)x += mod;
	}
	
	inline ll Add(ll x,ll y){
		Inc(x,y);return x;
	}
	
	inline ll Sub(ll x,ll y){
		Dec(x,y);return x;
	}
}
using namespace ModCalc;

ll sphi[D+M+5],pri[D/5+5];
ll n,pn,v2,v3;
bool isp[D+5];

inline ll F(ll x){
	x = x % mod * ((x + 1) % mod) % mod * v2 % mod;
	return x * x % mod;
}

inline void Eular(){
	sphi[1] = 1;
	for(rg ll i = 2;i <= D;i++)isp[i] = 1;
	for(rg ll i = 2;i <= D;i++){  
		if(isp[i])pri[++pn] = i,sphi[i] = i - 1;
		for(rg ll j = 1;i * pri[j] <= D;j++){
			isp[i*pri[j]] = 0;
			if(i % pri[j])sphi[i*pri[j]] = sphi[i] * sphi[pri[j]];
			else{
				sphi[i*pri[j]] = sphi[i] * pri[j];
				break;
			}
		}
	}
	for(rg ll i = 1;i <= D;i++)sphi[i] = sphi[i] * i % mod * i % mod;
	for(rg ll i = 2;i <= D;i++)Inc(sphi[i],sphi[i-1]); 
}

inline ll id(ll x){
	if(x <= D)return x;
	return D + n / x;
} 

inline ll sg(ll x){
	return (x % mod) * ((x + 1) % mod) % mod * (((x << 1) | 1) % mod) % mod * v2 % mod * v3 % mod;
}

inline ll sumf(ll x){
	if(sphi[id(x)])return sphi[id(x)];
	ll ans = F(x),L,R = 1;
	while(R < x){
		L = R + 1,R = x / (x / L);
		Dec(ans,Sub(sg(R),sg(L-1)) * sumf(x/L) % mod);
	}
	return sphi[id(x)] = ans;	
}

int main(){
	scanf("%lld%lld",&mod,&n);
	Eular();
	v2 = (mod + 1) >> 1;
	v3 = mod % 3 == 1 ? ((mod << 1) | 1) / 3 : (mod + 1) / 3;
	ll ans = 0,L,R = 0;
	while(R < n){
		L = R + 1;
		R = n / (n / L);
		Inc(ans,Sub(sumf(R),sumf(L-1)) * F(n/L) % mod);
	}
	cout << ans << endl;
	return 0;
}