题解 P5170 【【模板】类欧几里得算法】

· · 题解

emmmm……虽说这家伙叫“类欧几里得”，但实际上它和欧几里得算法的关系真的不是很大……

欧几里得算法，其实就是辗转相除法求 \gcd

辗转相除法怎么写的？是用递归计算的。

于是类欧几里得也是用递归去求，除此之外基本没有关系了……

类欧几里得算法怎么用？我们来看看这个柿子：

\sum_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor

这个柿子带有下取整符号，直接算不大好算，怎么办呢？我们分类讨论

首先先设这一坨东西为 f(a,b,c,n) ，毕竟我们待会要递归计算嘛

  1. 当 a=0 时，就是 (n+1)\lfloor\dfrac{b}{c}\rfloor

  2. 当 a\ge c 或 b \ge c 时，我们就可以把整除的部分单独分离出来

\begin{aligned}&\quad f(a,b,c,n)\\&=\sum_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor\\&=\sum_{i=0}^n(\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor+i\lfloor\dfrac{a}{c}\rfloor+\lfloor\dfrac{b}{c}\rfloor)\\&=\sum_{i=0}^n\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor+\dfrac{n(n+1)}{2}\lfloor\dfrac{a}{c}\rfloor+(n+1)\lfloor\dfrac{b}{c}\rfloor\\&=f(a\bmod c,b\bmod c,c,n)+\dfrac{n(n+1)}{2}\lfloor\dfrac{a}{c}\rfloor+(n+1)\lfloor\dfrac{b}{c}\rfloor\end{aligned}
  1. 当 a<c 且 b<c 时，我们设 m=\lfloor\dfrac{an+b}{c}\rfloor，并且做几（十）步玄学的转化：
\begin{aligned}&\quad f(a,b,c,n)\\&=\sum_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor\\&=\sum_{i=0}^n\sum_{j=1}^m[j\le\lfloor\dfrac{ai+b}{c}\rfloor]\\&=\sum_{i=0}^n\sum_{j=0}^{m-1}[j+1\le\lfloor\dfrac{ai+b}{c}\rfloor]\\&=\sum_{i=0}^n\sum_{j=0}^{m-1}[cj+c\le ai+b]\\&=\sum_{i=0}^n\sum_{j=0}^{m-1}[cj+c< ai+b+1]\\&=\sum_{j=0}^{m-1}\sum_{i=0}^n[cj+c< ai+b+1]\\&=\sum_{j=0}^{m-1}\sum_{i=0}^n[ai>cj+c-b-1]\\&=\sum_{j=0}^{m-1}\sum_{i=0}^n[i>\lfloor\dfrac{cj+c-b-1}{a}\rfloor]\\&=\sum_{j=0}^{m-1}(n-\lfloor\dfrac{cj+c-b-1}{a}\rfloor)\\&=mn-\sum_{j=0}^{m-1}\lfloor\dfrac{cj+(c-b-1)}{a}\rfloor\\&=mn-\; f(c,c-b-1,a,m-1)\end{aligned}

于是乎，我们就把这玩意整成了一个递归的柿子，可以在 O(\log n) 的时间内求出来啦！

f(a,b,c,n)=\begin{cases}(n+1)\lfloor\dfrac{b}{c}\rfloor,a=0\\f(a\bmod c,b\bmod c,c,n)+\dfrac{n(n+1)}{2}\lfloor\dfrac{a}{c}\rfloor+(n+1)\lfloor\dfrac{b}{c}\rfloor,a\ge c \text{或} b\ge c\\mn-\; f(c,c-b-1,a,m-1),a<c \text{且} b<c\end{cases}

经过我们不懈的努力，我们终于求出了 f 函数！

那么……这个呢？？？

g(a,b,c,n)=\sum_{i=0}^ni\lfloor\dfrac{ai+b}{c}\rfloor h(a,b,c,n)=\sum_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor^2

