类欧几里得算法小结

基础版

类欧几里得算法被用来处理下面这种式子的值

[f(a,b,c,n)=sum_{i=0}^n lfloor frac{ai+b}{c} floor ]

求值的话我们考虑分类讨论

（1） 若(ageq c)或(b geq c)
注意到(lfloor frac{ac}{b} floor=lfloor frac{(a mod b)*c}{b} floor+clfloor frac{a}{b} floor)，在(ageq c)或(bgeq c)时，由于(lfloorfrac{a}{c} floor)或(lfloorfrac{b}{c} floorgeq 1)，所以利用上式化简可以得到

[egin{aligned} f(a,b,c,n)=&sum_{i=0}^n(frac{(a mod c)*i+(b mod c)}{c}+i lfloor frac{a}{c} floor+lfloorfrac{b}{c} floor\ =& f(a mod c,b mod c,c,n)+frac{n(n+1)}{2}lfloorfrac{a}{c} floor+(n+1)lfloorfrac{b}{c} floor end{aligned} ]

（2）若(a<c)且(b<c)
此时使用上面的式子化简就变得没有意义了，我们需要考虑点其它的东西

考虑我们所计算的式子的实际意义：一条方程为(ax+cy+b=0)的直线，在这条直线下方且满足(xgeq 0),(y>0)的整点的个数，于是我们可以枚举(y)来进行化简
记(m=lfloor frac{an+b}{c} floor)

[egin{aligned} f(a,b,c,n)=&sum_{i=0}^nsum_{j=1}^m[lfloorfrac{ai+b}{c}geq j]\ =&sum_{i=0}^n sum_{j=0}^{m-1}[lfloorfrac{ai+b}{c} floorgeq j+1] end{aligned} ]

注意到下取整遇到大于号时可以直接丢掉，化简一下逻辑判断式

[egin{aligned} ai+b&geq jc+c\ i&geqfrac{jc+c-b}{a}\ i&>lfloor frac{jc+c-b-1}{a} floor end{aligned} ]

带回原式

[egin{aligned} f(a,b,c,n)=&sum_{i=0}^nsum_{j=0}^{m-1}[i>lfloorfrac{jc+c-b-1}{a} floor]\ =&sum_{j=0}^{m-1}sum_{i=0}^n[i>lfloorfrac{jc+c-b-1}{a} floor]\ =&sum_{j=0}^{m-1}(n-lfloorfrac{jc+c-b-1}{a} floor)\ =&n*m-f(c,c-b-1,a,m-1) end{aligned} ]

这样我们就得到了两个十分重要的递推式，注意到它的形式类似于欧几里得算法求(gcd)，所以被称作类欧几里得算法，时间复杂度和欧几里得算法几乎一致

接下来就是边界条件，这里的边界条件我写的比较多

（1）若(a=0)，则(f(a,b,c,n)=(n+1)lfloorfrac{b}{c} floor)

（2）若(n=0)，则(f(a,b,c,n)=lfloorfrac{b}{c} floor)

（3）若(n=1)，则(f(a,b,c,n)=lfloorfrac{a+b}{c} floor+lfloorfrac{b}{c} floor)

（4）若(c<0)，则(f(a,b,c,n)=f(-a,-b,-c,n))

（5）若(a<0)或(b<0)，利用最上面的第一条式子将(a,b)取正即可，即(f(a,b,c,n)=f(a mod c+c,b mod c+c,c,n)+frac{n(n+1)}{2}(lfloorfrac{a}{c} floor-1)+(n+1)(lfloorfrac{b}{c} floor-1))

反正最后一般就是用到(a=0)的情况了

ll work(ll a,ll b,ll c,ll n)
{
    if (!a) return (n+1)*(b/c);
    if (!n) return b/c;
    if (n==1) return (a+b)/c+b/c;
    if (c<0) return work(-a,-b,-c,n);
    ll d=abs(gcd(gcd(a,b),c));
    a/=d;b/=d;c/=d;
    if ((a>=c) || (b>=c)) return work(a%c,b%c,c,n)+n*(n+1)/2*(a/c)+(n+1)*(b/c);
    if ((a<0) || (b<0)) return work(a%c+c,b%c+c,c,n)+n*(n+1)/2*(a/c-1)+(n+1)*(b/c-1);
    ll m=(a*n+b)/c;
    return n*m-work(c,c-b-1,a,m-1);
}

扩展版

扩展版主要是两个式子

[g(a,b,c,n)=sum_{i=0}^n ilfloorfrac{ai+b}{c} floor ]

[h(a,b,c,n)=sum_{i=0}^nlfloorfrac{ai+b}{c} floor^2 ]

这两个式子的求解是互相关联的，这一点将在下面的推导过程中体现出来，并且它们的推导和(f)有很大的相似之处

求(g)

(g(a,b,c,n))较为简单，我们先推导它

（1）若(ageq c)或(bgeq c)
和上面一样的化简，注意到此时(lfloorfrac{ai}{c} floor)的前面还有一个系数(i)，因此按照上面的展开的话得到的系数是一个平方和，同样的(lfloorfrac{b}{c} floor)的系数也由常数(1)变成了等差数列求和，这里就直接给出结果

[g(a,b,c,n)=g(a mod c,b mod c,c,n)+frac{n(n+1)(2n+1)}{6}lfloorfrac{a}{c} floor+frac{n(n+1)}{2}lfloorfrac{b}{c} floor ]

（2） 若(a<c)且(b<c)

继续将(lfloorfrac{ai+b}{c} floor)展开

[egin{aligned} g(a,b,c,n)=&sum_{i=0}^nisum_{j=1}^m[lfloorfrac{ai+b}{c} floorgeq j]\ =&sum_{i=0}^nisum_{j=0}^{m-1}[lfloorfrac{ai+b}{c} floorgeq j+1]\ =&sum_{i=0}^nisum_{j=0}^{m-1}[i>frac{jc+c-b-1}{a}]\ =&sum_{j=0}^{m-1}sum_{i=0}^ni[i>frac{jc+c-b-1}{a}] end{aligned} ]

为了满足最后的逻辑表达式，(i)的取值在([lfloorfrac{jc+c-b-1}{a} floor+1,n])之中，所以这也是个等差数列，利用求和公式继续化简

[egin{aligned} g(a,b,c,n)=&sum_{j=0}^{m-1}frac{(n+1+lfloorfrac{jc+c-b-1}{a} floor)(n-lfloorfrac{jc+c-b-1}{a} floor)}{2}\ =&frac{1}{2}sum_{j=0}^{m-1}((n+1)+lfloorfrac{jc+c-b-1}{a} floor)(n-frac{jc+c-b-1}{a} floor)\ =&frac{1}{2}(mn(n+1)-sum_{j=0}^{m-1}(lfloorfrac{jc+c-b-1}{a} floor+lfloorfrac{jc+c-b-1}{a} floor^2))\ =&frac{1}{2}(mn(n+1)-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1)) end{aligned} ]

