前几日做了名为"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 lfloorfrac{ai+b}{c}
floor)
其中，a,b,c,n都是给定的整数,时间复杂度为(O( log(max(a,c)) ))
下面给出算法的原理以及步骤:
首先，把这个式子的参数化小一点
(f(a,b,c,n)=sum_{i=1}^n lfloorfrac{ai+b}{c}
floor
=sum_{i=1}^n lfloorfrac{clfloorfrac{a}{c}
floor i+(a\,mod\,c)i+b\,mod\,c+clfloor frac{b}{c}
floor}{c}
floor
=sum_{i=1}^n lfloorfrac{(a\,mod\,c)i+b\,mod\,c}{c}
floor+sum_{i=1}^n lfloor frac{a}{c}
floor i+ lfloor frac{b}{c}
floor)
第三个等号的依据是
当(y)是一个整数，(x)是一个实数，则有
(lfloor x+y
floor =lfloor x
floor + y)
其中，第二个求和式可以用等差数列公式计算，结果为
(frac{1}{2} lfloor frac{a}{c}
floor n^2+frac{1}{2}(lfloor frac{a}{c}
floor+2lfloorfrac{b}{c}
floor)n)
用(f)的形式简写，可以得到
等式1：(f(a,b,c,n)=f(a\,mod\,c,b\,mod\,c,c,n)+frac{1}{2} lfloor frac{a}{c}
floor n^2+frac{1}{2}(lfloor frac{a}{c}
floor+2lfloorfrac{b}{c}
floor)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 lfloorfrac{b}{c} floor=lfloorfrac{b}{c} floor n)
这种情况对应着递归求解的终点，在下面的一般情况的推导中，会介绍递归算法的递推方程
2.一般情况：
(a>0)
(sum_{i=1}^n lfloor frac{ai+b}{c}
floor)
(=sum_{i=1}^n sum_{j=1}^{lfloorfrac{ai+b}{c}
floor} 1)
(=sum_{i=1}^n sum_{j=1}^{lfloorfrac{an+b}{c}
floor} (jle lfloorfrac{ai+b}{c}
floor))
其中(xle y)为逻辑表达式，逻辑表达式为真时值为1，为假时值为0
第二层求和的上限从(lfloorfrac{ai+b}{c}
floor)变成了(lfloorfrac{an+b}{c}
floor)
这是为了让两层求和上限不相关，可以交换两个求和号的顺序
(=sum_{j=1}^{lfloorfrac{an+b}{c}
floor} sum_{i=1}^n (jle lfloorfrac{ai+b}{c}
floor))
下面，我们只需要变换(jle lfloorfrac{ai+b}{c}
floor)，使求和式容易计算
注意：我们对对该不等式只能进行等价变换，非等价变换会导致多算或漏算，切记
(jle lfloorfrac{ai+b}{c}
floor)
(Leftrightarrow cjle ai+b)
依据：当x,y为整数,z为正整数,(xle lfloor frac{y}{z}
floorLeftrightarrow xz le y)
(Leftrightarrow ai ge cj-b)
(Leftrightarrow i ge frac{cj-b}{a})
(Leftrightarrow i ge lceil frac{cj-b}{a}
ceil)
依据：当x,y为整数,z为正整数,(xge frac{y}{z}Leftrightarrow xge lceil frac{y}{z}
ceil)
(Leftrightarrow i ge lceil frac{cj-b}{a}
ceil)
(Leftrightarrow i ge lfloor frac{cj-b+a-1}{a}
floor)
依据：当x为整数,y为正整数,(lceil frac{x}{y}
ceil=lfloor frac{x+y-1}{y}
floor)，注意：这个依据不是很显然，读者可以自己证明一下
所以，原求和式
(=sum_{j=1}^{lfloorfrac{an+b}{c}
floor} sum_{i=1}^n i ge lfloor frac{cj-b+a-1}{a}
floor)
对于(i ge lfloor frac{cj-b+a-1}{a}
floor)
由(a,b<c) , (1 le j le lfloorfrac{an+b}{c}
floor)可知
(lfloor frac{cj-b+a-1}{a}
floor ge lfloor frac{c-b+a-1}{a}
floor = 1+lfloor frac{c-b-1}{a}
floor ge 1)
(lfloor frac{cj-b+a-1}{a}
floor le lfloor frac{clfloorfrac{an+b}{c}
floor-b+a-1}{a}
floor le lfloor frac{an+b-b+a-1}{a}
floor =n)
即(1 lelfloor frac{cj-b+a-1}{a} floor le n)
这一点非常重要，这确保当(i)可以取到(lfloorfrac{cj-b+a-1}{a} floor le i le n) 的所有取值。
读者可以想一想如果(lfloor frac{cj-b+a-1}{a} floor)可能大于(n)或小于(1)时会发生什么
所以，原求和式
(=sum_{j=1}^{lfloorfrac{an+b}{c}
floor} n-lfloor frac{cj-b+a-1}{a}
floor+1)
(=sum_{j=1}^{lfloorfrac{an+b}{c} floor} n-lfloor frac{cj-b-1}{a} floor)
(=n(lfloorfrac{an+b}{c} floor)-sum_{j=1}^{lfloorfrac{an+b}{c} floor}lfloor frac{cj-b-1}{a} floor=n(lfloorfrac{an+b}{c} floor)-f(c,-b-1,a,lfloorfrac{an+b}{c} floor))
即有等式2：(f(a,b,c,n)=(lfloorfrac{an+b}{c} floor)n-f(c,-b-1,a,lfloorfrac{an+b}{c} floor))
这里可以看出左右两式(c,a)的位置互换了
因为(c>a)
所以可以用等式1处理(f(c,-b-1,a,lfloorfrac{an+b}{c} floor))
得(f(c,-b-1,a,lfloorfrac{an+b}{c} floor)=f(c\,mod\,a,(-b-1)\,mod\,a,a,lfloorfrac{an+b}{c} floor)+frac{1}{2} lfloor frac{c}{a} floor lfloorfrac{an+b}{c} floor^2+frac{1}{2}(lfloor frac{c}{a} floor+2lfloorfrac{-b-1}{a} floor)lfloorfrac{an+b}{c} floor)
即等式3：
(f(a,b,c,n)=(lfloorfrac{an+b}{c} floor)nfrac{1}{2} lfloor frac{c}{a} floor lfloorfrac{an+b}{c} floor^2-frac{1}{2}(lfloor frac{c}{a} floor+2lfloorfrac{-b-1}{a} floor)lfloorfrac{an+b}{c} floor-f(c\,mod\,a,(-b-1)\,mod\,a,a,lfloorfrac{an+b}{c} floor))
这一等式中
左式中(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)=lfloorfrac{b}{c} floor 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} floor n^2+frac{1}{2}(lfloor frac{a}{c} floor+2lfloorfrac{b}{c} floor)n) （等式1）
(a<c)时
(f(a,b,c,n)=(lfloorfrac{an+b}{c} floor)n-f(c,-b-1,a,lfloorfrac{an+b}{c} floor)) （等式2）
由于等式3是由等式1和等式2通过代换推出的，递归时没有必要进行等式3的计算。
这道题要求
(sum_{i=1}^n [(pi)\,mod\,q])
即(sum_{i=1}^npi-qsum_{i=1}^n lfloor frac{pi}{q} floor)
(=pfrac{n(n+1)}{2}-qf(p,0,q,n))
递归求解(f(p,0,q,n))即可算出答案
特别提醒：在编程时，负数参与的整除和模运算会有一些问题
比如对于计算机来说，(-5/3)的值为(-1)，而(lfloorfrac{-5}{3} floor=-2)
模的情况也是比较复杂，读者可以自行测试