「算法笔记」类欧几里德算法

类欧几里德算法基于欧几里德算法,本质是通过取模缩小问题规模,递归求解。


为了简化代码,定义以下函数,下文不再重复:

$$ \begin{aligned} S0(n) &= \sum_{i = 0} ^ n i ^ 0 = n + 1 \\ S1(n) &= \sum_{i = 0} ^ n i ^ 1 = \frac{n(n + 1)}{2} \\ S2(n) &= \sum_{i = 0} ^ n i ^ 2 = \frac{n(n+1)(2n+1)}{6} \end{aligned} $$

1. $f(a,b,c,n)$

定义

$$ f(a,b,c,n)=\sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor $$

推导

  1. 当 $a\ge c$ 或 $b\ge c​$ 时:

$$ \begin{align*} f(a,b,c,n) & = \sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor \\ & = \sum_{i=0}^n \left(i\left\lfloor\frac{a}{c}\right\rfloor+\left\lfloor\frac{b}{c}\right\rfloor\right)+\sum_{i=0}^n \left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor\tag 1 \\ & = \frac{n(n+1)}{2}\left\lfloor\frac{a}{c}\right\rfloor+(n+1)\left\lfloor\frac{b}{c}\right\rfloor + f(a\bmod c,b\bmod c,c,n)\tag 2 \end{align*} $$

式子 $(1)​$:我们把 $a​$ 拆成 $\left\lfloor\frac{a}{c}\right\rfloor\cdot c​$ 和 $a\bmod c​$ 两部分,$b​$ 同理。就可以将原式拆成两项了。

式子 $(2)$:第一项直接 $\mathcal O(1)$ 计算,第二项递归计算

  1. 当 $a < c​$ 且 $b < c​$ 时:

我们将 $\frac{ai+b}{c}​$ 看作是一条线段 $y=\frac{ax+b}{c}(x\in [0,n])​$,那么原式的本质就是求线段下方、$x​$ 轴上方有多少个整点(包括线段,不包括 $x​$ 轴)。

设 $m=\left\lfloor\frac{an+b}{c}\right\rfloor$,枚举所有点并判断是否在线段下方。

$$ \begin{align*} f(a,b,c,n) & = \sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor \\ & = \sum_{i=0}^n \sum_{j=1}^m \left[j\le\frac{ai+b}{c}\right]\tag 1 \\ & = \sum_{i=0}^n \sum_{j=1}^m \left[ai > cj-b-1\right]\tag 2 \\ & = \sum_{j=1}^m \sum_{i=0}^n \left[i\ge \left\lfloor\frac{cj-b-1}{a}\right\rfloor+1\right]\tag 3 \\ & = \sum_{j=1}^m \left(n-\left\lfloor\frac{cj-b-1}{a}\right\rfloor\right)\tag 4 \\ & = nm-\sum_{j=0}^{m-1} \left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\tag 5 \\ & = nm-f(c,c-b-1,a,m-1)\tag 6 \end{align*} $$

式子 $(1)​$:转化为线段下方的整点个数。

式子 $(2)​$:移项后再不等式右侧减去 $1​$ 使得 $\ge​$ 变为 $>​$ 号。

式子 $(3)$:将 $>$ 一个实数转化为 $\ge$ 一个整数。我们发现如果式子 $(2)$ 中不进行转化,那么就会出现 $i\ge \frac{cj-b}{a}$,这个式子是不能直接上取整的(比如 $\lceil3\rceil=4​$ 显然是一个反例)。

式子 $(4)​$:我们把 $\sum_{i=1}^n​$ 给去掉,直接变为一个整数。

式子 $(5)$:将 $\sum$ 给拆开,并改变求和下界和上界

式子 $(6)$:直接将 $\sum$ 转化为 $f$ 函数,递归计算

  1. 边界条件 $a=0$:答案为 $(n+1)\left\lfloor\frac{b}{c}\right\rfloor$。

伪代码

int f(int a, int b, int c, int n) {
    if (!a) {
        return b / c * S0(n);
    }
    if (a >= c || b >= c) {
        return a / c * S1(n) + b / c * S0(n) + f(a % c, b % c, c, n);
    }
    int m = (a * n + b) / c;
    return n * m - f(c, c - b - 1, a, m - 1);
}

2. $g(a,b,c,n)$

定义

$$ g(a,b,c,n)=\sum_{i=0}^n i\left\lfloor\frac{ai+b}{c}\right\rfloor $$

推导

  1. 当 $a\ge c​$ 或 $b\ge c​$ 时:

