类欧几里得算法

类欧几里得算法是用 $O(\log n)$ 时间求形如 $\sum _{i=0}^n\lfloor \frac{ai+b}{c}\rfloor$ 的式子的答案的。$a,c>0,b\ge 0$

类欧几里得算法

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

我们可以把 $a,b$ 对 $c$ 取模。
$$\begin{array}{rl}f(n,a,b,c)& =\sum _{i=0}^n\lfloor \frac{(c\lfloor \frac{a}{c}\rfloor +a\mathrm{mod}c)i+c\lfloor \frac{b}{c}\rfloor +b\mathrm{mod}c}{c}\rfloor \\ & =\sum _{i=0}^n\lfloor \frac{a}{c}\rfloor i+\lfloor \frac{b}{c}\rfloor +\lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor \\ & =\lfloor \frac{a}{c}\rfloor \frac{n(n+1)}{2}+(n+1)\lfloor \frac{b}{c}\rfloor +f(n,a\mathrm{mod}c,b\mathrm{mod}c,c)\end{array}$$

下面的 $a,b$ 都是模过 $c$ 的。

考虑继续化简变换。设 $m=\lfloor \frac{an+b}{c}\rfloor$

$$\begin{array}{rl}f(n,a,b,c)& =\sum _{i=0}^n\lfloor \frac{ai+b}{c}\rfloor \\ & =\sum _{i=0}^n\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}1\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j+1\le \lfloor \frac{ai+b}{c}\rfloor \right]\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j+1\le \frac{ai+b}{c}\right]\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[cj+c-b-1<ai\right]\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[\lfloor \frac{cj+c-b-1}{a}\rfloor <i\right]\\ & =\sum _{j=0}^{m-1}\left(n-\lfloor \frac{cj+c-b-1}{a}\rfloor \right)\\ & =nm-f(m-1,c,c-b-1,a)\end{array}$$

这样，我们就得到一个递推式。发现 $a,c$ 位置调换了，而且 $a$ 还模了 $c$。这类似于辗转相除法，故时间复杂度是 $O(\log n)$ 的。

其实这个算法只有时间复杂度和欧几里得扯得上关系。

下面是参考实现代码

ll f(ll n,ll a,ll b,ll c)
{
    if(!a) return (n+1)*(b/c);
    if(a>=c||b>=c) return f(n,a%c,b%c,c)+n*(n+1)/2*(a/c)+(n+1)*(b/c);
    return (a*n+b)/c*n-f((a*n+b)/c-1,c,c-b-1,a);
}

拓展

求 $\sum _{i=0}^{n}i\lfloor \frac{ai+b}{c}\rfloor ,\sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2$

设 $g(n,a,b,c)=\sum _{i=0}^{n}i\lfloor \frac{ai+b}{c}\rfloor ,h(n,a,b,c)=\sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2$

求 $g$：首先 $a,b$ 对 $c$ 取模。类似得到
$$\begin{array}{rl}g(n,a,b,c)& =\sum _{i=0}^{n}i\lfloor \frac{(c\lfloor \frac{a}{c}\rfloor +a\mathrm{mod}c)i+c\lfloor \frac{b}{c}\rfloor +b\mathrm{mod}c}{c}\rfloor \\ & =\sum _{i=0}^n\left(i^2\lfloor \frac{a}{c}\rfloor +i\lfloor \frac{b}{c}\rfloor +i\lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor \right)\\ & =\lfloor \frac{a}{c}\rfloor \frac{n(n+1)(2n+1)}{6}+\lfloor \frac{b}{c}\rfloor \frac{n(n+1)}{2}+g(n,a\mathrm{mod}c,b\mathrm{mod}c,c)\end{array}$$

设 $p=\lfloor \frac{cj+c-b-1}{a}\rfloor$

$$\begin{array}{rl}g(n,a,b,c)& =\sum _{i=0}^{n}i\lfloor \frac{ai+b}{c}\rfloor \\ & =\sum _{i=0}^n\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}i\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]i\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[j+1\le \frac{ai+b}{c}\right]i\\ & =\sum _{j=0}^{m-1}\sum _{i=0}^n\left[p<i\right]i\\ & =\sum _{j=0}^{m-1}\frac{(n+p+1)(n-p)}{2}\\ & =\frac{1}{2}\sum _{j=0}^{m-1}\left(n^2+n-p^2-p\right)\\ & =\frac{1}{2}\left(nm(n+1)-\sum _{j=0}^{m-1}{p}^2-\sum _{j=0}^{m-1}p\right)\\ & =\frac{1}{2}\left(nm(n+1)-h(m-1,c,c-b-1,a)-f(m-1,c,c-b-1,a)\right)\end{array}$$

发现这里 $g$ 函数的递归调用了 $h$，接下来求 $h$。$a,b$ 还是要对 $c$ 取模。

