JZOJ 4872 太阳神

首先取补集。
即枚举 \(i = \dfrac a{\gcd(a,b)},j = \dfrac b{\gcd(a,b)},d = \gcd(a,b)\),条件即 \(\gcd(i,j) = 1,dij \le n\)
固定了 \(i,j\) 后,\(d\) 的取值只有 \(\genfrac\lfloor\rfloor{}{}n{ij}\) 种。

然后有 \[ \newcommand\ffrac[2]{\genfrac\lfloor\rfloor{}{}{ #1 }{ #2 }} \begin{align*} & \sum\limits_{i=1}^n \sum\limits_{j=1}^n [\gcd(i,j)=1]\ffrac n{ij} \\ = & \sum\limits_{i=1}^n \sum\limits_{j=1}^n \ffrac n{ij} \sum\limits_{d|\gcd(i,j)} \mu(d) \\ = & \sum\limits_{d=1}^n \mu(d) \sum\limits_{i=1}^{\ffrac nd} \sum\limits_{j=1}^{\ffrac nd} \ffrac n{d^2ij} \end{align*} \]

修正上界,原式变为 \[ \sum\limits_{d=1}^{\lfloor\sqrt n\rfloor} \mu(d) \sum\limits_{i=1}^{\ffrac n{d^2}} \sum\limits_{j=1}^{\ffrac n{d^2i}} \ffrac n{d^2ij} \]

\(f(n) = \sum\limits_{i=1}^n \sum\limits_{j=1}^{\ffrac ni} \ffrac n{ij}\)
发现后面的和式很有趣,如果能求出 \(g(n) = \sum\limits_{i=1}^n \ffrac ni\) 就可以数论分块了。

观察发现 \(g\) 的定义,发现 \(g(n) = \sum\limits_{i=1}^n \sigma_0(i)\)
于是可以类似杜教筛,先筛出 \(n^{2/3}\) 以内的 \(g\) 值,然后在 \(n^{2/3}\) 以外的 \(g\) 值暴力数论分块做。
可以发现这一部分很像杜教筛所以复杂度也不会超过 \(O(n^{2/3})\)(人家还没递归)
其实可以令 \(\sigma_0 * \mu\) 得到 \(\textbf 1\) 这样就可以了两遍杜教筛带走了(但是常数过大)

然后对 \(f\) 数论分块,求出 \(\sum\limits_{d=1}^{\lfloor\sqrt n\rfloor} \mu(d)f(\ffrac n{d^2})\) 即可。
别忘了取补集回去。

代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <cstdio>
#include <cstring>
using namespace std;
const long long N = 1e10;
const int MX = 1e7;
const int mod = 1e9 + 7;
long long n;
int vis[MX + 5],cnt,prime[MX + 5],mu[MX + 5],d[MX + 5],mn[MX + 5],ans;
int mem[MX + 5];
inline int id(long long x)
{
return n / x;
}
int g(long long n)
{
if(n <= MX)
return d[n];
if(~mem[id(n)])
return mem[id(n)];
int ret = 0;
for(register long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret + (long long)(r - l + 1) * (n / l) % mod) % mod;
}
return mem[id(n)] = ret;
}
int f(long long n)
{
int ret = 0;
for(register long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret + (long long)(r - l + 1) * g(n / l) % mod) % mod;
}
return ret;
}
int main()
{
freopen("ra.in","r",stdin),freopen("ra.out","w",stdout);
memset(mem,-1,sizeof mem),mu[1] = d[1] = 1;
for(register int i = 2;i <= MX;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = mod - 1,d[i] = mn[i] = 2;
for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
{
d[i * prime[j]] = d[i] / mn[i] * (mn[i] + 1),mn[i * prime[j]] = mn[i] + 1;
break;
}
else
d[i * prime[j]] = d[i] * 2,mn[i * prime[j]] = 2,mu[i * prime[j]] = mod - mu[i];
}
d[i] = (d[i - 1] + d[i]) % mod;
}
scanf("%lld",&n);
for(register int i = 1;(long long)i * i <= n;++i)
ans = (ans + (long long)mu[i] * f(n / i / i) % mod) % mod;
ans = ((long long)(n % mod) * (n % mod) % mod - ans + mod) % mod;
printf("%d\n",ans);
}