数论分块入门

七海

发布于：2021年7月8日

《退役人的自我救赎系列》——其一，也就是基础的数论知识。

整除分块

求解 \(\sum_{i=1}^n \frac{n}i\) 的值，\(n\) 是 \(10^{10}\) 级别。

暴力是 \(\mathcal{O}(n)\) 的，死了；但是 \(\lfloor \frac{n}i\rfloor\) 的结果是根据块状分布的，且最多只会有 \(2\sqrt{n}\) 种不同的值。因此，我们可以利用这个性质得到一个 \(\mathcal{O}(\sqrt{n})\) 的算法。

直接地说，\(\lfloor \frac{n}i\rfloor\) 的值在一段连续的区间内具有相同的值，且这个区间具有右端点 \(\lfloor \frac{n}{\lfloor\frac{n}i\rfloor}\rfloor\)。

for (int l = 1, r; l <= n; l = r + 1) {
    r = min(n / (n / l), n);
    procedure(l, r);	// [l, r] 的值是 n / l
}

上述的写法需要注意的就是 n / l 可能为 0 导致除法除零，一般来说需要特判一下。

CQOI2007：余数求和

那来点模板题。求的是模数和，看起来和上面的整除分块没有什么关系，但是： \[ \sum_{i=1}^n k\mod i = \sum_{i=1}^n k - \lfloor\frac{k}i\rfloor i = nk - \sum_{i=1}^n \lfloor\frac{k}i\rfloor i \] 这样处理之后，右边的部分仍然不是我们熟悉的整出分块形式，所以还需要进一步处理；

我们已经知道了 \(\sum_{i=1}^n \lfloor\frac{k}i\rfloor\) 在一定范围内具有相同的值，那么在块内我们令 \(T_j = \lfloor\frac{k}j\rfloor\)，那么就可以进行如下形式的化简: \[ \sum_{i=1}^n \lfloor\frac{k}i\rfloor i = \sum_{i=1}^n T_i\cdot i = \sum_j\sum_{i=L}^R T_i\cdot \sum_{i=L}^R i \] 这样，在每一个分块内，左边是整除分块，右边是一个公差为 1 的等差数列；块内只需要一个定值和一个很好求的和，就可以分块来做了：

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    auto n = read(), k = read();
    auto ans = n * k;
    for (llong l = 1, r; l <= n; l = r + 1) {
        if (k / l) r = min(k / (k / l), n);
        else r = n;
        ans -= (r - l + 1) * (l + r) / 2 * (k / l);
    }
    cout << ans << endl;
    return 0;
}

好了，现在你已经学会了数论分块了！

清华集训2012：模积和

题目是要求了 \(i\ne j\)，所以我们需要转化为下面的形式： \[ \sum_{i=1}^n\sum_{j=1}^m[i\ne j](n\ \text{mod}\ i)(m\ \text{mod}\ j) \\ = \sum_{i=1}^n\sum_{j=1}^m(n\ \text{mod}\ i)(m\ \text{mod}\ j) - \sum_{i=1}^{\min(n, m)}(n\ \text{mod}\ i)(m\ \text{mod}\ i) \] 然后，先处理上述式子的第一项；最直接的做法就是拆开化简： \[ {align} \begin{align} 原式 &= \sum_{i=1}^n\sum_{j=1}^m (n-\lfloor\frac{n}i\rfloor i)(m - \lfloor\frac{m}j\rfloor j) \\ &= \sum_{i=1}^n\sum_{j=1}^m nm - \lfloor\frac{n}i\rfloor im - \lfloor\frac{m}j\rfloor jn + \lfloor\frac{m}j\rfloor\lfloor\frac{n}i\rfloor i j \\ &= n^2m^2 - nm^2\sum_{i=1}^n\lfloor\frac{n}i\rfloor i-n^2m\sum_{j=1}^m\lfloor\frac{m}j\rfloor j + nm\sum_{i=1}^n\lfloor\frac{n}i\rfloor i\cdot\sum_{j=1}^m\lfloor\frac{m}j\rfloor j \end{align} \] 当然，也可以类似于上面求余数和的做法那样直接拆成： \[ \sum_{i=1}^n(n\ \text{mod}\ i) \cdot \sum_{j=1}^m(m\ \text{mod}\ j) \] 也可以拆出相同的结果。

总而言之，上面的式子就被拆成了一些和上面模数求和一样的形式的组合，只需要分别求出两个互不相关的部分的和就可以算出题目种式子的左边部分；接下来处理右边部分，我们先约定 \(k = \min(n, m)\)： \[ \begin{align} 原式 &= \sum_{i=1}^k (n-\lfloor\frac{n}i\rfloor i)(m - \lfloor\frac{m}i\rfloor i) \\ &= \sum_{i=1}^k(nm - m\lfloor\frac{n}i\rfloor i - n\lfloor\frac{m}i\rfloor i + \lfloor\frac{m}i\rfloor\lfloor\frac{n}i\rfloor i^2) \\ &= knm - m\sum_{i=1}^k\lfloor\frac{n}i\rfloor i - n\sum_{i=1}^k\lfloor\frac{m}i\rfloor i + \sum_{i=1}^k\lfloor\frac{n}i\rfloor \lfloor\frac{m}i\rfloor i^2 \end{align} \] 还是一样可以使用模数分块的方法求出。