好了我们做完了

个鬼啊，那个(h)是什么啊

求(h)

依然是按照上面的思路分类讨论

（1）若(ageq c)或(bgeq c)
由((a+b+c)^2=a^2+b^2+c^2+2ab+2ac+2bc)可得

[egin{aligned} h(a,b,c,n)=&sum_{i=0}^n(lfloorfrac{a mod c*i+b mod c}{c} floor+ilfloorfrac{a}{c} floor+lfloorfrac{b}{c} floor)^2\ =&h(a mod c,b mod c,c,n)+frac{n(n+1)(2n+1)}{6}lfloorfrac{a}{c} floor^2+(n+1)lfloorfrac{b}{c} floor^2+2lfloorfrac{a}{c} floor g(a mod c,b mod c,c,n)+2lfloorfrac{b}{c} floor f(a mod c,b mod c,c,n)+n(n+1)lfloorfrac{a}{c} floorlfloorfrac{b}{c} floor end{aligned} ]

（2）若(a<c)且(b<c)

首先给出一个变形

[n^2=n(n+1)-n=2(sum_{i=1}^ni)-n ]

于是

[egin{aligned} h(a,b,c,d)=&sum_{i=0}^n(2sum_{j=1}^{lfloorfrac{ai+b}{c} floor}j-lfloorfrac{ai+b}{c} floor)\ =&2sum_{j=0}^{m-1}(j+1)sum_{i=0}^n[lfloorfrac{ai+b}{c} floorgeq j+1]-f(a,b,c,n)\ =&2sum_{j=0}^{m-1}(j+1)sum_{i=0}^n[igeq lfloorfrac{jc+c-b-1}{a} floor]-f(a,b,c,n)\ =&2sum_{j=0}^{m-1}(j+1)(n-lfloorfrac{jc+c-b-1}{a} floor)-f(a,b,c,d)\ =&2nsum_{j=0}^{m-1}(j+1)-2sum_{j=0}^{m-1}jlfloorfrac{jc+c-b-1}{a} floor-2sum_{j=0}^{m-1}lfloorfrac{jc+c-b-1}{a} floor-f(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) end{aligned} ]

好了这回真的做完了

注意到我们需要在(f,g,h)之间交替求值，我们对每一对((a,b,c,n))存下它的(f,g,h)值即可

node solve(ll a,ll b,ll c,ll n)
{
    node ans;
    if (!a)
    {
        ans.f=(n+1)*(b/c)%maxd;
        ans.g=(n+1)*n%maxd*inv2%maxd*(b/c)%maxd;
        ans.h=(n+1)*(b/c)%maxd*(b/c)%maxd;
        return ans;
    }
    if ((a>=c) || (b>=c))
    {
        node tmp=solve(a%c,b%c,c,n);
        ans.f=(tmp.f+n*(n+1)%maxd*inv2%maxd*(a/c)%maxd+(n+1)*(b/c)%maxd)%maxd;
        ans.g=(tmp.g+n*(n+1)%maxd*(n*2+1)%maxd*inv6%maxd*(a/c)%maxd+n*(n+1)%maxd*inv2%maxd*(b/c)%maxd)%maxd;
        ans.h=(tmp.h+n*(n+1)%maxd*(n*2+1)%maxd*inv6%maxd*(a/c)%maxd*(a/c)%maxd+(n+1)%maxd*(b/c)%maxd*(b/c)%maxd+tmp.f*(b/c)%maxd*2%maxd+tmp.g*(a/c)%maxd*2%maxd+(a/c)*(b/c)%maxd*n%maxd*(n+1)%maxd)%maxd;
        return ans;
    }
    ll m=(a*n+b)/c;
    node tmp=solve(c,c-b-1,a,m-1);m%=maxd;
    ans.f=(n*m%maxd+maxd-tmp.f)%maxd;
    ans.g=(m*n%maxd*(n+1)%maxd+maxd-tmp.f+maxd-tmp.h)%maxd*inv2%maxd;
    ans.h=(n*m%maxd*(m+1)%maxd-tmp.g*2%maxd-tmp.f*2%maxd-ans.f)%maxd;
    return ans;
}

记得在最后ans.f=(ans.f+maxd)%maxd;ans.g=(ans.g+maxd)%maxd;ans.h=(ans.h+maxd)%maxd; 强制转为正数即可