类欧几里得学习笔记

作用

实际上是求一些关于整除的和式，比如

[sum_{i=0}^nlfloorfrac{ai+b}{c} floor ， sum_{i=0}^nlfloorfrac{ai+b}{c} floor^2 ， sum_{i=0}^nilfloorfrac{ai+b}{c} floor ]

（这里的三个例子来自luoguP5170）

虽然叫这个名字，但其实跟 gcd 没什么关系，只是复杂度和复杂度分析相同。所以复杂度是 (O(logmax(a,c))) 。

第一个式子实际上是求直线下方在第一象限内的整点数

算法思想

没啥思想，就是推式子，运用一些基本的数学技巧。

推上面那三个式子需要的 trick：

[lfloorfrac{ai+b}{c} floor=lfloorfrac{(amod c)i+(b mod c)}{c} floor + ilfloorfrac{a}{c} floor+lfloorfrac{b}{c} floor \ X^2=2sum_{d=0}^nd-X ]

第一个式子非常好证明，只需要把 a 和 b拆成 kc+r 的形式即可。

第二个式子是通过添加一个求和来消除二次方，同理三次方也可以用类似方法消掉。之后可以通过交换求和顺序来化简。

这种东西分为两类来考虑： a 或 b 大于等于 c ；a 和 b 都小于 c。

(a ge c or b ge c)

运用第一个式子很容易从 (a,b,c,n) 递归到 (a%c,b%c,c,n) 。

(a,b < c)

此时发现 (lfloorfrac{ai+b}{c} floorle n) ，可以用这个缩小范围。令 (m=lfloorfrac{an+b}{c} floor) ：

[f(a,b,c,n) = sum_{i=0}^nlfloorfrac{ai+b}{c} floor \ = sum_{i=0}^nsum_{d=0}^{m-1}[lfloorfrac{ai+b}{c} floor ge d+1] \ = sum_{d=0}^{m-1}sum_{i=0}^n[i>lfloorfrac{cd+c-b-1}{a} floor] \ = sum_{d=0}^{m-1}n - lfloorfrac{cd+c-b-1}{a} floor\ = n*m - f(c,c-b-1,a,m-1) ]

重点其实在第一步，这个谓词很关键（实际上是 (X=sum_{d=1}^X1)）。后面的就很显然了。再推一个。

[g(a,b,c,n) = sum_{i=0}^nlfloorfrac{ai+b}{c} floor^2 \ = sum_{i=0}^n2sum_{d=1}^{lfloorfrac{ai+b}{c} floor}d - lfloorfrac{ai+b}{c} floor\ = 2sum_{d=0}^{m-1}d+1sum_{i=0}^n[lfloorfrac{ai+b}{c} floor ge d+1] - f(a,b,c,d) ]

这里用到了上述的第二个式子。其他的都类似，不再赘述。

应用

luoguP5170

其中

[f=sum_{i=0}^nlfloorfrac{ai+b}{c} floor ， g=sum_{i=0}^nlfloorfrac{ai+b}{c} floor^2 ，h= sum_{i=0}^nilfloorfrac{ai+b}{c} floor ]

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

typedef long long LL;
typedef long double LD;
typedef pair<int,int> pii;
typedef pair<LL,int> pli;
const int SZ = 1e6 + 10;
const LL INF = 1e15 + 10;
const int mod = 998244353;
const LD eps = 1e-8;

LL read() {
    LL n = 0;
    char a = getchar();
    bool flag = 0;
    while(a > '9' || a < '0') { if(a == '-') flag = 1; a = getchar(); }
    while(a <= '9' && a >= '0') { n = n * 10 + a - '0',a = getchar(); }
	if(flag) n = -n;
	return n;
}


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

struct LO {

    const LL inv2 = 499122177;
    const LL inv6 = 166374059;

    LL F1(LL n) {
        n %= mod;
        return n * (n+1) % mod * inv2 % mod;
    }

    LL F2(LL n) {
        n %= mod;
        return n * (n+1) % mod * (2*n+1) % mod * inv6 % mod;
    }

    node solve(LL a,LL b,LL c,LL n) {
        LL f,g,h;
        if(a == 0) {
            n %= mod;
            f = (b/c%mod) * (n+1) % mod;
            g = (b/c%mod) * (b/c%mod)%mod * (n+1) % mod;
            h = (b/c%mod) * F1(n) % mod;
        }
        else if(a>=c || b>=c) {
            node t = solve(a%c,b%c,c,n);
            f = (t.f + F1(n)*(a/c%mod)%mod + ((n+1)%mod)*(b/c%mod)%mod) % mod;
            g = (F2(n)*(a/c%mod)%mod*(a/c%mod)%mod + ((n+1)%mod)*(b/c%mod)%mod*(b/c%mod)%mod + t.g
                  + 2*F1(n)*(a/c%mod)%mod*(b/c%mod)%mod + 2*(a/c%mod)*t.h%mod + 2*(b/c%mod)*t.f%mod) % mod;
            h = (t.h + F2(n)*(a/c%mod)%mod + F1(n)*(b/c%mod)%mod) % mod;
        }
        else {
            LL m = (a*n+b)/c;
            node t = solve(c,c-b-1,a,m-1);
            m %= mod;
            f = (n*m%mod - t.f) % mod;
            g = (2*n*F1(m)%mod - 2*(t.h+t.f) - f) % mod;
            h = (F1(n)*m%mod - (t.g+t.f)*inv2%mod) % mod;
        }
        f += mod; f %= mod;
        g += mod; g %= mod;
        h += mod; h %= mod;
        return (node){f,g,h};
    }
}lo;

int main() {
    int T = read();
    while(T --) {
        LL n = read(),a = read(),b = read(),c = read();
        node ans = lo.solve(a,b,c,n);
        printf("%lld %lld %lld
",ans.f,ans.g,ans.h);
    }
}