$$ \begin{align*} g(a,b,c,n) & =\sum_{i=0}^n i\left\lfloor\frac{ai+b}{c}\right\rfloor \\ & =\sum_{i=0}^n \left(i^2\left\lfloor\frac{a}{c}\right\rfloor+i\left\lfloor\frac{b}{c}\right\rfloor\right)+\sum_{i=0}^n i\left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor\tag 1 \\ & =\frac{n(n+1)(2n+1)}{6}\left\lfloor\frac{a}{c}\right\rfloor+\frac{n(n+1)}{2}\left\lfloor\frac{b}{c}\right\rfloor+g(a\bmod c,b\bmod c,c,n)\tag 2 \end{align*} $$

式子 $(1)$:我们还是仿照 $f$ 的过程拆成两项然后分别求和。

式子 $(2)$:根据平方和公式 $\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}$,我们可以 $\mathcal O(1)$ 计算出第一项的值,第二项递归计算

  1. 当 $a < c​$ 且 $b < c​$ 时:

总体思路还是和 $f​$ 函数一样,利用其几何意义来求解,设 $m=\left\lfloor\frac{an+b}{c}\right\rfloor​$。

$$ \begin{align*} g(a,b,c,n) & = \sum_{i=0}^n i\left\lfloor\frac{ai+b}{c}\right\rfloor \\ & = \sum_{i=0}^n i\times \sum_{j=1}^m \left[j\le\frac{ai+b}{c}\right]\tag 1 \\ & = \sum_{i=0}^n i\times \sum_{j=1}^m \left[ai > cj-b-1\right]\tag 2 \\ & = \sum_{j=1}^m \sum_{i=0}^n i\left[i\ge \left\lfloor\frac{cj-b-1}{a}\right\rfloor+1\right]\tag 3 \\ & = \sum_{j=1}^m \frac{1}{2}\left(n+\left\lfloor\frac{cj-b-1}{a}\right\rfloor+1\right)\left(n-\left\lfloor\frac{cj-b-1}{a}\right\rfloor\right)\tag 4 \\ & = \frac{mn(n+1)}{2}-\frac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor-\frac{1}{2}\sum_{j=0}^{m-1}\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor^2 \tag 5 \\ & = \frac{mn(n+1)}{2}-\frac{1}{2}f(c,c-b-1,a,m-1)-\frac{1}{2}h(c,c-b-1,a,m-1) \tag 6 \end{align*} $$

式子 $(1)$:转化为线段下方的整点个数。

式子 $(2)$:移项后再不等式右侧减去 $1$ 使得 $\ge$ 变为 $>$ 号。

式子 $(3)$:依据还是和 $f$ 函数相同。

式子 $(4)​$:由于第二个 $\sum​$ 内就是对 $i\in \left[\left\lfloor\frac{cj-b-1}{a}\right\rfloor+1,n\right]​$ 求和,因此直接利用等差数列求和

式子 $(5)$:将等差数列求和公式继续展开,转化为 $3​$ 项。

式子 $(6)​$:第一项可以 $\mathcal O(1)​$ 计算,第二项和第三项发现是 $f​$ 和 $h​$ 函数,可以递归计算

  1. 边界条件 $a=0​$:答案为 $\frac{n(n+1)}{2}\left\lfloor\frac{b}{c}\right\rfloor​$。

伪代码

int g(int a, int b, int c, int n) {
    if (!a) {
        return b / c * S1(n);
    }
    if (a >= c || b >= c) {
        return a / c * S2(n) + b / c * S1(n) + g(a % c, b % c, c, n);
    }
    int m = (a * n + b) / c;
    return m * S1(n) - f(c, c - b - 1, a, m - 1) / 2 - h(c, c - b - 1, a, m - 1) / 2;
}

3. $h(a,b,c,n)$

定义

$$ h(a,b,c,n)=\sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor^2 $$

推导

  1. 当 $a\ge c​$ 或 $b\ge c​$ 时:

