BZOJ4176 Lucas的数论 【莫比乌斯反演 + 杜教筛】

题目

去年的Lucas非常喜欢数论题，但是一年以后的Lucas却不那么喜欢了。

在整理以前的试题时，发现了这样一道题目“求Sigma(f(i)),其中1<=i<=N”，其中 表示i的约数个数。他现在长大了，题目也变难了。
求如下表达式的值：

其中 表示ij的约数个数。
他发现答案有点大，只需要输出模1000000007的值。

输入格式

第一行一个整数n。

输出格式

一行一个整数ans，表示答案模1000000007的值。

输入样例

2

输出样例

8

提示

对于100%的数据n <= 10^9。

题解

这题推导和SDOI2015约数个数和那道题是一样的
只不过计算的方式有差别

这道题没有多组询问，而且n特别大【不能O(n)实现】，要用杜教筛

最后推出式子：
\[ans=\sum\limits_{d=1}^{n} \mu(d) (\sum\limits_{i=1}^{\lfloor \frac{n}{d} \rfloor} \lfloor \frac{\lfloor \frac{n}{d} \rfloor}{i} \rfloor)^2\]

如果我们记
\[sum(n) = \sum\limits_{i=1}^{n}\lfloor \frac{n}{i} \rfloor\]
那么式子可以写成：
\[ans=\sum\limits_{d=1}^{n} \mu(d) sum(\lfloor \frac{n}{d} \rfloor)^2\]

显然可以分块计算
因为\(n<=10^9\)，所以对于\(\mu\)的前缀和我们采用杜教筛，时间复杂度\(O(n^{\frac{2}{3}}logn)\)
对于\(sum(n)\)，我们内部也分块计算，时间复杂度\(O(\int_{0}^{\sqrt{n}} x^{\frac{1}{2}} dx) = O(\frac{2}{3} n^{\frac{3}{4}}) = O(n^{\frac{3}{4}})\)

所以总的复杂度\(O(n^{\frac{2}{3}}logn + n^{\frac{3}{4}})\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#include<cstring>
#include<algorithm>
#define LL long long int
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<' '; puts("");
using namespace std;
const int maxn = 1000005,maxm = 100005,INF = 1000000000,P = 1000000007;
typedef map<LL,LL> Map;
Map _mu;
Map::iterator it;
LL p[maxn],pi,mu[maxn],N,n;
int isn[maxn];
void init(LL n){
N = (LL)pow(n,2.0 / 3.0);
mu[1] = 1;
for (int i = 2; i < N; i++){
if (!isn[i]) p[++pi] = i,mu[i] = -1;
for (int j = 1; j <= pi && i * p[j] < N; j++){
isn[i * p[j]] = true;
if (i % p[j] == 0){
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = -mu[i];
}
}
for (int i = 1; i < N; i++) mu[i] = (mu[i - 1] + mu[i]) % P;
}
LL sum(LL x){
LL ans = 0;
for (int i = 1,nxt; i <= x; i = nxt + 1){
nxt = x / (x / i);
ans = (ans + (nxt - i + 1) * (x / i) % P) % P;
}
return ans;
}
LL S(LL n){
if (n < N) return mu
;
if ((it = _mu.find(n)) != _mu.end())
return it->second;
LL ans = 1;
for (int i = 2,nxt; i <= n; i = nxt + 1){
nxt = n / (n / i);
ans = (ans - (nxt - i + 1) * S(n / i) % P) % P;
}
return _mu
 = ans;
}
int main(){
cin >> n;
init(n);
LL ans = 0;
for (int i = 1,nxt; i <= n; i = nxt + 1){
nxt = n / (n / i);
LL tmp = sum(n / i);
tmp = tmp * tmp % P;
ans = ans + (S(nxt) - S(i - 1)) % P * tmp % P;
}
ans = (ans % P + P) % P;
cout << ans << endl;
return 0;
}