推完了 f 函数，相信各位已经有经验了：首先我们还是要分类讨论+取模

g(a,b,c,n)=g(a \bmod c,b\bmod c,c,n)+\dfrac{n(n+1)(2n+1)}{6}\lfloor\dfrac{a}{c}\rfloor+\dfrac{n(n+1)}{2}\lfloor\dfrac{b}{c}\rfloor

由于多乘了一个 i，原先的 n+1 变成了等差数列和，原来的等差数列和变成了平方和……

接下来就是单独研究 a<c 且 b<c 的情况了

还是原来的配方，还是熟悉的味道：把柿子转化成不等关系再转回来！

为了满足各位懒癌患者的需求，我们把上次最后算出来的那个 \lfloor\dfrac{cj+(c-b-1)}{a}\rfloor 设为 t

于是我们可以得到：

\begin{aligned}&\quad g(a,b,c,n)\\&=\sum_{i=0}^n\sum_{j=0}^{m-1}(i\times[j+1\le\lfloor\dfrac{ai+b}{c}\rfloor])\\&=\sum_{j=0}^{m-1}\sum_{i=0}^n(i\times[i>t])\\&=\sum_{j=0}^{m-1}\sum_{i=t+1}^ni\\&=\sum_{j=0}^{m-1}\dfrac{(t+1+n)(n-t)}{2}\\&=\dfrac{1}{2}\sum_{j=0}^{m-1}(n^2+n-t^2-t)\\&=\dfrac{1}{2}[mn(n+1)-\sum_{j=0}^{m-1}t^2-\sum_{j=0}^{m-1}t]\\&=\dfrac{1}{2}[mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)]\end{aligned}

这样我们就把 g 函数也推导出来啦

话说我们在推 g 函数的时候，最后一步竟然出现了 h 函数，真是奇妙啊QAQ

最后我们要推导 h 函数的递归形式咯！

车到山前必有路，推就完了！奥利给！

按照国际惯例，依然是先取模：

\begin{aligned}&\quad h(a,b,c,n)\\&=\sum_{i=0}^n\lfloor\dfrac{ai+b}{c}\rfloor^2\\&=\sum_{i=0}^n(\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor+i\lfloor\dfrac{a}{c}\rfloor+\lfloor\dfrac{b}{c}\rfloor)^2\\&=\sum_{i=0}^n[\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor^2+i^2\lfloor\dfrac{a}{c}\rfloor^2+\lfloor\dfrac{b}{c}\rfloor^2]+2\sum_{i=0}^n[\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor i\lfloor\dfrac{a}{c}\rfloor+\lfloor\dfrac{(a\bmod c)i+(b\bmod c)}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor+i\lfloor\dfrac{a}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor]\\&=h(a\bmod c,b\bmod c,c,n)+\dfrac{n(n+1)(2n+1)}{6}\lfloor\dfrac{a}{c}\rfloor^2+(n+1)\lfloor\dfrac{b}{c}\rfloor^2+2[\lfloor\dfrac{a}{c}\rfloor \;g(a\bmod c,b\bmod c,c,n)+\lfloor\dfrac{b}{c}\rfloor\;f(a\bmod c,b\bmod c,c,n)+\dfrac{n(n+1)}{2}\lfloor\dfrac{a}{c}\rfloor\lfloor\dfrac{b}{c}\rfloor]\end{aligned}

总之，就是将平方式利用三元的完全平方公式展开，再把它分成“平方”和“乘积的2倍”的形式，分别转化成递归形式。

在考虑 a<c 且 b<c 的情况之前，让我们来看这个柿子，非常好理解：

n^2=2\dfrac{n(n+1)}{2}-n=2\sum_{i=1}^n i-n

就这样，我们把平方中“积”的形式成功转化成了“和”的形式，就不怕出现 \sum\times\sum 的麻烦啦！

h(a,b,c,n)=\sum_{i=0}^n[2\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor} j-\lfloor\dfrac{ai+b}{c}\rfloor]=2\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor}j-f(a,b,c,n)

