这道题其实不需要杜教筛（甚至 \(O(m \log^2 m)\) 可过）。
然而我无聊写了个 \(O(m^{2/3} + \log p)\) 的解法，其中 \(p = 10^9 + 7\)。
设 \(X\) 表示数列的长度，则要求 \(E(X)\)。
有一个显然的结论是 \(E(X) = \sum\limits_{i=1}^{\infty} P(X \ge i)\)。
因为如果最后的长度为 \(x\)，则会在 \(i=1,2,\ldots,x\) 处都产生贡献，故一共有 \(x\) 次贡献，刚好为长度。
然后把 \(i=1\) 的提出来，因为序列长度一定不小于 \(1\)，所以 \(P(X \ge 1) = 1\)。
于是就要求 \(\sum\limits_{i=2}^{\infty} P(X \ge i)\)。
考虑容斥求出方案数，即用所有序列减去前 \(i - 1\) 个数互质的。
注意到只需考虑前 \(i - 1\) 个数，所以方案数也只需考虑前 \(i - 1\) 个数的，后面的方案数可以约分。
则接下来需要求有多少个长度为 \(i - 1\) 的，每个元素是 \([1,m]\) 内整数的序列满足所有元素 \(\gcd\) 为 \(1\)，设这个东西为 \(f(1)\)。
考虑 \(F(x) = \sum\limits_{x|d} f(d)\)，即所有元素的 \(\gcd\) 是 \(x\) 的倍数的序列个数，也即所有元素都是 \(x\) 的倍数的序列个数，则显然这样的序列有 \(\left\lfloor\frac mx\right\rfloor^{i-1}\) 个。
反演一下有 \(f(1) = \sum\limits_{i=1}^m \mu(i) \left\lfloor\frac mx\right\rfloor^{i-1}\)。 总结一下，就是 \[
\sum\limits_{i=2}^{\infty} P(x \ge i) = \sum\limits_{i=2}^{\infty} \left(1 - \frac1{m^{i-1}} \sum\limits_{j=1}^m \mu(j) \left\lfloor\frac mj\right\rfloor^{i-1} \right)
\]
有点不爽，把 \(i\) 移个位。 \[ \sum\limits_{i=2}^{\infty} \left(1 - \frac1{m^{i-1}} \sum\limits_{j=1}^m \mu(j) \left\lfloor\frac mj\right\rfloor^{i-1} \right) = \sum\limits_{i=1}^{\infty} \left(1 - \frac1{m^i} \sum\limits_{j=1}^m \mu(j) \left\lfloor\frac mj\right\rfloor^i \right) \]
涉及无穷级数求和，这个常数不能要，否则不好交换求和号，考虑分离出 \(j=1\) 的项与之抵消。
则 \[
\sum\limits_{i=1}^{\infty} \left(1 - \frac1{m^i} \sum\limits_{j=1}^m \mu(j) \left\lfloor\frac mj\right\rfloor^i \right) = \sum\limits_{i=1}^{\infty} \frac1{m^i} \sum\limits_{j=2}^m \mu(j) \left\lfloor\frac mj\right\rfloor^i
\]
交换一下求和号 \[ \sum\limits_{i=1}^{\infty} \frac1{m^i} \sum\limits_{j=2}^m \mu(j) \left\lfloor\frac mj\right\rfloor^i = \sum\limits_{j=2}^m \mu(j) \sum\limits_{i=1}^{\infty}\left(\frac{\left\lfloor\frac mj\right\rfloor}m\right)^i \]
然后根据一些基本数学知识（雾 \[ \sum\limits_{j=2}^m \mu(j) \sum\limits_{i=1}^{\infty}\left(\frac{\left\lfloor\frac mj\right\rfloor}m\right)^i = \sum\limits_{j=2}^m \mu(j) \frac{\left\lfloor\frac mj\right\rfloor}m \frac1{1-\frac{\left\lfloor\frac mj\right\rfloor}m} = \sum\limits_{j=2}^m \mu(j) \frac{\left\lfloor\frac mj\right\rfloor}{m - \left\lfloor\frac mj\right\rfloor} \]
然后！
当然可以线性筛素数和线性求逆元，直接算。
然而我选择数论分块 + 杜教筛 + 离线处理逆元。
复杂度是 \(O(m^{2/3} + \log p)\) 的，如果每次求逆元就是 \(O(m^{2/3} + \sqrt m\log p)\)。
代码： 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
65
66
using namespace std;
const int N = 1e5;
const int MX = 2300;
const int LIM = 320;
const int mod = 1e9 + 7;
inline int fpow(int a,int b)
{
int ret = 1;
for(;b;b >>= 1)
(b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
return ret;
}
int n;
int vis[MX + 5],cnt,prime[MX + 5],mu[MX + 5];
int tot,c[LIM * 2 + 5],a[LIM * 2 + 5],prod[LIM * 2 + 5],iprod[LIM * 2 + 5];
int smu[MX + 5];
inline int &mem(int x)
{
return smu[n / x];
}
int calc(int n)
{
if(n <= MX)
return mu[n];
if(~mem(n))
return mem(n);
int ret = 1;
for(register int l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
ret = (ret - (long long)(r - l + 1) * calc(n / l) % mod + mod) % mod;
}
return mem(n) = ret;
}
int ans = 1;
int main()
{
mu[1] = 1,memset(smu,-1,sizeof smu);
for(register int i = 2;i <= MX;++i)
{
if(!vis[i])
mu[prime[++cnt] = i] = mod - 1;
for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
mu[i * prime[j]] = (mod - mu[i]) % mod;
}
mu[i] = (mu[i] + mu[i - 1]) % mod;
}
scanf("%d",&n),prod[0] = 1;
for(register int l = 2,r;l <= n;l = r + 1)
{
r = n / (n / l);
c[++tot] = (long long)(calc(r) - calc(l - 1) + mod) * (n / l) % mod,prod[tot] = (long long)prod[tot - 1] * (a[tot] = n - n / l) % mod;
}
iprod[tot] = fpow(prod[tot],mod - 2);
for(register int i = tot;i;--i)
iprod[i - 1] = (long long)iprod[i] * a[i] % mod;
for(register int i = 1;i <= tot;++i)
ans = (ans - (long long)c[i] * iprod[i] % mod * prod[i - 1] % mod + mod) % mod;
printf("%d\n",ans);
}