$$ \begin{align*} h(a,b,c,n) & =\sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor^2 \\ & = \sum_{i=0}^n \left(i\left\lfloor\frac{a}{c}\right\rfloor+\left\lfloor\frac{b}{c}\right\rfloor+\left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor\right)^2 \tag 1 \\ & = \sum_{i=0}^n i^2\left\lfloor\frac{a}{c}\right\rfloor^2+\left\lfloor\frac{b}{c}\right\rfloor^2+\left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor^2+2i\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor+2i\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor+2\left\lfloor\frac{b}{c}\right\rfloor\left\lfloor\frac{a\bmod c\times i+b\bmod c}{c}\right\rfloor \tag 2\\ & = \frac{n(n+1)(2n+1)}{6}\left\lfloor\frac{a}{c}\right\rfloor^2+(n+1)\left\lfloor\frac{b}{c}\right\rfloor^2+n(n+1)\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor+h(a\bmod c,b\bmod c,c,n)+2\left\lfloor\frac{a}{c}\right\rfloor g(a\bmod c,b\bmod c,c,n)+2\left\lfloor\frac{b}{c}\right\rfloor f(a\bmod c,b\bmod c,c,n) \tag 3 \end{align*} $$

式子 $(1)$:把原式的 $a$ 和 $b$ 均拆成两个部分。

式子 $(2)​$:暴力展开平方项。

式子 $(3)$:将展开式中能 $\mathcal O(1)$ 直接计算的提出,剩下的转化为 $f,g,h​$ 函数递归计算

  1. 当 $a < c$ 且 $b < c$ 时:

我们试图把平方项用另一种形式表示出来:$n^2=2\sum_{i=1}^n i-n$。我们就根据这个公式进行转化(其中 $m$ 的值依旧是 $\left\lfloor\frac{an+b}{c}\right\rfloor​$。

$$ \begin{align*} h(a,b,c,n) & =\sum_{i=0}^n \left\lfloor\frac{ai+b}{c}\right\rfloor^2 \\ & = \sum_{i=0}^n \left(2\times\sum_{j=1}^{\left\lfloor\frac{ai+b}{c}\right\rfloor} j-\left\lfloor\frac{ai+b}{c}\right\rfloor\right) \tag 1 \\ & = 2\times \sum_{j=1}^m \sum_{i=0}^n j\times\left[j\le \frac{ai+b}{c}\right]-f(a,b,c,n) \tag 2 \\ & = 2\times \sum_{j=0}^{m-1} (j+1)\left(n-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)-f(a,b,c,n)\tag 3 \\ & = 2\times\sum_{j=0}^{m-1} \left(n(j+1)-j\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor-\left\lfloor\frac{cj+c-b-1}{a}\right\rfloor\right)-f(a,b,c,n) \tag 4 \\ & = nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n)\tag 5 \end{align*} $$

式子 $(1)$:我们按照上述公式将 $\left\lfloor\frac{ai+b}{c}\right\rfloor^2$ 进行变形。

式子 $(2)$:变换求和顺序,将 $j$ 的上界限定在 $\left\lfloor\frac{ai+b}{c}\right\rfloor$ 即 $\frac{ai+b}{c}$ 范围内。上述公式中多出来的一项转化为 $f​$ 函数。

式子 $(3)$:改变 $j$ 的上下界,并按照前文的方法(详见 $f$ 或 $g$ 函数的这部分内容)确定 $i$ 的范围。

式子 $(4)​$:将求和式展开。

式子 $(5)​$:将能够直接计算的提出,其余递归计算

  1. 边界条件 $a=0$:答案为 $(n+1)\left\lfloor\frac{b}{c}\right\rfloor^2$。

伪代码

int h(int a, int b, int c, int n) {
    if (!a) {
        return sqr(b / c) * S0(n);
    }
    if (a >= c || b >= c) {
        return sqr(a / c) * S2(n) + sqr(b / c) * S0(n) + (a / c) * (b / c) * 2 * S1(n)
        + h(a % c, b % c, c, n) + 2 * (a / c) * g(a % c, b % c, c, n) + 2 * (b / c) * f(a % c, b % c, c, n);
    }
    int m = (a * n + b) / c;
    return n * m * S0(m) - f(a, b, c, n) - 2 * g(c, c - b - 1, a, m - 1) - 2 * f(c, c - b - 1, a, m - 1);
}

整体求解

如果我们要对 $f,g,h$ 三个函数同时求解,显然是没法记忆化的。

注意到每个函数中都调用了 $(a\bmod c,b\bmod c,c,n)$ 和 $(c,c-b-1,a,m-1)$,那么我们考虑将 $3$ 个函数在同一个函数中同时求解

时间复杂度:$\mathcal O(\log n)$

代码如下:

#include <cstdio>

typedef long long LL;

const int mod = 998244353;
const LL i2 = 499122177, i6 = 166374059;

struct Ans {
    LL f, g, h;
};

