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

发表于 | | 阅读次数:
字数统计: 3.1k 字
类欧几里德算法基于欧几里德算法，本质是通过取模缩小问题规模，递归求解。

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

定义

推导

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

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

    式子 $(2)$：第一项直接 $O(1)$ 计算，第二项递归计算。

  2. 当 $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$，枚举所有点并判断是否在线段下方。

    式子 $(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$ 函数，递归计算。

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

伪代码

1
2
3
4
5
6
 
int f(int a,int b,int c,int n) {
    if(!a) return b/c*(n+1);
    if(a>=c||b>=c) return a/c*n*(n+1)/2+b/c*(n+1)+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)$

定义

推导

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

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

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

  2. 当 $a<c​$ 且 $b<c​$ 时：

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

    式子 $(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)​$：第一项可以 $O(1)​$ 计算，第二项和第三项发现是 $f​$ 和 $h​$ 函数，可以递归计算。

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

伪代码

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

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

定义

推导

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

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

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

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

  2. 当 $a<c$ 且 $b<c$ 时：

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

    式子 $(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)​$：将能够直接计算的提出，其余递归计算。

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

伪代码

1
2
3
4
5
6
 
int h(int a,int b,int c,int n) {
    if(!a) return sqr(b/c)*(n+1);
    if(a>=c||b>=c) return sqr(a/c)*n*(n+1)*(2*n+1)/6+sqr(b/c)*(n+1)+(a/c)*(b/c)*n*(n+1)+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*(m+1)-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$ 个函数在同一个函数中同时求解。

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

代码如下：

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
 
#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 S2(LL n) {
    return n*(n+1)%mod*i2%mod;
}
LL S3(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=(n+1)*B%mod;
        now.g=S2(n)*B%mod;
        now.h=(n+1)*sqr(B)%mod;
        return now;
    }
    if(a>=c||b>=c) {
        Ans nxt=solve(a%c,b%c,c,n);
        now.f=S2(n)*A%mod+(n+1)*B%mod+nxt.f;
        now.g=S3(n)*A%mod+S2(n)*B%mod+nxt.g;
        now.h=S3(n)*sqr(A)%mod+(n+1)*sqr(B)%mod+2*S2(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*(n+1)%mod-nxt.f-nxt.h)*i2%mod;
    now.h=n*m%mod*(m+1)%mod-now.f-2*(nxt.f+nxt.g);
    now.f%=mod,now.g%=mod,now.h%=mod;
    return now;
}
int main() {
    int T;
    for(scanf("%d",&T);T--;) {
        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;
}

习题

附：高次方和求法

首先，我们有如下公式：

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

根据

我们可以得到

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

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

对于平方求和，证明如下：

再证明一下立方求和：

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

如果用代码实现的话，可以采用拉格朗日插值（详见「算法笔记」拉格朗日插值）求出 $k$ 次方和，时间复杂度最优可以达到 $O(k)$。

0%