类欧几里得算法

前几日做了名为"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)

模的情况也是比较复杂，读者可以自行测试

注：笔者在写完这篇文章的一个月后做了一道相关题目，WA了10次，最后发现就是整除和模没有处理好，说明这一点很容易错。理论上整除和模都要自己手写，除非能保证某个算法绝对不会出现负数参与整除或模。