类欧几里得算法小结

基本定义

$f(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor$
$g(a,b,c,n)=\sum_{i=0}^ni\lfloor\frac{ai+b}{c}\rfloor$
$h(a,b,c,n)=\sum_{i=0}^n\lfloor\frac{ai+b}{c}\rfloor^2$
$m=\lfloor\frac{an+b}{c}\rfloor$
为了方便，接下来下取整都不写下取整，出现除法就是整除

推f

当a>=c或b>=c的时候， 容易得到
$f(a,b,c,n)=(a/c)*n*(n+1)/2+(b/c)*(n+1)+f(a\%c,b\%c,c,n)$
如果a和b均小于c呢？
$f(a,b,c,n)=\sum_{i=0}^n\sum_{j=1}^m[(ai+b)/c>=j]$
这条式子可以这样理解，通过几何意义。
f的式子的意义其实是一条倾斜直线与x=0和x=n两条铅垂线组成的直角梯形内整点数量（不包括纵坐标为0的点，包括在线上的点）。
接下来我们这条式子就是在枚举一条纵线，统计上面整点的个数，我们现在改为枚举一条横线，统计横线上整点的个数，就是这条式子了。
$f(a,b,c,n)=\sum_{i=0}^n\sum_{j=0}^{m-1}[(ai+b)/c>=j+1]$
$f(a,b,c,n)=\sum_{i=0}^n\sum_{j=0}^{m-1}[ai>=cj+c-b]$
$f(a,b,c,n)=\sum_{i=0}^n\sum_{j=0}^{m-1}[ai>cj+c-b-1]$
$f(a,b,c,n)=\sum_{i=0}^n\sum_{j=0}^{m-1}[i>(cj+c-b-1)/a]$
$f(a,b,c,n)=\sum_{j=0}^m(n-(cj+c-b-1)/a)$
$f(a,b,c,n)=nm-f(c,c-b-1,a,m-1)$
继续递归处理。
可以发现只看第一和第三维表示状态的话，每次由(a,c)变成(c,a%c)，因此复杂度与欧几里得算法一致。

推g

a>=c或b>=c的时候，容易得到
$g(a,b,c,n)=(a/c)*n*(n+1)*(2n+1)/6+(b/c)*n*(n+1)/2+g(a\%c,b\%c,c,n)$
注意这里0到n的二次方和公式就是$n*(n+1)*(2n+1)/6$
接下来我们来推a和b均小于c的情况。
首先先和f一样把下取整那玩意变换一下
$g(a,b,c,n)=\sum_{i=0}^ni\sum_{j=1}^m[(ai+b)/c>=j]$
$g(a,b,c,n)=\sum_{i=0}^ni\sum_{j=0}^{m-1}[i>(cj+c-b-1)/a]$
$g(a,b,c,n)=1/2*\sum_{j=0}^{m-1}(n+1+(cj+c-b-1)/a)*(n-(cj+c-b-1)/a)$
这里看的有点难受，f的时候转化的其实是0次方和，而g是1次方和，所以这里用了个求和公式。
$g(a,b,c,n)=1/2*\sum_{j=0}^{m-1}n(n+1)-(cj+c-b-1)/a-[(cj+c-b-1)/a]^2$
拆出来了
$g(a,b,c,n)=1/2*[n(n+1)m-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1)]$
看来g比较麻烦，我们还要会推h。

推h

a>=c或b>=c的时候，h也是很麻烦的……
三项式的平方，拆出来大概是这样
$h(a,b,c,n)=(a/c)^2*n(n+1)(2n+1)/6+(b/c)^2*(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)$
接下来a和b均小于c的情况，我们要转化一下思路了。
如何获得一个n^2？
$n^2=2*\frac{n(n+1)}{2}-n=2\sum_{i=0}^ni-n$
有了思路我们来推h
$h(a,b,c,n)=\sum_{i=0}^n(2*\sum_{j=1}^{(ai+b)/c}j-(ai+b)/c)$
可以想到交换主体。
$h(a,b,c,n)=2*\sum_{j=0}^{m-1}(j+1)*\sum_{i=0}^n[(ai+b)/c>=j+1]-f(a,b,c,n)$
$h(a,b,c,n)=2*\sum_{j=0}^{m-1}(j+1)*\sum_{i=0}^n[i>(cj+c-b-1)/a]-f(a,b,c,n)$
$h(a,b,c,n)=2*\sum_{j=0}^{m-1}(j+1)*(n-(cj+c-b-1)/a)-f(a,b,c,n)$
$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)$

实现

我的实现是直接返回f、g和h分别算出来的值，这样是比较简单的。
边界是a为0，返回三个0。
upd:
好像边界我写错了，大家随手推推吧。
写得挺丑，加了的强制转换其实没用……

dong likegcd(ll a,ll b,ll c,ll n){
    if (!a){
        gjx.f=gjx.g=gjx.h=0;
        return gjx;
    }
    if (a>=c||b>=c){
        zlt=likegcd(a%c,b%c,c,n);
        gjx.f=((ll)(a/c)*n%mo*(n+1)%mo*ni2%mo+(ll)(b/c)*(n+1)%mo)%mo;
        (gjx.f+=zlt.f)%=mo;
        gjx.g=((ll)(a/c)*n%mo*(n+1)%mo*(2*n+1)%mo*ni6%mo+(ll)(b/c)*(n+1)%mo*n%mo*ni2%mo)%mo;
        (gjx.g+=zlt.g)%=mo;
        gjx.h=(ll)(a/c)*(a/c)%mo*n%mo*(n+1)%mo*(2*n+1)%mo*ni6%mo;
        (gjx.h+=(ll)(b/c)*(b/c)%mo*(n+1)%mo)%=mo;
        (gjx.h+=(ll)(a/c)*(b/c)%mo*n%mo*(n+1)%mo)%=mo;
        (gjx.h+=(ll)2*(a/c)%mo*zlt.g%mo)%=mo;
        (gjx.h+=(ll)2*(b/c)%mo*zlt.f%mo)%=mo;
        (gjx.h+=zlt.h)%=mo;
        return gjx;
    }
    ll m=((ll)a*n+(ll)b)/(ll)c;
    zlt=likegcd(c,c-b-1,a,m-1);
    gjx.f=((ll)n*m%mo-(ll)zlt.f)%mo;
    gjx.g=((ll)(n+1)*n%mo*m%mo-(ll)zlt.f-(ll)zlt.h)%mo;
    gjx.g=(ll)gjx.g*ni2%mo;
    gjx.h=((ll)n*m%mo*(m+1)%mo-(ll)2*zlt.g%mo-(ll)2*zlt.f%mo-(ll)gjx.f)%mo;
    return gjx;
}