算法数学

发布日期: 2022-11-22

更新日期: 2022-11-22

文章字数: 2k

数论分块学习笔记

设函数 $f,g : \mathbb{N} \to \mathbb{N}$ ，并且

前缀和 $S_r-S_{l-1}$ 可以在 $\mathcal{O}(1)$ 计算，其中 $S_i = \sum_{j=1}^{i} f(j),~S_0 = 0$ 。
$g(i)$ 形如 $\left\lfloor\dfrac{n}{i}\right\rfloor$ ，有时候也可以形如 $\left\lfloor\dfrac{m}{ki+b}\right\rfloor$ ，总之取值具有阶梯状分布的特点。

数论分块可以在 $\mathcal{O}(\sqrt{n})$ 的复杂度内计算形如 $\sum_{i=1}^n f(i) g(i)$ 的和式。

数学原理

证明详见参考文献[1]，以及这里一段可以直接跳过，见算法流程

引理1：

$\forall a, b, c \in \mathbb{Z},~\left\lfloor\frac{a}{b c}\right\rfloor=\left\lfloor\frac{\left\lfloor\frac{a}{b}\right\rfloor}{c}\right\rfloor$

引理2：

$\forall n \in \mathbb{N}_{+},~\left|\textstyle \left\{\left\lfloor\frac{n}{i}\right\rfloor \mid i \leq n\right\}\right| \leq\lfloor 2 \sqrt{n}\rfloor$

数论分块结论：

对于常数 $n$，使得 $\left\lfloor\frac ni\right\rfloor=\left\lfloor\frac nj\right\rfloor$ 成立的最大的满足 $i\leq j\leq n$ 的 $j$ 的值为 $\left\lfloor\dfrac n{\lfloor\frac ni\rfloor}\right\rfloor$。

即值 $\left\lfloor\dfrac ni\right\rfloor$ 所在的块的右端点为 $\left\lfloor\dfrac n{\lfloor\frac ni\rfloor}\right\rfloor$ ，并且根据引理2可知块数不超过 $\mathcal{O}(\sqrt{n})$ 。

算法讲解

我并不想像其他人的文章一样直接抛出算法流程，我更想用一种通俗解释来描述它。

刚刚的引理1就是个辅助的结论，引理2和主要结论是最重要的。

观察 $g(i) = \left\lfloor\frac{n}{ki}\right\rfloor$ 的图像，可以发现它呈阶梯状分布。

举个例子，$g(x) = \left\lfloor\frac{10}{x}\right\rfloor$ 的图像如下图：（负半轴不合法，无视即可）

显然它呈阶梯状分布。而数论分块做的就是快速统计每个块并计算贡献。

而块的数量，由引理2可知为 $\mathcal{O}(\sqrt{n})$ ，并且主要结论做到了快速统计，因此复杂度是 $\mathcal{O}(\sqrt{n})$ 的。

算法流程

算法流程：

$\begin{array}{ll} 1 & \texttt{获取}\,f(i)\,\texttt{函数的前缀和，记为}\,s(i)\,\texttt{。}\\ 2 & l\gets 1\\ 3 & r\gets 0\\ 4 & \textsf{result}\gets 0\\ 5 & \textbf{while }l\leq n\textbf{ do}:\\ 6 & \qquad r\gets\left\lfloor\dfrac n{\lfloor\frac nl\rfloor}\right\rfloor\\ 7 & \qquad \textsf{result}\gets (S_r-S_{l-1})\times\left\lfloor\dfrac nl\right\rfloor\\ 8 & \qquad l\gets r+1\\ 9 & \textbf{end while}\\ \end{array}$

代码实现

模板题：洛谷P3935 Calculating

答案要求计算形如 $F(r) - F(l-1)$ 的式子（具体推导详见洛谷P3935 Calculating 题解）

$F(n) = \sum_{i=1}^n\left\lfloor\frac{n}{i}\right\rfloor$

首先这里的 $S_i = i$ ，因为 $f(i) = 1$ 。

模板代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
const int mod = 998244353;
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
void add(int &x,int y) { (x += y) >= mod ? x -= mod : 0; }
#define N ((int)())

int solve(int n)
{
    int res = 0;
    for(int l = 1, r = 0; l <= n; l = r + 1) // 注意这里不是 l++ 
    {
        r = n / (n / l); // 右端点
        // 增量为 [S(r) - S(l-1)] * (n / l) ， 后面的是 n / l ，这是算法规定的。
        add(res, (r - l + 1) * (n / l) % mod);
    }
    return res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int L,R; cin >> L >> R;
    cout << (solve(R) - solve(L - 1) + mod) % mod << endl;
    return 0;
}

例题

方便起见，这里指讲怎么算最后的答案，推导详见对应的题解。

例题1：来自洛谷P2261 [CQOI2007]余数求和题解

计算下列式子：

