Codeforces 1139D Steps to One

这道题其实不需要杜教筛（甚至 O(mlog2m) 可过）。
然而我无聊写了个 O(m2/3+logp) 的解法，其中 p=109+7。

设 X 表示数列的长度，则要求 E(X)。
有一个显然的结论是 E(X)=∞∑i=1P(X≥i)。
因为如果最后的长度为 x，则会在 i=1,2,…,x 处都产生贡献，故一共有 x 次贡献，刚好为长度。

然后把 i=1 的提出来，因为序列长度一定不小于 1，所以 P(X≥1)=1。
于是就要求 ∞∑i=2P(X≥i)。

考虑容斥求出方案数，即用所有序列减去前 i−1 个数互质的。
注意到只需考虑前 i−1 个数，所以方案数也只需考虑前 i−1 个数的，后面的方案数可以约分。
则接下来需要求有多少个长度为 i−1 的，每个元素是 [1,m] 内整数的序列满足所有元素 gcd 为 1，设这个东西为 f(1)。
考虑 F(x)=∑x|df(d)，即所有元素的 gcd 是 x 的倍数的序列个数，也即所有元素都是 x 的倍数的序列个数，则显然这样的序列有 ⌊mx⌋i−1 个。
反演一下有 f(1)=m∑i=1μ(i)⌊mx⌋i−1。 总结一下，就是 ∞∑i=2P(x≥i)=∞∑i=2(1−1mi−1m∑j=1μ(j)⌊mj⌋i−1)

有点不爽，把 i 移个位。 ∞∑i=2(1−1mi−1m∑j=1μ(j)⌊mj⌋i−1)=∞∑i=1(1−1mim∑j=1μ(j)⌊mj⌋i)

涉及无穷级数求和，这个常数不能要，否则不好交换求和号，考虑分离出 j=1 的项与之抵消。
则 ∞∑i=1(1−1mim∑j=1μ(j)⌊mj⌋i)=∞∑i=11mim∑j=2μ(j)⌊mj⌋i

交换一下求和号 ∞∑i=11mim∑j=2μ(j)⌊mj⌋i=m∑j=2μ(j)∞∑i=1(⌊mj⌋m)i

然后根据一些基本数学知识（雾 m∑j=2μ(j)∞∑i=1(⌊mj⌋m)i=m∑j=2μ(j)⌊mj⌋m11−⌊mj⌋m=m∑j=2μ(j)⌊mj⌋m−⌊mj⌋

然后！
当然可以线性筛素数和线性求逆元，直接算。
然而我选择数论分块 + 杜教筛 + 离线处理逆元。
复杂度是 O(m2/3+logp) 的，如果每次求逆元就是 O(m2/3+√mlogp)。

代码：

#include <cstdio>
#include <cstring>
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);
}