题意
给出\(n(1\leq n \leq 10^9)\),求\(\sum_{i=1}^n\sum_{j=1}^n\sigma(ij)\),其中\(\sigma(n)\)表示\(n\)的约数之和。
balabala
交了两道杜教筛的的板子题(51nod 1239, 1244)就看到了这题，然后不会搞，然后看题解看了一天一夜终于彻底搞明白一发A掉了。。。感觉学到了很多，写个博客整理一下，如有错请指出。
技能需求
数论函数与线性筛
莫比乌斯反演（也可以当成容斥去理解）
狄利克雷卷积
杜教筛
强大的数学推导功力（我没有）
先贴参考链接：
题目：http://www.51nod.com/Challenge/Problem.html#problemId=1220
这位这题的题解写的很好: https://www.cnblogs.com/hzoier/p/6937717.html
杜教筛：https://www.cnblogs.com/peng-ym/p/9446555.html 这位写的一看就懂
题解
推推推推推推推，我们记答案为\(Ans\)
1
我们一上来就需要用到一个公式，就是下面这个
\[
\sigma(ij)=\sum_{a|i}\sum_{b|j}[gcd(a,b)=1]\frac{aj}{b}
\]
感觉网上那群写博客的都是一眼就能看懂这是啥意思，但我却想了几个小时。。。所以我写详细一点
我们有\(\sigma(n)=\prod_{i=1}^k(p_i^0+p_i^1+...+p_i^{e_i})\),其中不同的\(p\)是出现在\(n\)的质因数分解式中不同的质因子，所以我们设\(i=p_1^{a_1}p_2^{a_2}...p_k^{a_k},j=p_1^{b_1}p_2^{b_2}...p_k^{b_k}(a,b\geq 0)\)
则\(\sigma(ij)=\prod_{m=1}^k(p_m^0+p_m^1+...+p_m^{a_m+b_m})\)
现在要想想上面那个公式的含义了，我们不妨考虑\(k=2\)的情况，假设我们选了\(a=p_1^{x_1},b=p_2^{x_2}\)，这时\(gcd(a,b)=1\),那么当\(x_1\)从\(1\)往\(a_1\)增大时，\(x_2=0\),则\(\frac{aj}{b}\)中\(p_1\)的幂次范围是\(1+b_1\)到\(a_1+b_1\);反过来，当\(x_2\)从\(1\)往\(b_1\)增大时，\(x_1=0\),则\(\frac{aj}{b}\)中\(p_1\)的幂次范围是\(b_1-1\)到\(0\).再加上\(x_1=x_2=0\)的情况，我们正好不重不漏的选出了\(p_1\)的所有可能选择。
回头考虑\(\sigma(ij)=\prod_{m=1}^k(p_m^0+p_m^1+...+p_m^{a_m+b_m})\),如果我们把这\(k\)个因式乘出来，变成一个多项式，从这个多项式的角度的某一项去考虑，它有\((a_1+b_1)(a_2+b_2)...(a_k+b_k)\)种选择，选取多项式\(i\)的哪一项不会对选取多项式\(j\)产生任何影响，这就是它们”互相独立“的含义。
2
下面这一段主要思路就是不停的交换求和符号，把整除条件改成枚举条件以化简式子。举个例子，我们要计算下面这个式子：
\[
f(n)=\sum_{i=1}^n\sum_{j|i}i^2j
\\=\sum_{j=1}^n\sum_{j|i}i^2j
\\=\sum_{j=1}^n\sum_{a=1}^{\lfloor\frac{n}{j}\rfloor}(aj)^2j
\\=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2i^3
\\=\sum_{i=1}^ni^3\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2
\]
回到正题：
\[
Ans=\sum_{i=1}^n\sum_{j=1}^n\sum_{a|i}\sum_{b|j}[gcd(a,b)=1]\frac{aj}{b}
\\=\sum_{a=1}^n\sum_{b=1}^n\sum_{p=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{b}\rfloor}[gcd(a,b)=1]\frac{aqb}{b}
\\=\sum_{a=1}^n\sum_{b=1}^n\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{b}\rfloor}[gcd(a,b)=1]aj
\\=\sum_{a=1}^na\sum_{b=1}^n[gcd(a,b)=1]\sum_{i=1}^{\lfloor\frac{n}{a}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{b}\rfloor}j
\\=\sum_{a=1}^na\lfloor\frac{n}{a}\rfloor\sum_{b=1}^n[gcd(a,b)=1]\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2}
\]
接下来还要搞的话只能在\([gcd(a,b) = 1]\)上做文章了
3
公式:\([gcd(a,b)=1]=\sum_{d=1}^{min(a,b)}\mu(d)[d|gcd(a,b)]\)
原理你可以当成莫比乌斯反演去搞，也可以当成容斥原理来想，（因为这两个本质上是一个东西----huchi）当\(gcd=1\)时，右边等于\(1\)，当\(gcd>1\)时，以\(gcd=12=2^2*3\)为例，右边\(=\mu(1)+\mu(2)+\mu(3)+\mu(4)+\mu(6)+\mu(12)\)，去掉为\(0\)的值，按质因数的个数分个组,\(=(\mu(1))+(\mu(2)+\mu(3))+(\mu(6))=C_2^01^0(-1)^2+C_2^11^1(-1)^1+C_2^21^2(-1)^0=(1-1)^2=0\)
以此类推，如果\(gcd\)有\(x\)个质因子那么就是\((1-1)^x\)，这个公式是这么来的，如果需要严谨的证明的话我不会或许应该有人写过
4
继续。。
\[
Ans=\sum_{a=1}^na\lfloor\frac{n}{a}\rfloor\sum_{b=1}^n\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2}\sum_{d=1}^n\mu(d)[d|gcd(a,b)]
\]
老套路，继续交换求和符号，这次要连着\([d|gcd(a,b)]\)一起考虑。\(d|gcd(a,b)\)等价于\(d|a,d|b\).
\[
Ans=\sum_{d=1}^n\sum_{d|a}\sum_{d|b}a\lfloor\frac{n}{a}\rfloor\frac{\lfloor\frac{n}{b}\rfloor(\lfloor\frac{n}{b}\rfloor+1)}{2}\mu(d)
\\=\sum_{d=1}^n\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}pd\lfloor\frac{n}{pd}\rfloor\frac{\lfloor\frac{n}{pd}\rfloor(\lfloor\frac{n}{pd}\rfloor+1)}{2}\mu(d)
\\=\sum_{d=1}^nd\mu(d)\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}p\lfloor\frac{n}{pd}\rfloor\sum_{q=1}^{\lfloor\frac{n}{d}\rfloor}\frac{\lfloor\frac{n}{qd}\rfloor(\lfloor\frac{n}{qd}\rfloor+1)}{2}
\]
咦，后两个求和符号那两个式子好像是独立的，那就不妨设
\(f(n)=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\)
\(=\sum_{i=1}^n\sum_{j=1}^ni[i|j]\) 因为\(\lfloor\frac{n}{i}\rfloor\)意思就是\(1\)到\(n\)中\(i\)的倍数的个数。
\(=\sum_{j=1}^n\sum_{i=1}^ni[i|j]\)
\(=\sum_{j=1}^n\sigma(j)\).
\(g(n)=\sum_{i=1}^n\frac{\lfloor\frac{n}{i}\rfloor(\lfloor\frac{n}{i}\rfloor+1)}{2}\) 之前好像推出过这个东西，我们再把它变回去试试
\(=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j\)
\(=\sum_{j=1}^n\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}j\)
\(=\sum_{j=1}^n\sum_{i=1}^{n}j[j|i]\)
\(=f(n)\).
我们把\(f(n)\)记成\(S_{\sigma}(n)\)吧.
因此有
\[
Ans=\sum_{d=1}^nd\mu(d)(S_{\sigma}(\lfloor\frac{n}{d}\rfloor))^2
\]
5
套杜教筛公式的时间到啦
\(\sigma\)是积性函数，可以线筛可以杜教筛，做法很多，我是卷上了个\(\mu\)得到\(id\)来筛的。
记\(f(d)=d\mu(d)\),\(f\)也是个积性函数，可以线筛可以杜教筛，卷上\(id\)得到\(e\),然后杜教筛起来就完事了。
最后求\(Ans\)的时候分块求，需要啥前缀和了就\(Calc\)啥，反正都是\(O(n^{\frac{2}{3}})\)我也不会分析
代码
\(main\)里是主分块，\(Calcmu\)是\(\mu\)的筛，\(Calcf\)是\(f\)的筛，\(Calcs\)是\(\sigma\)的筛
\(359ms\)
#include <bits/stdc++.h>
#define MAXN 1000000
#define LL long long
using namespace std;
const int Mx = 10000;
const int inv2 = 500000004;
const int mod = 1000000007;
LL gcd(LL a, LL b)
{
return b?gcd(b,a%b):a;
}
int prime[MAXN + 8];
bool notprime[MAXN + 8];
int mu[MAXN + 8];
LL fsum[MAXN + 8], ssum[MAXN + 8], dsum[MAXN + 8];
int ind = 0;
void getprime()
{
mu[1] = 1;
ssum[1] = dsum[1] = 1;
for (int i = 2; i <= MAXN; ++i)
{
if(!notprime[i])
{
prime[++ind] = i;
mu[i] = -1;
ssum[i] = dsum[i] = i + 1;
}
for (int j = 1; j <= ind && i * prime[j] <= MAXN; ++j)
{
notprime[i * prime[j]] = 1;
if(i % prime[j] == 0)
{
mu[i * prime[j]] = 0;
ssum[i * prime[j]] = ssum[i] / dsum[i] * (dsum[i] * prime[j] + 1); dsum[i * prime[j]] = dsum[i] * prime[j] + 1;
break;
}
mu[i * prime[j]] = -mu[i];
ssum[i * prime[j]] = ssum[i] * (1 + prime[j]); dsum[i * prime[j]] = 1 + prime[j];
}
}
for (int i = 1; i <= MAXN; ++i)
{
fsum[i] = (fsum[i - 1] + i * mu[i] % mod + mod) % mod;
ssum[i] = (ssum[i - 1] + ssum[i] % mod) % mod;
mu[i] = (mu[i - 1] + mu[i] + mod) % mod;
}
}
unordered_map<int,LL> _musum, _fsum, _ssum;
LL Calcmu(int n)
{
if (n <= MAXN)
return mu[n];
auto it = _musum.find(n);
if (it != _musum.end())
return it->second;
int last;
LL ret = 1;
for (int i = 2; i <= n; i = last + 1)
{
last = n / (n / i);
ret = (ret - (1LL * last - i + 1) * Calcmu(n / i) % mod + mod) % mod;
}
return _musum[n] = ret;
}
LL Calcf(int n)
{
if (n <= MAXN)
return fsum[n];
auto it = _fsum.find(n);
if (it != _fsum.end())
return it->second;
int last;
LL ret = 1;
for (int i = 2; i <= n; i = last + 1)
{
last = n / (n / i);
ret = (ret - ((1LL * last + i) % mod * (last - i + 1) % mod * inv2 % mod * Calcf(n / i) % mod) + mod) % mod;
}
return _fsum[n] = ret;
}
LL Calcs(int n)
{
if (n <= MAXN)
return ssum[n];
auto it = _ssum.find(n);
if (it != _ssum.end())
return it->second;
int last;
LL ret = 1LL * n * (n + 1) % mod * inv2 % mod;
for (int i = 2; i <= n; i = last + 1)
{
last = n / (n / i);
ret = (ret - ((Calcmu(last) - Calcmu(i - 1) + mod) % mod * Calcs(n / i) % mod) + mod) % mod;
}
return _ssum[n] = ret;
}
int main()
{
getprime();
LL n, last;
cin >> n;
LL ans = 0;
for (int i = 1; i <= n; i = last + 1)
{
last = n / (n / i);
LL cc = Calcs(n / i);
ans = (ans + (Calcf(last) - Calcf(i - 1) + mod) % mod * cc % mod * cc % mod) % mod;
}
cout << ans << endl;
return 0;
}
原文地址：https://www.cnblogs.com/zhugezy/p/11312301.html