$n k-\sum_{i=1}^n i\left\lfloor\frac{k}{i}\right\rfloor$

数据范围：$1 \le n,k \le 10^9$ 。

分析：这个就很板子，主要是这个前缀和。

根据等差数列的性质不难发现 $S_r - S_{l-1}$ 就是 $(r-l + 1) \times \dfrac{l + r}{2}$ 。

注意当 $k \le n$ 时只要枚举到 $k$ ，因为后面 $\left\lfloor\frac{k}{i}\right\rfloor$ 全为 $0$ ，所以不需要枚举（而且枚举后代码会出错）

时间复杂度 $\mathcal{O}(\sqrt{n})$

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
int n,k;
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n >> k;int res=n*k;
    for(int l=1,r; l<=min(k,n); l=r+1)
    {
        r=min(k/(k/l),n);
        res-=(k/l)*(r-l+1)*(l+r)/2;
    }
    cout << res << endl;
    return 0;
}

例题2：来自 CF1485C Floor and Mod 题解

多组询问，计算下列式子

$\sum_{b=1}^y\left\lfloor\frac{x}{b+1}\right\rfloor$

数据范围：$1\le Q \le 100,~1 \leq x, y \leq 10^9$ 。

分析：

遗憾的是，数论分块不能直接算这样的 $b+1$ 。

咋办，换元（改写以下式子）呗

$\sum\limits_{b=2}^{y+1}\left\lfloor\dfrac{x} {b}\right\rfloor$

好啦，没啦。只要令 $l=2$ 跑就好啦。

不过这道题是分类讨论（~~有点根号分治的味道~~），所以建议直接查看题解。

时间复杂度 $\mathcal{O}(Q\sqrt{x})$

#define int long long
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int Q,a,b,x,y,res;
    cin >> Q;
    while(Q--)
    {
        cin >> x >> y;
        a=1,b=1,res=0;
        for(; b*b+b-1<x&&b<=y; b++) // 这里不用管，这是原题分类讨论的地方
            res+=(b*b+b-1)/(b+1);
        for(int l=b+1,r; l<=min(x,y+1); l=r+1) // 其实主要就是改了下式子（这里代码可读性很差）
        {
            r=min(x/(x/l),y+1);
            res+=x/l*(r-l+1); // 你会发现这里没啥区别
        }
        cout << res << '\n';
    }
    return 0;
}

相信你已经对数论分块有初步的了解了

例题3：来自洛谷P3327 [SDOI2015]约数个数和题解

计算下列式子

$\sum_{d=1}^n \mu(d) \sum_{x=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\left\lfloor\frac{n}{d x}\right\rfloor \sum_{y=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\left\lfloor\frac{m}{d y}\right\rfloor$

数据范围：多组数据，$1\le T,n,m \le 50000$ 。

分析：

这个式子和板子差距有亿点大，但是写法还是类似的

注意到这里的 $dx,dy$ 其实你把它的 $d$ 放上去就是可以算的 $\left\lfloor\frac{\left\lfloor\frac{n}{d}\right\rfloor}{x}\right\rfloor$ 了

可以预处理一下这个东西为 $g\left(\left\lfloor\frac{n}{d}\right\rfloor\right)$ ，然后每组询问查表就好了。

~~感觉这么说比较抽象，直接看详细题解可能会比较好（大雾~~

接下来假设已经知道 $g$ 的所有值了，然后怎么做呢？

对于前面的那个求和，因为后面的那个 $g$ 也构成了阶梯状的值，很好，我们继续数论分块。

时间复杂度 $\mathcal{O}(T\sqrt{n} + n\sqrt{n})$

代码：（早期代码+注释+删掉了一些不影响阅读的东西）

#define int long long
#define N (int)(5e4+15)
int prime[N],sum[N],pcnt,mu[N],g[N];
bool ck[N];
int solve(int n)
{
    int res=0;
    for(int l=1,r; l<=n; l=r+1)
    {
        r=(n/(n/l));
        res+=(r-l+1)*(n/l);
    }
    return res;
}
void Mobius() // 筛法求 μ(n)
{
    mu[1]=1;
    for(int i=2; i<=N-5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=N-5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;
            if(i%prime[j])
            {
                mu[pos]=-mu[i];
            }else
            {
                mu[pos]=0;
                break;
            }
        }
    }
    for(int i=1; i<=N-5; i++) sum[i]+=sum[i-1]+mu[i]; // 做了一个 μ 的前缀和
    for(int i=1; i<=N-5; i++) g[i]=solve(i); // 预处理了模板题的那个式子
}
int Q,n,m;
signed main()
{
    Mobius();
    read(Q);
    while(Q--)
    {
        read(n);read(m);
        if(n>m)swap(n,m);
        int res=0;
        for(int l=1,r; l<=n; l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            res+=(sum[r]-sum[l-1])*g[n/l]*g[m/l];
        }
        write(res);pc('\n');
    }
    return 0;
}