51Nod1220：约数之和
题意：
(d(k))表示(k)所有约数的和。
比如说(d(6)=1+2+3+6=12)。
定义：(S(N)=sum_{i=1}^Nsum_{j=1}^Nd(i*j))。
给出(Nleq 10^9)，求(S(N))。
思路：
我们知道：(sigma(n))表示约数个数和，且为积性函数。
我们需要求(sigma(ij))，但是(sigma)不是完全积性函数。
这里有个结论：
[sigma(ij)=sum_{x|i}sum_{y|j}[gcd(x,y)=1]frac{yi}{x}
]
代入原式：
[sum_{i=1}^Nsum_{j=1}^Nsum_{x|i}sum_{y|j}[gcd(x,y)=1]frac{yi}{x}
]
反演：
[sum_{i=1}^Nsum_{j=1}^Nsum_{x|i}sum_{y|j}frac{yi}{x}sum_{d|gcd(x,y)}mu(d)\sum_{i=1}^Nsum_{j=1}^Nsum_{x|i}sum_{y|j}frac{yi}{x}sum_{d|x,d|y}mu(d)
]
改变枚举项为(xd,yd)，接着修改为(id,jd)：
[sum_{d=1}^Nmu(d)sum_{i=1}^Nsum_{j=1}^Nsum_{xd|i}sum_{yd|j}frac{yi}{x}\sum_{d=1}^Ndmu(d)sum_{i=1}^frac{N}{d}sum_{j=1}^frac{N}{d}sum_{x|i}sum_{y|i}frac{yi}{x}
]
分配一下：
[sum_{d=1}^Ndmu(d)(sum_{i=1}^frac{N}{d}sum_{x|i}frac{i}{x})(sum_{j=1}^frac{N}{d}sum_{y|i}y)
]
那这后面不就是约数的平方嘛。
[sum_{d=1}^Ndmu(d)sum_{i=1}^frac{N}{d}sigma^2(i)
]
显然后面可以数论分块一下，前面可以杜教筛一下。
先看看怎么样怎样杜教筛这个式子：
设(f=mu id,g=id)，那么有：
[(f*g)(n)=sum_{d|n}(mu(d)d)frac{n}{d}=nsum_{d|n}mu(d)=[n=1]
]
又有：
[g(1)S(n)=sum_{i=1}^n(f*g)(i)-sum_{i=2}^ng(i)S(frac{n}{i})
]
可以得到：
[S(n)=1-sum_{i=2}^niS(frac{n}{i})
]
之后要求这个：
[sum_{i=1}^nsigma(i)
]
怎么计算呢？
需要有点逆向思维，枚举约数不好弄就枚举倍数，所以有：
[sum_{i=1}^nsigma(i)=sum_{i=1}^nifrac{n}{i}
]
此时理清楚就可以写代码了。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9+7;
const int maxn = 3e6;
const int inv2 = 500000004;
bool vis[maxn+10];
ll mu[maxn+10];
int primes[maxn+10], cnt;
void init(int n)
{
mu[1] = 1;
for(int i = 2; i <= n; i++)
{
if(!vis[i])
{
primes[++cnt] = i;
mu[i] = -1;
}
for(int j = 1; primes[j] <= n/i; j++)
{
vis[primes[j]*i] = 1;
if(i % primes[j] == 0) break;
else mu[i*primes[j]] = -mu[i];
}
}
for(int i = 1; i <= n; i++)
{
mu[i] = i*mu[i]%mod;
mu[i] = (mu[i]+mu[i-1])%mod;
}
}
unordered_map<ll, ll> Smu;
ll getSmu(ll n)
{
if(n <= maxn) return mu[n];
if(Smu[n]) return Smu[n];
ll res = 1;
for(ll l = 2, r; l <= n; l = r+1)
{
r = n/(n/l);
res -= (r-l+1)%mod*(r+l)%mod*inv2%mod*getSmu(n/l)%mod;
res = ((res%mod)+mod)%mod;
}
return Smu[n] = res;
}
ll n;
unordered_map<ll, ll> G;
ll getG(ll n)
{
if(G[n]) return G[n];
ll res = 0;
for(ll l = 1, r; l <= n; l = r+1)
{
r = n/(n/l);
res = (res + (r-l+1)%mod*(r+l)%mod*inv2%mod*(n/l)%mod)%mod;
}
return G[n] = res;
}
int main()
{
init(maxn);
cin >> n;
ll res = 0;
for(ll l = 1, r; l <= n; l = r+1)
{
r = n/(n/l);
res = (res + ((getSmu(r)-getSmu(l-1))%mod+mod)%mod*getG(n/l)%mod*getG(n/l)%mod)%mod;
}
res = ((res%mod)+mod)%mod;
cout << res << endl;
return 0;
}