【算法】类欧几里得算法 ——quhengyi11

类欧几里得算法能够用 $logn​$的时间进行如下求和
$$\sum_{i = 0}^n \left \lfloor \frac{ai+b}c \right \rfloor$$

我们首先记

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

$S_2(n)=\sum_{i=1}^n i=\frac{n(n+1)}2$

$S_3(n)=\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}6$

当 $a\geq c \;||\;b \geq 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 \lfloor \frac{\left (\lfloor \frac{a}{c} \right \rfloor \times c+a \% c)i+(\left\lfloor\frac{b}{c}\right\rfloor\times c+b\%c)}{c} \right \rfloor\\ =&S_2(n)\times\left\lfloor\frac{a}{c}\right\rfloor+(n+1)\times\left\lfloor\frac{b}{c}\right\rfloor+f(a\%c,b\%c,c,n) \end{align}$$

当 $a<c \;\&\&\;b < c$，记 $m=\frac{an+b}c$

$$\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\lfloor\frac{ai+b}{c}\right\rfloor\geq j]\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[\left\lfloor\frac{ai+b}{c}\right\rfloor\geq j+1]\quad&(1)\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[{ai}\geq jc+c-b]\quad &(2)\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[ai> jc+c-b-1]\quad&(3)\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[i> \frac{jc+c-b-1}a]\quad&(4)\\ =&\sum_{j=0}^{m-1}n-\left\lfloor\frac{jc+c-b-1}a\right\rfloor\\ =&nm-\sum_{j=0}^{m-1}\left\lfloor\frac{jc+c-b-1}a\right\rfloor\\ =&nm-f(c,c-b-1,a,m – 1) \end{align}$$

这里我们解释一些细节

首先由式子 $(1)->(2)$为啥不能这样写？

$$\begin{align} =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[\left\lfloor\frac{ai+b}{c}\right\rfloor> j]\quad&(1′)\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[{ai}> jc-b]\quad &(2’)\\ \end{align}$$
我们可以假象一条数轴，整除 $d$这个操作，相当于将这条数轴切成一段一段的，每一段的优先级是一样的，并且每一段长度都为 $d$，而且每个数 $n$的优先级就是 $\left\lfloor\frac{n}{d}\right\rfloor$，显然对于一个 $n$，它所在段的最小的数就是 $\left\lfloor\frac{n}{d}\right\rfloor \times d$，而 $(1)$式到 $(2)$式代表的意义就是将 $j$映射到它在整除 $d$意义下，第 $j$段的最小值，显然只要 $\left\lfloor\frac{ai+b}{c}\right\rfloor$所在段 $\geq j$即对答案有贡献，是与原式子等价的

而 $(1′)->(2′)$就咕咕咕了，当 $j=\left\lfloor\frac{ai+b}{c}\right\rfloor$的时候，它显然对 $(1′)$式是没有贡献的，而如果 $c$不整除 $ai+b$的时候，$ai+b>jc$是显然成立的，于是它又对 $(2′)$式有贡献。

相似的方法，我们可以证明下面这个式子是错的

$$\begin{align} =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[ai\geq jc+c-b]\quad&(3′)\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}[i\geq \frac{jc+c-b}a]\quad&(4′)\\ \end{align}$$

读者可以自己证明

（tip: 若 $ai$和 $jc+c-b$在一个块里并且 $jc+c-b\geq$ai，则 $(3′)$无贡献，$(4′)$有贡献）

（感觉这里自己理解复杂了啊？求更简单的理解方法 qwq）

然后你跑一遍递归就能求出来原式子了。

于是基础的类欧几里得算法就讲完了？

但是我们会止步于此吗？

(大雾)

咕咕咕，所以我们的模板还是没讲完 qwq

来我们把剩下两条式子也推一推

记

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

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

其实手法是相似的呢，我们先求一遍 $g(a,b,c,n)$

当 $a\geq c \;||\;b \geq c$
$$\begin{align} g(a,b,c,n)=&\sum_{i=0}^ni\left \lfloor \frac{ai+b}c \right \rfloor\\ =&S_3(n)\times\left\lfloor\frac{a}{c}\right\rfloor+S_2(n)\times\left\lfloor\frac{b}{c}\right\rfloor+g(a\%c,b\%c,c,n) \end{align}$$
当 $a<c \;\&\&\;b < c$，记 $m=\frac{an+b}c$