此外，提一下考研常用公式：\(\sum_{i=1}^n i^2 = \frac16n(n+1)(2n+1)\)

constexpr lll mod = 19940417;

lll fast_pow(lll a, lll b = mod - 2) {
    lll ret = 1;
    while (b) {
        if (b & 1) (ret *= a) %= mod;
        (a *= a) %= mod, b >>= 1;
    }
    return ret %= mod;
}

llong euler_phi(llong n) {
    auto ans = n;
    for (llong i = 2; i * i <= n; ++i)
        if (!(n % i)) {
            ans = ans / i * (i - 1);
            while (!(n % i)) n /= i;
        }
    if (n > 1) ans = ans / n * (n - 1);
    return ans;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    lll n = read(), m = read(), k = min(n, m);
    auto nm = n * m % mod, n2m = nm * n % mod;
    auto n2m2 = nm * nm % mod, nm2 = nm * m % mod;
    const auto calc = [](lll n, lll k) {
        lll ret = 0;
        minimize(n, k);
        for (lll l = 1, r; l <= n; l = r + 1) {
            if (k / l) r = min(k / (k / l), n);
            else r = n;
            auto sum = (r - l + 1) * (l + r) / 2 % mod;
            (ret += sum * (k / l) % mod) %= mod;
        }
        return ret %= mod;
    };
    const auto sum_mod = [&calc](lll n, lll k) {
        return (n * k % mod - calc(n, k) + mod) % mod;
    };
    const auto phi_mod = euler_phi(mod) - 1;
    const auto inv6 = fast_pow(6, phi_mod);
    const auto sum_i2 = [inv6](lll n) {
        return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
    };
    const auto lp = sum_mod(n, n) * sum_mod(m, m) % mod;
    auto knm = k * nm % mod, rp = knm;
    for (lll l = 1, r; l <= k; l = r + 1) {
        if (n / l && m / l)
            r = min({n / (n / l), m / (m / l), k});
        else r = k;
        auto si2 = (sum_i2(r) - sum_i2(l - 1) + mod) % mod;
        auto inm = (n / l) * (m / l) % mod;
        (rp += inm * si2 % mod) %= mod;
    }
    auto ckn = calc(k, n), ckm = calc(k, m);
    auto tmp = (m * ckn % mod + n * ckm % mod) % mod;
    rp = (rp + mod - tmp) % mod;
    cout << llong((lp + mod - rp) % mod) << endl;
    return 0;
}

需要注意的是，\(19940417 = 7\times2848631\)，并不是一个质数，所以不能用费马小定理求出模它意义下的逆元；但是 \(6\) 和它互质，故应当使用欧拉定理：

欧拉定理：若 \(\gcd(a, m) = 1\)，则 \(a^{\phi(m)} \equiv 1 \ (\text{mod} \ m)\)。
欧拉函数：\(\phi(n)\)，表示了小于等于 \(n\) 的正整数和 \(n\) 互质的整数的个数。
标准分解式：将质因数分解的结果按照大小，由小到大排列，并将相同质因数连乘以指数形式表示。
如果一个数字的标准分解式可以写成：\(n = p_1^{k_1}p_2^{k_2}\cdots p_r^{k_r}\)，那么： \[ \phi(n) = \prod_{i=1}^rp_i^{k_i-1}(p_i-1)=\prod_{p|n}p^{\alpha_p-1}(p-1)=n\prod_{p|n}(1-\frac1p) \] 可以使用上述公式计算欧拉函数。

求解单个欧拉函数值，可以使用 Pollard Rho 算法优化后根据定义求解；多个欧拉函数值则可以使用线性筛求解。

实战（？）

给了 \(x_1, x_2, y_1, y_2\)，求解： \[ \sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor+\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor)^2 \] 结果对 \(10^9+7\) 取模。

这是 CCPC 湘潭邀请赛 2021 的 C 题，无处补题。

首先还是直接把它们乘开： \[ \begin{align} 原式 =& (y_2-y_1+1)\sum_{i=x_1}^{x_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)^2 + (x_2-x_1+1)\sum_{j=y_1}^{y_2}(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor)^2 + \\ &2\sum_{i=x_1}^{x_2}\sum_{j=y_1}^{y_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor) \end{align} \] 可以看出前两项是一类的，第三项是另一类的；其中，前两项又可以进行这样的拆分： \[ 原式 = \sum_{i=x_1}^{x_2}\lfloor\frac{i}{x_1}\rfloor^2 + 2\sum_{i=x_1}^{x_2}\lfloor\frac{i}{x_1}\rfloor\lfloor\frac{x_2}i\rfloor + \sum_{i=x_1}^{x_2}\lfloor\frac{x_2}i\rfloor^2 \] 又可以分成三项；其中第三项是我们这里提到的整除分块问题，第一项是我们再熟悉不过的“分块”，可以直接计算；那么问题就变成了中间的一部分——在这里我想要说的是这并没有什么特殊的性质（或者是我没有发现），只需要分块再分块然后直接计算就可以了。

