类欧几里得算法

做到了一个$\rm cf$题需要用到这个东西我就去学了下…

以前觉得这玩意好毒瘤，不过其实也不难。

类欧几里得算法一般是用来快速计算以下几个函数的（当然有其他的作用但是一般不会出现，以后看到了在写吧）：

$$\begin{align*} f(a,b,c,n)&=\sum_{i=0}^{n}\lfloor\frac{ai+b}{c}\rfloor\\ g(a,b,c,n)&=\sum_{i=0}^{n}i\lfloor\frac{ai+b}{c}\rfloor\\ h(a,b,c,n)&=\sum_{i=0}^{n}\lfloor \frac{ai+b}{c}\rfloor ^2 \end{align*}$$

其复杂度为$O(\log n)$。

计算F

我们分情况讨论，先考虑$a\geqslant c$或者$b\geqslant c$的情况。

注意到根据整除的性质我们可以把式子写成这样：

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

注意到这里递归了，而且递归后一定是一个$a<c,b<c$的情况，那么此时用上面的展开就没用了。

注意到这个和式实质上是求一条直线下的整点个数，一开始我们是枚举$x$，然后算$y$，我们现在先枚举$y$，也就是横着的一条，可以得到：

$$\begin{align*} m&=\lfloor\frac{an+b}{c}\rfloor\\ f(a,b,c,n)&=\sum_{i=0}^{n}\sum_{j=1}^{m}[\lfloor\frac{ai+b}{c}\rfloor\geqslant j] \end{align*}$$

推一下后面那个不等式，因为是大于等于号所以整除可以去掉：

$$\begin{align*}\lfloor\frac{ai+b}{c}\rfloor&\geqslant j\\ ai+b&\geqslant jc\\ i&\geqslant \frac{jc-b}{a}\\ i&> \lfloor\frac{jc-b-1}{a}\rfloor\\ \end{align*}$$

那么式子就可以化简了：

$$\begin{align*} f(a,b,c,n)&=\sum_{j=0}^{m-1}(n-\lfloor\frac{jc+c-b-1}{a}\rfloor)\\ f(a,b,c,n)&=nm-f(c,c-b-1,a,m-1) \end{align*}$$

观察下递归时$a,c$的变化：$(a,c)\to (a\bmod c,a)\to (a,a\bmod c)$。

这就是欧几里得算法的递归过程，所以最终$a$一定会递归到$0$，此时式子可以直接计算，递归次数为$O(\log n)$。

计算G

和上面思路相同，当$a\geqslant c$或$b\geqslant c$时：

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

否则：

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

其中$f$我们已经会算了，这里还用到了$h$，所以我们在推推$h$。

计算H

当$a\geqslant c$或$b\geqslant c$时，思路一样，但是式子有点复杂注意别搞错：

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

其中函数后面参数都是$(a\bmod c,b\bmod c,c,n)$，我省略没写了。

第二种情况这时候就有点难算了，因为平方无法处理，这里有一个小技巧可以把平方拆开：

$$n^2=-n+2\sum_{i=0}^{n}i$$

那么就很简单了，一个整点$(x,y)$会被统计$y$的贡献，式子写出来就是：

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

代码

如果按照上面写三个函数互相调用显然是不优的，因为一个状态有可能被算很多次，当然可以采取记忆化，但是这样会多一个查找的$\log $。

有一个小技巧就是，我们可以写一个$\rm dfs$，返回三个值，分别是$f,g,h$，注意到这三个函数递归的时候参数都是一样的，所以这样只会递归$O(\log n)$次，并且每次递归都是$O(1)$算。

所以总复杂度$O(\log n)$。

代码长这样：

#include<bits/stdc++.h>
using namespace std;

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define lf double
#define ll long long 

#define pii pair<int,int >
#define vec vector<int >

#define pb push_back
#define mp make_pair
#define fr first
#define sc second

#define FOR(i,l,r) for(int i=l,i##_r=r;i<=i##_r;i++)

const int maxn = 1e6+10;
const int inf = 1e9;
const lf eps = 1e-8;
const int mod = 998244353;
const int inv2 = 499122177;
const int inv6 = 166374059;

struct data {int f,g,h;};

data calc(int a,int b,int c,int n) {
    data ans;
    int s1=1ll*n*(n+1)%mod*inv2%mod,s2=1ll*n*(n+1)%mod*(n*2+1)%mod*inv6%mod;
    if(!a) {
        ans.f=1ll*(n+1)*(b/c)%mod;
        ans.g=1ll*s1*(b/c)%mod;
        ans.h=1ll*(b/c)*(b/c)%mod*(n+1)%mod;
        return ans;
    }
    if(a>=c||b>=c) {
        data res=calc(a%c,b%c,c,n);
        ans.f=(0ll+res.f+1ll*s1*(a/c)%mod+1ll*(n+1)*(b/c)%mod)%mod;
        ans.g=(0ll+res.g+1ll*s2*(a/c)%mod+1ll*s1*(b/c)%mod)%mod;
        ans.h=(0ll+res.h+1ll*(a/c)*res.g*2%mod+1ll*(b/c)*res.f*2%mod)%mod;
        ans.h=(0ll+ans.h+1ll*s2*(a/c)%mod*(a/c)%mod+1ll*n*(n+1)%mod*(b/c)%mod*(a/c)%mod)%mod;
        ans.h=(ans.h+1ll*(n+1)*(b/c)%mod*(b/c)%mod)%mod;
        return ans;
    }int m=(1ll*a*n+b)/c;
    data res=calc(c,c-b-1,a,m-1);
    ans.f=(1ll*m*n%mod-res.f+mod)%mod;
    ans.g=(1ll*m*n%mod*(n+1)%mod-res.h-res.f)%mod;ans.g=(1ll*ans.g*inv2%mod+mod)%mod;
    ans.h=(0ll+1ll*m*n%mod*(m+1)%mod-2ll*res.g-2ll*res.f-ans.f)%mod;ans.h=(ans.h+mod)%mod;
    return ans;
}

int main() {
    int t;read(t);
    while(t--) {
        int a,b,c,n;read(n),read(a),read(b),read(c);
        data ans=calc(a,b,c,n);
        printf("%d %d %d\n",ans.f,ans.h,ans.g);
    }
    return 0;
}

看起来很复杂其实就是对着式子抄了一遍，提交至洛谷模板可$\rm AC$。