$$\begin{align} g(a,b,c,n)=&\sum_{i=0}^n \left\lfloor \frac{ai+b}c \right\rfloor\\ =&\sum_{i = 0}^n \sum_{j=0}^{m-1}i[i> \left\lfloor\frac{jc+c-b-1}a\right\rfloor]\\ =&\sum_{j=0}^{m-1} \frac 1 2(n+\left\lfloor\frac{jc+c-b-1}a\right\rfloor+1)(n – \left\lfloor\frac{jc+c-b-1}a\right\rfloor])\\ =&\frac 1 2\sum_{j=0}^{m-1}n^2-\left\lfloor\frac{jc+c-b-1}a\right\rfloor^2 + n – \left\lfloor\frac{jc+c-b-1}a\right\rfloor\\ =&\frac 1 2(n^2m-h(c,c-b-1,a,m-1)+nm-f(c,c-b-1,a,m-1))\\ =&\frac 1 2 [nm(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)] \end{align}$$

你会震惊地发现，这个 $g(a,b,c,n)$竟然还和 $h$有关系，不过不用惊讶，我们等等推完 $h$你就知道怎么处理了

现在我们来求 $h(a,b,c,n)$

当 $a\geq c \;||\;b \geq 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 \lfloor \frac{\left (\lfloor \frac{a}{c} \right \rfloor \times c+a \% c)i+(\left\lfloor\frac{b}{c}\right\rfloor\times c+b\%c)}{c} \right \rfloor^2\\ =&\sum_{i=0}^n(\left\lfloor\frac{a}{c}\right\rfloor i+\left\lfloor\frac{b}{c}\right\rfloor+\left\lfloor\frac{(a\%c)i+(b\%c)}{c}\right\rfloor)^2\\ =&\sum_{i=0}^n\left\lfloor\frac{a}{c}\right\rfloor^2i^2+\left\lfloor\frac{b}{c}\right\rfloor^2+\left\lfloor\frac{(a\%c)i+(b\%c)}{c}\right\rfloor^2+\\&2\left(\left\lfloor\frac{a}{c}\right\rfloor\left\lfloor\frac{b}{c}\right\rfloor i+\left\lfloor\frac{a}{c}\right\rfloor \left\lfloor\frac{(a\%c)i+(b\%c)}{c}\right\rfloor i + \left\lfloor\frac{b}{c}\right\rfloor\left\lfloor\frac{(a\%c)i+(b\%c)}{c}\right\rfloor\right)\\ =&S_3(n)\left\lfloor\frac{a}{c}\right\rfloor^2+(n+1)\left\lfloor\frac{b}{c}\right\rfloor^2+h(a\%c,b\%c,c,n)+2S_2(n)\left\lfloor\frac{a}{c}\right\rfloor \left\lfloor\frac{b}{c}\right\rfloor\\&+2g(a\%c,b\%c,c,n)\left\lfloor\frac{a}{c}\right\rfloor+2f(a\%c,b\%c,c,n)\left\lfloor\frac{b}{c}\right\rfloor \end{align}$$

当 $a<c \;\&\&\;b < c$，记 $m=\frac{an+b}c$
$$\begin{align} h(a,b,c,n)=&\sum_{i=0}^n\left \lfloor \frac{ai+b}c \right \rfloor^2\ \end{align}$$

等等，如果用相同的手法的话，我们不应该令 $m=(\frac{an+b}c)$喵？但是酱紫做的话你会发现递归就不能保证 $m\leq n$了，然后你就咕咕咕了。

我们可以采用求和降次的手法 (就是 $n$次多项式的 $n$项求和为 $n+1$次，我们只需要多一次枚举就能降次)

在这题我们用的是

$$\left(2\sum_{i=1}^ni\right) – n = n^2$$

代入有