那么回到第一组式子的第三项；注意到 \(i\) 和 \(j\) 是互不相关的，所以直接分别计算然后乘起来就行： \[ 原式 = \sum_{i=x_1}^{x_2}(\lfloor\frac{i}{x_1}\rfloor+\lfloor\frac{x_2}i\rfloor)\cdot\sum_{j=y_1}^{y_2}(\lfloor\frac{j}{y_1}\rfloor+\lfloor\frac{y_2}j\rfloor) \] 只需要分别计算后求和，之后再在模意义下乘起来就得到了。

constexpr llong mod = 1e9 + 7;

constexpr lll fast_pow(lll a, lll b = mod - 2) {
    lll ret = 1;
    while (b) {
        if (b & 1) (ret *= a) %= mod;
        (a *= a) %= mod, b >>= 1;
    }
    return ret %= mod;
}

constexpr lll inv6 = fast_pow(6), inv2 = fast_pow(2);

lll sum_i2(lll n) {
    return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}

lll sum_i(lll n) {
    return (n + 1) * n % mod * inv2 % mod;
}

lll xi2(lll l, lll n, lll k) {
    lll ret = 0;
    minimize(n, k);
    for (lll r; l <= n; l = r + 1) {
        if (k / l) r = min(k / (k / l), n);
        else r = n;
        auto val = (k / l) * (k / l) % mod;
        (ret += (r - l + 1) * val % mod) %= mod;
    }
    return ret;
}

lll xi(lll l, lll n, lll k) {
    lll ret = 0;
    minimize(n, k);
    for (lll r; l <= n; l = r + 1) {
        if (k / l) r = min(k / (k / l), n);
        else r = n;
        (ret += (r - l + 1) * (k / l) % mod) %= mod;
    }
    return ret;
}

lll sum_ix(lll n, lll x) {
    auto m = n % x + 1, d = n / x;
    if (!d) return 0;
    lll ret = x * sum_i(d - 1) % mod;
    (ret += m * d % mod) %= mod;
    return ret;
}

lll sum_ix2(lll n, lll x) {
    auto m = n % x + 1, d = n / x;
    if (!d) return 0;
    lll ret = x * sum_i2(d - 1) % mod;
    auto d2 = d * d % mod;
    (ret += m * d2 % mod) %= mod;
    return ret;
}

lll ix(lll l, lll r, lll k) {
    return (sum_ix(r, k) + mod - sum_ix(l - 1, k)) % mod;
}

lll ix2(lll l, lll r, lll k) {
    return (sum_ix2(r, k) + mod - sum_ix2(l - 1, k)) % mod;
}

lll ix_xi(lll l, lll n, lll x1, lll x2) {
    lll ret = 0;
    minimize(n, x2);
    for (lll r; l <= n; l = r + 1) {
        if (x2 / l) r = min(x2 / (x2 / l), n);
        else r = n;
        auto tmp = ix(l, r, x1);
        (ret += tmp * (x2 / l) % mod) %= mod;
    }
    return ret;
}

signed main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
#endif
    for (auto T = read(); T --;) {
        auto x1 = read(), x2 = read(), y1 = read(), y2 = read();
        auto dy = y2 - y1 + 1, dx = x2 - x1 + 1;
        auto part1 = ix2(x1, x2, x1) + xi2(x1, x2, x2) + ix_xi(x1, x2, x1, x2) * 2;
        auto part2 = ix2(y1, y2, y1) + xi2(y1, y2, y2) + ix_xi(y1, y2, y1, y2) * 2;
        auto part3L = (ix(x1, x2, x1) + xi(x1, x2, x2)) % mod;
        auto part3R = (ix(y1, y2, y1) + xi(y1, y2, y2)) % mod;
        auto part3 = part3L * part3R % mod * 2 % mod;
        auto part12 = (part1 %= mod) * dy % mod + (part2 %= mod) * dx % mod;
        auto ans = (part12 + part3) % mod;
        cout << (llong)ans << endl;
    }
    return 0;
}

已通过五组测试样例。

后记

参考资料

更新于：2021年7月11日

算法

ACM

数论

NTT 简单入门

《退役人的自我救赎系列》——其二，也就是差不多算是初等数论的知识。一句话简介：NTT 即快速数论变换，是一种可以在 \(n\log n\) 的时间内完成多项式乘法的算法——的一部分。前置知识...

计算机网络实验报告（三合一）

本报告涉及到的源代码的位置：https://github.com/ma-hunter/cn_exp 班级姓名学号软件工程 18XX 班 XXX U2018XXXXX ...