接下来我们只需要着重化简前边那一块就行了，以下的 m 还是那个 m，t 还是那个 t

还是通过不等关系的限制来转化，但是因为这次 j 的范围与 i 有关，转化会稍微麻烦些：

\begin{aligned}&\quad\sum_{i=0}^n\sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor}j\\&=\sum_{i=0}^n\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}(j+1)\\&=\sum_{i=0}^n\sum_{j=0}^{m-1}(j+1)[j< \lfloor\frac{ai+b}{c}\rfloor]\\&=\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[j< \lfloor\frac{ai+b}{c}\rfloor]\\&=\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n[i>t]\\&=\sum_{j=0}^{m-1}(j+1)(n-t)\\&=\sum_{j=0}^{m-1}(j+1)n-\sum_{j=0}^{m-1}(j+1)t\\&=\dfrac{nm(m+1)}{2}-\sum_{j=0}^{m-1}jt-\sum_{j=0}^{m-1}t\\&=\dfrac{nm(m+1)}{2}-g(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)\end{aligned}

于是，我们可以得到 h 函数的递归形式：

h(a,b,c,n)=nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n) \huge\color{orangered}{\text{终\quad于\quad推\quad完\quad了\quad！！！！！}}

最后就是代码实现的方面，

  1. 我们注意到 g 函数和 h 函数它们俩是互相递归计算的，所以我们需要把它们同时计算出来！这里我用了一个结构体来存储计算结果，大家可以来参考一下。

  2. 参照欧几里得算法（辗转相除法求 \gcd）可以得到这个算法的时间复杂度也是 O(\log n) 。

  3. 还有还有，细节决定成败。

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
struct ccf{
    ll f,g,h;
};
const ll mod=998244353ll,inv2=499122177ll,inv6=166374059ll;
ccf wyxakioi(ll a,ll b,ll c,ll n){
    ll s1=(n%mod*(n+1)%mod)*inv2%mod,s2=((n%mod*(n+1)%mod)*(2*n+1)%mod)*inv6%mod;
    if(a==0){
        ll F=(n+1)*(b/c)%mod;
        ll G=s1*(b/c)%mod;
        ll H=((n+1)*(b/c)%mod)*(b/c)%mod;
        return (ccf){F%mod,G%mod,H%mod};
    }
    else if(a>=c||b>=c){
        ccf qwq=wyxakioi(a%c,b%c,c,n);
        ll F=qwq.f+s1*(a/c)%mod+(n+1)*(b/c)%mod;
        ll G=qwq.g+s2*(a/c)%mod+s1*(b/c)%mod;
        ll H=qwq.h+(s2*(a/c)%mod)*(a/c)%mod+(n+1)*(b/c)%mod*(b/c)%mod+2ll*((a/c)*qwq.g%mod+(b/c)*qwq.f%mod+(s1*(a/c)%mod)*(b/c)%mod)%mod;
        return (ccf){F%mod,G%mod,H%mod};
    }
    else{
        ll m=(a*n+b)/c;
        ccf qwq=wyxakioi(c,c-b-1,a,m-1);
        ll F=(n*m)%mod-qwq.f;
        ll G=(((n*m)%mod)*(n+1)%mod-qwq.h+mod-qwq.f)%mod*inv2%mod;
        ll H=((n*m)%mod)*(m+1)%mod-2ll*qwq.g%mod+mod-2ll*qwq.f%mod+mod-F%mod;
        return (ccf){F%mod,G%mod,H%mod};
    }
}
int main(){
    int t;
    scanf("%d",&t);
    while(t--){
        ll a,b,c,n;
        scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
        ccf ans=wyxakioi(a,b,c,n);
        printf("%lld %lld %lld\n",(ans.f+mod)%mod,(ans.h+mod)%mod,(ans.g+mod)%mod);
    }
}

嘀！又有一道黑题进账啦！