$$\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\sum_{j=1}^{m}j\right)-{\left \lfloor \frac{ai+b}c \right \rfloor}\\ =&\sum_{i=0}^n \left(2\sum_{j=0}^{m-1}j+1\right)-{\left \lfloor \frac{ai+b}c \right \rfloor}\\ =&2\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[\left \lfloor \frac{ai+b}c \right \rfloor \geq j+1]-\sum_{i=0}^n{\left \lfloor \frac{ai+b}c \right \rfloor}\\ =&2\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[ai \geq jc+c-b]-f(a,b,c,n)\\ =&2\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[i > \ \left\lfloor \frac{jc+c-b-1} a\right\rfloor]-f(a,b,c,n)\\ =&2\sum_{j=0}^{m-1}(j+1)(n-\left\lfloor \frac{jc+c-b-1} a\right\rfloor)-f(a,b,c,n)\\ =&2\sum_{j=0}^{m-1}nj+n-j\left\lfloor \frac{jc+c-b-1} a\right\rfloor-\left\lfloor \frac{jc+c-b-1} a\right\rfloor-f(a,b,c,n)\\ =&2nS_2(m-1)+2nm-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n)\\ =&2nS_2(m)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n) \end{align}$$

然后我们的式子就推完了，在实现方面，因为 $g$和 $h$都有互相推导的关系，所以我们可以存一个结构体，在结构体里把三个函数算出来然后递归做

具体实现看代码吧 qwq

#include<bits/stdc++.h>
#define fo(i, a, b) for (Re int i = (a); i <= (b); ++i)
#define fd(i, a, b) for (Re int i = (a); i >= (b); --i)
#define edge(i, u) for (Re int i = head[u], v = e[i].v; i; i = e[i].nxt, v = e[i].v)
#define Re register
#define pb push_back
#define F first
#define S second
#define ll long long
#define lowbit(x) (x & -x)
#define mod 998244353
#define N 15
#define inv2 499122177
#define inv6 166374059
ll n, a, b, c, T;
struct node{
    ll f, g, h;
};
inline ll S2 (ll x) {return (x * (x + 1) >> 1) % mod;}
inline ll S3 (ll x) {return x * (x + 1) % mod * (x << 1 | 1) % mod * inv6 % mod;}
inline ll add (ll x, ll ad) {return (x + ad) % mod;}
inline ll sqr (ll x) {return x * x % mod;}
inline node calc (ll a, ll b, ll c, ll n)
{
    node ret;
    if (a == 0)
    {
        ret.f = (b / c) * (n + 1) % mod;
        ret.g = S2(n) * (b / c) % mod;
        ret.h = (b / c) * ret.f % mod;
//        printf("1:%d %d %d %d    %d %d %d\n", a, b, c, n, ret.f, ret.g, ret.h);
        return ret;
    }
    if (a >= c || b >= c)
    {
        node tmp = calc(a % c, b % c, c, n);
        ll ac = a / c, bc = b / c;
        ret.f = add(add(tmp.f, ac * S2(n) % mod), bc * (n + 1) % mod);
        ret.g = add(add(tmp.g, ac * S3(n) % mod), bc * S2(n) % mod);
        ret.h = add(add(add(add(add(tmp.h, sqr(ac) * S3(n) % mod), sqr(bc) * (n + 1) % mod), 2 * ac * bc % mod * S2(n) % mod), 2 * ac * tmp.g % mod), 2 * bc * tmp.f % mod);
 //       printf("2:%d %d %d %d    %d %d %d\n", a, b, c, n, ret.f, ret.g, ret.h);
        return ret;
    }
    ll m = (a * n + b) / c;
    node tmp = calc(c, c - b - 1, a, m - 1);
    ret.f = add(n * m % mod, -tmp.f);
    ret.g = add(n * m % mod * (n + 1) % mod, -tmp.h - tmp.f) * inv2 % mod;
    ret.h = add(n * m % mod * (m + 1) % mod, -2 * (tmp.f + tmp.g) - ret.f);
 //       printf("3:%d %d %d %d    %d %d %d\n", a, b, c, n, ret.f, ret.g, ret.h);
    return ret;
}
int main ()
{
    scanf("%lld", &T);
    while (T--)
    {
        scanf("%lld %lld %lld %lld", &n, &a, &b, &c);
        node ans = calc(a, b, c, n);
        printf("%lld %lld %lld\n", (ans.f + mod) % mod, (ans.h + mod) % mod, (ans.g + mod) % mod);
    }
    return 0;
}

由于这篇 blog 的公式太多了，如果有错的话请在评论处指出谢谢 qwq