$$\begin{array}{rl}h(n,a,b,c)& =\sum _{i=0}^n{\lfloor \frac{(c\lfloor \frac{a}{c}\rfloor +a\mathrm{mod}c)i+c\lfloor \frac{b}{c}\rfloor +b\mathrm{mod}c}{c}\rfloor }^2\\ & =\sum _{i=0}^n{\left(i\lfloor \frac{a}{c}\rfloor +\lfloor \frac{b}{c}\rfloor +\lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor \right)}^2\\ & =\sum _{i=0}^n\left(i^2{\lfloor \frac{a}{c}\rfloor }^2+{\lfloor \frac{b}{c}\rfloor }^2+{\lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor }^2+2i\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor +2i\lfloor \frac{a}{c}\rfloor \lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor +2\lfloor \frac{b}{c}\rfloor \lfloor \frac{(a\mathrm{mod}c)i+b\mathrm{mod}c}{c}\rfloor \right)\\ & ={\lfloor \frac{a}{c}\rfloor }^2\frac{n(n+1)(2n+1)}{6}+{\lfloor \frac{b}{c}\rfloor }^2(n+1)+\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor n(n+1)+h(n,a\mathrm{mod}c,b\mathrm{mod}c,c)+2\lfloor \frac{a}{c}\rfloor g(n,a\mathrm{mod}c,b\mathrm{mod}c,c)+2\lfloor \frac{b}{c}\rfloor f(n,a\mathrm{mod}c,b\mathrm{mod}c,c)\end{array}$$

虽然挺复杂，但还是能表示。

继续求值，和上面差不多，但是有点不一样。

如果要把 ${\lfloor \frac{ai+b}{c}\rfloor }^2$ 像上面一样转化，就变成 $\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}1$，有两个求和符号，不好处理，考虑换一种表示方式。下面的 $m,p$ 和上面的是一样的。

由于 $n^2=2\cdot \frac{n(n+1)}{2}-n=2\sum _{i=1}^{n}i-n$

所以
$$\begin{array}{rl}h(n,a,b,c)& =\sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2\\ & =\sum _{i=0}^n\left(2\sum _{j=1}^{\lfloor \frac{ai+b}{c}\rfloor }j-\lfloor \frac{ai+b}{c}\rfloor \right)\\ & =2\sum _{i=0}^n\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}(j+1)-f(n,a,b,c)\end{array}$$

只看前面
$$\begin{array}{rl}& =2\sum _{i=0}^n\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor -1}(j+1)\\ & =2\sum _{j=0}^{m-1}(j+1)\sum _{i=0}^n\left[j<\lfloor \frac{ai+b}{c}\rfloor \right]\\ & =2\sum _{j=0}^{m-1}(j+1)\sum _{i=0}^n\left[p<i\right]\\ & =2\sum _{j=0}^{m-1}(j+1)(n-p)\\ & =2\sum _{j=0}^{m-1}\left(n(j+1)-pj-p\right)\\ & =nm(m+1)-2g(m-1,c,c-b-1,a)-2f(m-1,c,c-b-1,a)\end{array}$$

即 $h(n,a,b,c)=n{m}^2-2g(m-1,c,c-b-1,a)-f(m-1,c,c-b-1,a)$

$f,g,h$ 相互调用，但是里面参数变化是一致的。所以可以同时递归求 $f,g,h$，时间复杂度为 $O(\log n)$。

例题 P5170 【模板】类欧几里得算法

参考代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353,inv2=499122177,inv6=166374059;
ll n,a,b,c;
struct node
{
    ll f,g,h;
};
node solve(ll n,ll a,ll b,ll c)
{
    if(!a) return {b/c*(n+1)%mod,b/c*n%mod*(n+1)%mod*inv2%mod,(b/c)*(b/c)%mod*(n+1)%mod};
    if(a>=c||b>=c){
        node x=solve(n,a%c,b%c,c);
        x.h=((a/c)*(a/c)%mod*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod+(b/c)*(b/c)%mod*(n+1)%mod+(a/c)*(b/c)%mod*n%mod*(n+1)%mod+x.h+a/c*2*x.g+b/c*2*x.f)%mod;
        x.f=(x.f+a/c*n%mod*(n+1)%mod*inv2%mod+b/c*(n+1))%mod;
        x.g=(x.g+a/c*n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod+b/c*n%mod*(n+1)%mod*inv2)%mod;
        return x;
    }
    ll m=(a*n+b)/c;
    node x=solve(m-1,c,c-b-1,a),ans;
    ans.f=(n*m%mod-x.f+mod)%mod;
    ans.g=inv2*(n*m%mod*(n+1)%mod-x.h-x.f+2*mod)%mod;
    ans.h=(n*m%mod*m%mod-2*x.g-x.f+4*mod)%mod;
    return ans;
}
int main()
{
    int t;
    cin>>t;
    while(t--){
        scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
        node ans=solve(n,a,b,c);
        printf("%lld %lld %lld\n",ans.f,ans.h,ans.g);
    }
}