LL sqr(LL x) {
    return x * x % mod;
}
LL S0(LL n) {
    return n + 1;
}
LL S1(LL n) {
    return n * (n + 1) % mod * i2 % mod;
}
LL S2(LL n) {
    return n * (n + 1) % mod * (n + n + 1) % mod * i6 % mod;
}
Ans solve(LL a, LL b, LL c, LL n) {
    Ans now;
    LL A = a / c, B = b / c;
    if (!a) {
        now.f = S0(n) * B % mod;
        now.g = S1(n) * B % mod;
        now.h = S0(n) * sqr(B) % mod;
        return now;
    }
    if (a >= c || b >= c) {
        Ans nxt = solve(a % c, b % c, c, n);
        now.f = S1(n) * A % mod + S0(n) * B % mod + nxt.f;
        now.g = S2(n) * A % mod + S1(n) * B % mod + nxt.g;
        now.h = S2(n) * sqr(A) % mod + S0(n) * sqr(B) % mod + 2 * S1(n) * A % mod * B % mod
        + (B + B) * nxt.f % mod + (A + A) * nxt.g % mod + nxt.h;
        now.f %= mod, now.g %= mod, now.h %= mod;
        return now;
    }
    LL m = (a * n + b) / c;
    Ans nxt = solve(c, c - b - 1, a, m - 1);
    now.f = n * m % mod - nxt.f;
    now.g = (m * n % mod * S0(n) % mod - nxt.f - nxt.h) * i2 % mod;
    now.h = n * m % mod * S0(m) % mod - now.f - 2 * (nxt.f + nxt.g);
    now.f %= mod, now.g %= mod, now.h %= mod;
    return now;
}
int main() {
    int m;
    scanf("%d", &m);
    for (int i = 1; i <= m; i++) {
        LL a, b, c, n;
        scanf("%lld%lld%lld%lld", &n, &a, &b, &c);
        Ans ans = solve(a, b, c, n);
        ans.f = (ans.f % mod + mod) % mod;
        ans.g = (ans.g % mod + mod) % mod;
        ans.h = (ans.h % mod + mod) % mod;
        printf("%lld %lld %lld\n", ans.f, ans.g, ans.h);
    }
    return 0;
}

高次方和求法

首先,我们有如下公式:

$$ \begin{aligned} &\sum_{i=1}^n i^1=\frac{n(n+1)}{2}\\ &\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}{6}\\ &\sum_{i=1}^n i^3=\frac{n^2(n+1)^2}{4}\\ &\sum_{i=1}^n i^4=\cdots \end{aligned} $$

接下来我们就分析一下这些公式是怎么得到的。

根据

$$ \binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k} $$

我们可以得到

$$ \sum_{i=1}^n \binom{i}{k}=\binom{n+1}{k+1} $$

具体证明只需要将右侧式子按照组合数的基本公式 $\binom{n}{k}=\binom{n-1}{k-1}+\binom{n-1}{k}​$ 暴力展开。

那么我们可以根据这个结论得到:

$$ \begin{align*} \sum_{i=1}^n i^1 & = \sum_{i=1}^n i=\sum_{i=1}^n \binom{i}{1}=\binom{n+1}{2} \\ & =\frac{n(n+1)}{2} \end{align*} $$

对于平方求和,证明如下:

$$ \begin{align*} \sum_{i=1}^n i^2 &= \sum_{i=1}^n \left[i(i-1)+i\right]=2\sum_{i=1}^n \binom{i}{2}+\sum_{i=1}^n \binom{i}{1}=2\binom{n+1}{3}+\binom{n+1}{2} \\ &=\frac{n(n+1)(2n+1)}{6} \end{align*} $$

再证明一下立方求和:

$$ \begin{align*} \sum_{i=1}^n i^3&=\sum_{i=1}^n \left[i(i-1)(i-2)+3i(i-1)+i\right]=6\sum_{i=1}^n\binom{i}{3}+6\sum_{i=1}^n\binom{i}{2}+\sum_{i=1}^n\binom{i}{1}\\ &=6\binom{n+1}{4}+6\binom{n+1}{3}+\binom{n+1}{2} \\ &=\frac{n^2(n+1)^2}{4} \end{align*} $$

更高的次数也能按照以上规则写出求和式!

如果用代码实现的话,可以采用拉格朗日插值(详见「算法笔记」拉格朗日插值)求出 $k$ 次方和,时间复杂度最优可以达到 $\mathcal O(\min(k, n, P))$,其中 $n$ 是数字个数,$P$ 是模数。


习题

发表评论