类欧几里得算法
前几日做了名为"It's a Mod, Mod, Mod, Mod World"(膜世界)的题
给定\(t\)组\(p,n,q\),求解
\(\sum_{i=1}^n [(pi)\,mod\,q]\)，其中\(1 \le p,n,q \le 10^6\),\(1 \le t \le 10^4\),限时\(1s\)
注意，这个\(mod\)是求和号里面的，不能对求和整体取模
这个题简直是开幕雷击，完全不知道怎么下手，有兴趣的读者可以思考半个小时。后来我百度了题目名，才知道要用到“类欧几里得算法”。
欧几里德算法又称辗转相除法，用于计算最大公因数（我们把a,b的最大公因数简称\(gcd(a,b)\)），计算\(gcd(a,b)\)的复杂度为\(O( log(max(a,b)) )\)，非常快。
同样的，在“类欧几里得算法”当中，也用到了类似辗转相除的操作，那么最终的复杂度也是log级别的，解决上述问题非常合适
最简单的"类欧几里得算法"可以用于求解下列函数的值
\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{ai+b}{c} \rfloor\)
其中，a,b,c,n都是给定的整数,时间复杂度为\(O( log(max(a,c)) )\)
下面给出算法的原理以及步骤:
首先，把这个式子的参数化小一点
\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{ai+b}{c}\rfloor
=\sum_{i=1}^n \lfloor\frac{c\lfloor\frac{a}{c}\rfloor i+(a\,mod\,c)i+b\,mod\,c+c\lfloor \frac{b}{c} \rfloor}{c} \rfloor
=\sum_{i=1}^n \lfloor\frac{(a\,mod\,c)i+b\,mod\,c}{c} \rfloor+\sum_{i=1}^n \lfloor \frac{a}{c} \rfloor i+ \lfloor \frac{b}{c} \rfloor\)
第三个等号的依据是
当\(y\)是一个整数，\(x\)是一个实数，则有
\(\lfloor x+y \rfloor =\lfloor x \rfloor + y\)
其中，第二个求和式可以用等差数列公式计算，结果为
\(\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\)
用\(f\)的形式简写，可以得到
等式1：\(f(a,b,c,n)=f(a\,mod\,c,b\,mod\,c,c,n)+\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\)
以上操作让后得到的\(a\,mod\,c\),\(b\,mod\,c\)均小于c，所以在下面的讨论中，不妨假设\(a,b<c\)，这是因为只要我们能求出满足\(a,b<c\)的\(f\)函数的值，就能求出所有的\(f\)函数的值
下面分两种情况研究\(f\)函数
1.特殊情况:
\(a=0\)
\(f(a,b,c,n)=\sum_{i=1}^n \lfloor\frac{b}{c} \rfloor=\lfloor\frac{b}{c} \rfloor n\)
这种情况对应着递归求解的终点，在下面的一般情况的推导中，会介绍递归算法的递推方程
2.一般情况：
\(a>0\)
\(\sum_{i=1}^n \lfloor \frac{ai+b}{c} \rfloor\)
\(=\sum_{i=1}^n \sum_{j=1}^{\lfloor\frac{ai+b}{c}\rfloor} 1\)
\(=\sum_{i=1}^n \sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} (j\le \lfloor\frac{ai+b}{c}\rfloor)\)
其中\(x\le y\)为逻辑表达式，逻辑表达式为真时值为1，为假时值为0
第二层求和的上限从\(\lfloor\frac{ai+b}{c}\rfloor\)变成了\(\lfloor\frac{an+b}{c}\rfloor\)
这是为了让两层求和上限不相关，可以交换两个求和号的顺序
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} \sum_{i=1}^n (j\le \lfloor\frac{ai+b}{c}\rfloor)\)
下面，我们只需要变换\(j\le \lfloor\frac{ai+b}{c}\rfloor\)，使求和式容易计算
注意：我们对对该不等式只能进行等价变换，非等价变换会导致多算或漏算，切记
\(j\le \lfloor\frac{ai+b}{c}\rfloor\)
\(\Leftrightarrow cj\le ai+b\)
依据：当x,y为整数,z为正整数,\(x\le \lfloor \frac{y}{z} \rfloor\Leftrightarrow xz \le y\)
\(\Leftrightarrow ai \ge cj-b\)
\(\Leftrightarrow i \ge \frac{cj-b}{a}\)
\(\Leftrightarrow i \ge \lceil \frac{cj-b}{a}\rceil\)
依据：当x,y为整数,z为正整数,\(x\ge \frac{y}{z}\Leftrightarrow x\ge \lceil \frac{y}{z} \rceil\)
\(\Leftrightarrow i \ge \lceil \frac{cj-b}{a}\rceil\)
\(\Leftrightarrow i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
依据：当x为整数,y为正整数,\(\lceil \frac{x}{y}\rceil=\lfloor \frac{x+y-1}{y}\rfloor\)，注意：这个依据不是很显然，读者可以自己证明一下
所以，原求和式
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} \sum_{i=1}^n i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
对于\(i \ge \lfloor \frac{cj-b+a-1}{a}\rfloor\)
由\(a,b<c\) , \(1 \le j \le \lfloor\frac{an+b}{c}\rfloor\)可知
\(\lfloor \frac{cj-b+a-1}{a}\rfloor \ge \lfloor \frac{c-b+a-1}{a}\rfloor = 1+\lfloor \frac{c-b-1}{a}\rfloor \ge 1\)
\(\lfloor \frac{cj-b+a-1}{a}\rfloor \le \lfloor \frac{c\lfloor\frac{an+b}{c}\rfloor-b+a-1}{a}\rfloor \le \lfloor \frac{an+b-b+a-1}{a}\rfloor =n\)
即\(1 \le\lfloor \frac{cj-b+a-1}{a}\rfloor \le n\)
这一点非常重要，这确保当\(i\)可以取到\(\lfloor\frac{cj-b+a-1}{a}\rfloor \le i \le n\) 的所有取值。
读者可以想一想如果\(\lfloor \frac{cj-b+a-1}{a}\rfloor\)可能大于\(n\)或小于\(1\)时会发生什么
所以，原求和式
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} n-\lfloor \frac{cj-b+a-1}{a}\rfloor+1\)
\(=\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor} n-\lfloor \frac{cj-b-1}{a}\rfloor\)
\(=n(\lfloor\frac{an+b}{c}\rfloor)-\sum_{j=1}^{\lfloor\frac{an+b}{c}\rfloor}\lfloor \frac{cj-b-1}{a}\rfloor=n(\lfloor\frac{an+b}{c}\rfloor)-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)
即有等式2：\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)
这里可以看出左右两式\(c,a\)的位置互换了
因为\(c>a\)
所以可以用等式1处理\(f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\)
得\(f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)=f(c\,mod\,a,(-b-1)\,mod\,a,a,\lfloor\frac{an+b}{c}\rfloor)+\frac{1}{2} \lfloor \frac{c}{a} \rfloor \lfloor\frac{an+b}{c}\rfloor^2+\frac{1}{2}(\lfloor \frac{c}{a} \rfloor+2\lfloor\frac{-b-1}{a}\rfloor)\lfloor\frac{an+b}{c}\rfloor\)
即等式3：
\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n\frac{1}{2} \lfloor \frac{c}{a} \rfloor \lfloor\frac{an+b}{c}\rfloor^2-\frac{1}{2}(\lfloor \frac{c}{a} \rfloor+2\lfloor\frac{-b-1}{a}\rfloor)\lfloor\frac{an+b}{c}\rfloor-f(c\,mod\,a,(-b-1)\,mod\,a,a,\lfloor\frac{an+b}{c}\rfloor)\)
这一等式中
左式中\(f\)函数第一个参数和第三个参数分别是\(a,c\),右式中\(f\)的第一个参数和第三个参数分别是\(c\,mod\,a,a\)
而欧几里德算法（辗转相除法）在运算时也是利用
\(gcd(a,b)=gcd(b\,mod\,a,a)\)（假设第一个参数<第二个参数）
也就是说在递归求解过程中，类欧几里德算法的第1,3参数变化和欧几里得算法的两个参数变化的方式完全一致
最终的参数状态就是第一个参数变成0，对应着“特殊情况”，标志着递归的结束
故求解时间复杂度和欧几里德算法一致，是\(O(log(max(a,c)))\)的
最后总结一下递推方程：
\(f(0,b,c,n)=\lfloor\frac{b}{c} \rfloor n\)
\(a \ge c\)时
\(f(a,b,c,n)=f(a\,mod\,c,b\,mod\,c,c,n)+\frac{1}{2} \lfloor \frac{a}{c} \rfloor n^2+\frac{1}{2}(\lfloor \frac{a}{c} \rfloor+2\lfloor\frac{b}{c}\rfloor)n\) （等式1）
\(a<c\)时
\(f(a,b,c,n)=(\lfloor\frac{an+b}{c}\rfloor)n-f(c,-b-1,a,\lfloor\frac{an+b}{c}\rfloor)\) （等式2）
由于等式3是由等式1和等式2通过代换推出的，递归时没有必要进行等式3的计算。
这道题要求
\(\sum_{i=1}^n [(pi)\,mod\,q]\)
即\(\sum_{i=1}^npi-q\sum_{i=1}^n \lfloor \frac{pi}{q} \rfloor\)
\(=p\frac{n(n+1)}{2}-qf(p,0,q,n)\)
递归求解\(f(p,0,q,n)\)即可算出答案
特别提醒：在编程时，负数参与的整除和模运算会有一些问题
比如对于计算机来说，\(-5/3\)的值为\(-1\)，而\(\lfloor\frac{-5}{3}\rfloor=-2\)
模的情况也是比较复杂，读者可以自行测试
