P3768 莫比乌斯反演 + 杜教筛

最新推荐文章于 2021-09-17 22:01:47 发布

敲代码的欧文

最新推荐文章于 2021-09-17 22:01:47 发布

阅读量224

点赞数

本文链接：https://blog.csdn.net/Irving0323/article/details/104425044

版权

题目传送门

题意：

计算 $\left ( \sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j) \right ) \; mod \; p$ 。

是素数。

数据范围： $1 \leqslant n \leqslant 10^{10} , 5 \cdot 10^8 \leqslant p \leqslant 1.1 \cdot 10^9$ 。

题解：

这道题主要是考杜教筛。

这道题是这篇博客的强化版，就接着说了。

$\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j) &= \sum_{d=1}^{n}d^3 \sum_{t=1}^{\left \lfloor \frac{n}{d} \right \rfloor} t^2\mu(t) (\sum_{i=1}^{\left \lfloor \frac{n}{td} \right \rfloor}i)^2 \end{aligned}$

设 $L(n) = (\sum_{i=1}^{n}i)^2$ 。

$\begin{aligned} \sum_{i=1}^{n} \sum_{j=1}^{n} lcm(i,j) &= \sum_{d=1}^{n}d^3 \sum_{t=1}^{\left \lfloor \frac{n}{d} \right \rfloor} t^2\mu(t) L(\left \lfloor \frac{n}{td} \right \rfloor) \end{aligned}$

可以整除分块套整除分块求解，现在的复杂度是 O(n) ，还需要优化，考虑常用套路换元。

令 T=td 。

$\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j) &= \sum_{T=1}^{n} L(\left \lfloor \frac{n}{T} \right \rfloor) \sum_{d\mid T}d^3(\frac{T}{d})^2\mu(\frac{T}{d}) \\ &= \sum_{T=1}^{n}T^2L(\left \lfloor \frac{n}{T} \right \rfloor) \sum_{d\mid T}d\mu(\frac{T}{d}) \end{aligned}$

然后有一个常用的连接欧拉函数和莫比乌斯函数的式子： $\varphi(n) = \sum_{d\mid n}d\mu(\frac{n}{d})$ 。

$\begin{aligned} \sum_{i=1}^{n}\sum_{j=1}^{n}ijgcd(i,j) &= \sum_{T=1}^{n}T^2\varphi(T)L(\left \lfloor \frac{n}{T} \right \rfloor) \end{aligned}$

设 $f(n) = \sum_{i=1}^{n}i^2\varphi(n)$ 。

然后考虑整除分块去做，是可以 O(1) 计算的，主要就是计算的前缀和。

因为前缀和必须低于线性计算，所以我们考虑杜教筛。

杜教筛常用套路：为了求 $\sum_{i=1}^{n} f(i)$ ，我们构造积性函数，然后令 h = f *g ，是狄利克雷卷积。

然后如果 $\sum_{i=1}^{n} h(i)$ 很容易计算的话，那我们就可以求出来 $\sum_{i=1}^{n} f(i)$ 了。

设 $S(n)=\sum_{i=1}^{n} f(i)$ 。

杜教筛公式： $g(1)S(n) = \sum_{i = 1}^{n}h(i) - \sum_{d=2}^{n}g(d)\cdot S(\left\lfloor \frac{n}{d} \right\rfloor)$

然后我们怎么构造积性函数呢？

先把写出来吧。

$\begin{aligned} h(n) &= f * g\\ &= \sum_{d\mid n} f(d) g(\frac{n}{d})\\ &= \sum_{d\mid n} d^2 \varphi(n) g(\frac{n}{d})\\ \end{aligned}$

然后其实是一个套路，我们考虑消除 d^2 的方法去构造。

令 g(n) = n^2 。

$\begin{aligned} h(n) &= n^2 \sum_{d\mid n} \varphi(n)\\ \end{aligned}$

然后又是一个常用式子： $\sum_{d\mid n}\varphi(d) = n$ 。

然后 h(n) = n ^ 3 。

$g(1)S(n) = \sum_{i = 1}^{n}h(i) - \sum_{d=2}^{n}g(d)\cdot S(\left\lfloor \frac{n}{d} \right\rfloor)$

转化为

$S(n) = \sum_{i = 1}^{n}h(i) - \sum_{d=2}^{n}d^2\cdot S(\left\lfloor \frac{n}{d} \right\rfloor)$

然后我们就可以递归求解 S(n) 了。

有两个需要优化的地方：

（1）你需要线性筛 5e6 以内的 S(n) ，这样会少递归很多次。

（2） map 存储递归求解出的 S(n) ，这样也会少递归很多次。

感受：

第一次做杜教筛，感觉有点像高中学的化学方程式配平。

代码：

#include<bits/stdc++.h>
using namespace std ;
typedef long long ll ;
typedef pair<int , int> pli ;
const int maxn = 5e6 + 5 ;
ll mod ;
ll six , two ;
map<ll , ll> s ; 
int phi[maxn] ;
ll pre[maxn] ;
int num = 5e6 ;
void eular(int n)
{
   phi[1] = 1 ;
   for(int i = 2 ; i <= n ; i ++)  phi[i] = i ;
   for(int i = 2 ; i <= n ; i ++)
     if(phi[i] == i)
       for(int j = i ; j <= n ; j += i)
         phi[j] = phi[j] / i * (i - 1) ; 
   for(int i = 1 ; i <= n ; i ++)
     pre[i] = pre[i - 1] + ll(i) * i % mod * phi[i] % mod ,
     pre[i] %= mod ;
}
ll qpow(ll a , ll b)
{
   ll ans = 1 ; 
   a %= mod ;
   while(b)
   {
     if(b & 1)  ans = (ans * a) % mod ;
     b >>= 1 , a = (a * a) % mod ;
   }
   return ans % mod ;
}
ll inv(ll n)
{
	return qpow(n , mod - 2) % mod ;
}
ll sum(ll n)
{
	n %= mod ;
	ll c = n * (n + 1) % mod ;
	c *= two , c %= mod ;
	return c * c % mod ;
}
ll sum2(ll n)
{
	n %= mod ;
	ll c = n * (n + 1) % mod * (2 * n + 1) % mod ;
	c *= six , c %= mod ;
	return c ;
}
ll dls(ll n)
{
	if(n <= num)  return pre[n] ;
	if(s[n])  return s[n] ;
	//cout << n << '\n' ;
	ll ans = sum(n) ;
	for(ll l = 2 , r ; l <= n ; l = r + 1)
	{
		r = n / (n / l) ;
		ll c = (sum2(r) - sum2(l - 1) + mod) % mod ;
		ans -= c * dls(n / l) ;
		ans %= mod ;
	}
	ans = (ans + mod) % mod ;
	s[n] = ans ;
	return ans ;
}
void solve(ll n)
{
	ll ans = 0 ;
	for(ll l = 1 , r ; l <= n ; l = r + 1)
	{
		r = n / (n / l) ;
		ll c = (dls(r) - dls(l - 1) + mod) % mod ;
		ans += c * sum(n / l) ;
		ans %= mod ;
	}
	printf("%lld\n" , ans) ;
}
int main()
{
	ll n , m ;
	scanf("%lld%lld" , &mod , &n) ;
	eular(num) ;
	six = inv(6ll) ;
	two = inv(2ll) ;
	solve(n) ;
	return 0 ;
}