类欧几里得算法的推导

一般形式

$f(a,b,c,n)=\sum_{i=0}^n \lfloor \frac{ai+b}{c}\rfloor$
给定a,b,c,n，求f(a,b,c,n)
扩展：
$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$

f

首先从最简单的f开始推导。

$f(a,b,c,n)=\sum_{i=0}^n \lfloor \frac{ai+b}{c}\rfloor$

当a≥c或b≥c时
其实就相当于把上式的分数变成 $\lfloor \frac{ai}{c}+\frac{b}{c}\rfloor$，然后两个变成真分数。得到：
$f(a,b,c,n)=f(a$%$c,b$%$c,c,n)+\frac{n(n+1)}{2}\lfloor \frac{a}{c}\rfloor+(n+1)\lfloor\frac{b}{c}\rfloor$

当a < c且b < c时
令$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≥j]$
这个可以理解为一条直线在横坐标为0——n的整数时，越过第一象限的点数量。
继续推：
$f(a,b,c,n)=\sum_{i=0}^n \sum_{j=0}^{m-1}[\lfloor \frac{ai+b}{c}\rfloor≥j+1]$
把分数去掉
$f(a,b,c,n)=\sum_{i=0}^n \sum_{j=0}^{m-1}[ai≥jc+c-b]$
$f(a,b,c,n)=\sum_{i=0}^n \sum_{j=0}^{m-1}[ai>jc+c-b-1]$
两边除掉a
$f(a,b,c,n)=\sum_{i=0}^n \sum_{j=0}^{m-1}[i>\frac{jc+c-b-1}{a}]$
左边只有i，这样就舒服多了。下一步：
$f(a,b,c,n)=\sum_{j=0}^{m-1}(n-\lfloor\frac{jc+c-b-1}{a}\rfloor)$
观察到括号里面是一个等差数列和一个f形式的求和，得到：
$f(a,b,c,n)=nm-f(c,c-b-1,a,m-1)$
可以发现，第1、3个系数从(a,c)变成(c,a%c)。所以就叫它类欧几里得。
边界条件：a=0

g

$g(a,b,c,n)=\sum_{i=0}^n i\lfloor \frac{ai+b}{c}\rfloor$
首先依然是考虑a≥c或b≥c的情况。因为抽出来有一个平方和，所以得到：
$g(a,b,c,n)=g(a$%$c,b$%$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$

当a < c且b < c时，根据f的思路，可以得到：
$g(a,b,c,n)=\sum_{i=0}^ni \sum_{j=1}^{m}[\lfloor \frac{ai+b}{c}\rfloor≥j]$
$g(a,b,c,n)=\sum_{i=0}^n i\sum_{j=0}^{m-1}[i>\frac{jc+c-b-1}{a}]$
对于任意一个j，为了使表达式为真，i至少为$\lfloor\frac{jc+c-b-1}{a}\rfloor+1$，又因为上限为n，可以得到一个等差数列，所以：
$g(a,b,c,n)=\sum_{j=0}^{m-1}\frac{(n-\lfloor\frac{jc+c-b-1}{a}\rfloor)(n+\lfloor\frac{jc+c-b-1}{a}\rfloor+1)}{2}$
拆开来得到：
$g(a,b,c,n)=\frac{1}{2}\sum_{j=0}^{m-1}(n(n+1)-\lfloor\frac{jc+c-b-1}{a}\rfloor-\lfloor\frac{jc+c-b-1}{a}\rfloor^2)$
最终得到：
$g(a,b,c,n)=\frac{mn(n+1)-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1)}{2}$

式子还和h有关，所以还要推h

h

$h(a,b,c,n)=\sum_{i=0}^n \lfloor \frac{ai+b}{c}\rfloor^2$
当a≥c或b≥c，拆开来得到：
$h(a,b,c,n)=h(a$%$c,b$%$c,c,n)+\frac{n(n+1)(2n+1)}{6}\lfloor\frac{a}{c}\rfloor^2+(n+1)\lfloor\frac{b}{c}\rfloor^2+2\lfloor\frac{b}{c}\rfloor f(a$%$c,b$%$c,c,n)+2\lfloor\frac{a}{c}\rfloor g(a$%$c,b$%$c,c,n)+\lfloor\frac{a}{c}\rfloor\lfloor\frac{b}{c}\rfloor n(n+1)$

当a < c且b < c时
为了方便，要把$n^2$拆一下，得到：$n^2=2*\frac{n(n+1)}{2}-n=2\sum_{i=0}^ni-n$
可以得到：
$h(a,b,c,n)=\sum_{i=0}^n 2\sum_{j=1}^{\frac{ai+b}{c}}j-{\frac{ai+b}{c}}$
把两部分分开，得到：
$h(a,b,c,n)=2\sum_{j=0}^{m-1}(j+1)\sum_{i=0}^n [\lfloor\frac{ai+b}{c}\rfloor≥j+1]-f(a,b,c,n)$
中间的按照上面的思路推吧，最终得到：
$h(a,b,c,n)=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)$

一道求g裸题的代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const int mo=1e9+7,inv2=500000004,inv6=166666668;

typedef long long LL;

int a,b,c,l,r;

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

data calc(int a,int b,int c,LL n)
{
    data tmp;
    if (!a)
    {
        tmp.f=tmp.g=tmp.h=0;
        return tmp;
    }
    if (a>=c || b>=c)
    {
        tmp=calc(a%c,b%c,c,n);
        n%=mo;
        tmp.h=(tmp.h+n*(n+1)%mo*(2*n+1)%mo*inv6%mo*(a/c)%mo*(a/c)%mo
        +(n+1)*(b/c)%mo*(b/c)%mo
        +(LL)2*(a/c)*tmp.g%mo
        +(LL)2*(b/c)*tmp.f%mo
        +n*(n+1)%mo*(a/c)%mo*(b/c))%mo;
        tmp.f=(tmp.f+n*(n+1)/2%mo*(a/c)+(n+1)*(b/c))%mo;
        tmp.g=(tmp.g+n*(n+1)%mo*(2*n+1)%mo*inv6%mo*(a/c)+n*(n+1)/2%mo*(b/c))%mo;
        return tmp;
    }
    LL m=((LL)a*n+b)/c;
    data nxt=calc(c,c-b-1,a,m-1);
    n%=mo; m%=mo;
    tmp.f=((n*m-nxt.f)%mo+mo)%mo;
    tmp.g=(LL)((n*(n+1)%mo*m-nxt.f-nxt.h)%mo+mo)*inv2%mo;
    tmp.h=((m*(m+1)%mo*n-(LL)2*(nxt.g+nxt.f)%mo-tmp.f)%mo+mo)%mo;
    return tmp;
}

int main()
{
    freopen("task.in","r",stdin); freopen("task.out","w",stdout);
    scanf("%d%d%d%d%d",&a,&c,&b,&l,&r);
    printf("%d\n",(calc(a,b,c,r).g-calc(a,b,c,l-1).g+mo)%mo);
    return 0;
}