类欧几里得算法

基础

设

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

其中，\(a,b,c,n\) 为常数，如何 \(O(\log n)\) 求解。

当 \(a\geq c\) 或者 \(b\geq c\) ，可以将 \(a,b\) 对 \(c\) 取模以简化问题：

\[\begin{align} f(a,b,c,n) &= \sum_{i=0}^{n}{\left\lfloor \frac{ai+b}{c} \right\rfloor}\\ &=\sum_{i=0}^{n}{\left\lfloor \frac{(\lfloor \frac{a}{c}\rfloor c+a \bmod c)i+(\lfloor \frac{b}{c}\rfloor c+b \bmod c)}{c} \right\rfloor}\\ &=\frac{(1+n)n}{2}·\left\lfloor \frac{a}{c}\right\rfloor+(n+1)·\left\lfloor \frac{b}{c}\right\rfloor+\sum_{i=0}^{n}{ \left\lfloor \frac{(a \bmod c)i+(b \bmod c)}{c} \right\rfloor}\\ &=\frac{(1+n)n}{2}·\left\lfloor \frac{a}{c}\right\rfloor+(n+1)·\left\lfloor \frac{b}{c}\right\rfloor+f(a \bmod c ,b \bmod c,c,n)\\ \end{align} \]

那么，问题就转化为 \(a<c,b<c\) 的情况了。观察式子，发现只有 \(i\) 这个变量，需要从 \(i\) 下手推。在推求和式子中一个常见的技巧：条件和贡献的放缩与转化。具体而言，在原式中，\(0\leq i \leq n\) 是条件，而 \(\left\lfloor \frac{ai+b}{c} \right\rfloor\) 是对总和的贡献。

要加快一个求和式的计算过程，关键在于贡献合并计算。但这个式子的贡献并不好合并，因此考虑转化，得到另一个形式的和式：

\[f(a,b,c,n) = \sum_{i=0}^{n}{\left\lfloor \frac{ai+b}{c} \right\rfloor}=\sum_{i=0}^{n}{\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor-1}{1}} \]

现在增加了一个变量 \(j\) ，既然不好算 \(i\) 的贡献，就想办法计算 \(j\) 的贡献。为了得到一个和 \(j\) 有关的贡献式，采用数论中常用的变换方法：限制转移，即交换 \(i,j\) 的求和算子，用 \(n\) 限制 \(j\) 的上界。

\[\Rightarrow \sum_{j=0}^{\left\lfloor \frac{an+b}{c} \right\rfloor-1}{\sum_{i=0}^{n}{\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor \right]}} \]

这样做的目的是让 \(j\) 摆脱 \(i\) 的限制，现在 \(i,j\) 都被 \(n\) 限制。接着，对贡献式进行一些放缩处理。

先将取整符号拿掉：

\[j< \left\lfloor \frac{ai+b}{c} \right\rfloor \Leftrightarrow j+1 \leq \left\lfloor \frac{ai+b}{c} \right\rfloor \Leftrightarrow j+1 \leq \frac{ai+b}{c} \]

再做一些变换：

\[j+1 \leq \frac{ai+b}{c} \Leftrightarrow (j+1)c \leq ai+b \Leftrightarrow jc+c-b-1 < ai \]

最后一步，向下取整：

\[jc+c-b-1 < ai \Leftrightarrow \left\lfloor \frac{jc+c-b-1}{a} \right\rfloor <i \]

这一步的重要意义在于把变量 \(i\) 消掉了，令 \(m=\left\lfloor \frac{an+b}{c} \right\rfloor\) ，那么原式可以化为：

\[\begin{align} f(a,b,c,n) &= \sum_{j=0}^{m-1}{\sum_{i=0}^{n}{\left[ i>\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor \right]}}\\ &= \sum_{j=0}^{m-1}{\left(n-\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor \right)}\\ &= nm-\sum_{j=0}^{m-1}{\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor}\\ &= nm-f(c,c-b-1,a,m-1)\\ &= n\left\lfloor \frac{an+b}{c} \right\rfloor-f(c,c-b-1,a,\left\lfloor \frac{an+b}{c} \right\rfloor-1) \end{align} \]

可以发现，这是一个递归的式子。边界是：

• \(n=0\) 时，直接返回 \(\lfloor \frac{b}{c} \rfloor\)
• \(a=0\) 时，返回 \((n+1)\times\lfloor \frac{b}{c}\rfloor\)

综上所述：

\[f(a,b,c,n)= \begin{cases} \lfloor \frac{b}{c} \rfloor & n=0\\ (n+1)\times \lfloor \frac{b}{c} \rfloor & a=0\\ \frac{(1+n)·n}{2}\times \left\lfloor \frac{a}{c}\right\rfloor+(n+1)\times \left\lfloor \frac{b}{c}\right\rfloor+f(a \% c ,b\%c,c,n) & a\geq c\ or\ b\geq c\\ n\times \left\lfloor \frac{an+b}{c} \right\rfloor-f(c,c-b-1,a,\left\lfloor \frac{an+b}{c} \right\rfloor-1) & a<c\ and \ b<c \end{cases} \]

代码：

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

例题：

AtCoder Regular Contest 111 E - Simple Math 3

https://atcoder.jp/contests/arc111/tasks/arc111_e

拓展

求：

\[g(a,b,c,n)=\sum_{i=0}^{n}{i\left\lfloor \frac{ai+b}{c} \right\rfloor}\\ h(a,b,c,n)=\sum_{i=0}^{n}{\left\lfloor \frac{ai+b}{c} \right\rfloor^2} \]

g 函数求解

首先，考虑取模：

\[g(a,b,c,n)=g(a\%c,b\%c,c,n)+\frac{n(n+1)(2n+1)}{6} \left\lfloor \frac{a}{c} \right\rfloor+\frac{n(n+1)}{2}\left\lfloor \frac{b}{c} \right\rfloor \]

接下来，考虑 \(a<c,b<c\) 的情况，令 \(m=\left\lfloor \frac{an+b}{c} \right\rfloor\)，有：

\[g(a,b,c,n)=\sum_{i=0}^{n}{i\left\lfloor \frac{ai+b}{c} \right\rfloor}=\sum_{j=0}^{m-1}{\sum_{i=0}^{n}{i·\left[j<\left\lfloor \frac{ai+b}{c} \right\rfloor \right]}} \]

令 \(t=\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor\)，同理有：

\[\begin{align} &= \sum_{j=0}^{m-1}{\sum_{i=0}^{n}{[i>t]·i}}\\ &= \sum_{j=0}^{m-1}{\frac{1}{2}(t+1+n)(n-t)}\\ &= \frac{1}{2} \left(mn(n+1)-\sum_{j=0}^{m-1}{t^2}-\sum_{j=0}^{m-1}{t} \right)\\ &= \frac{1}{2} \left(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1) \right) \end{align} \]

综上所述：

\[g(a,b,c,n)=\begin{cases} 0 & n=0\\ \frac{n(n+1)}{2}\times \left\lfloor \frac{b}{c} \right\rfloor & a=0\\ g(a\%c,b\%c,c,n)+\frac{n(n+1)(2n+1)}{6} \left\lfloor \frac{a}{c} \right\rfloor+\frac{n(n+1)}{2}\left\lfloor \frac{b}{c} \right\rfloor & a\geq c\ or\ b\geq c \\ \frac{1}{2} \left(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1) \right) & a<c\ and\ b<c\\ \end{cases} \]

h 函数求解

首先取模：

\[\begin{align} h(a,b,c,n) &= h(a\%c,b\%c,c,n)\\ &+ 2 \left\lfloor \frac{b}{c} \right\rfloor f(a\%c,b\%c,c,n)+2 \left\lfloor \frac{a}{c} \right\rfloor g(a\%c,b\%c,c,n)\\ &+ \left\lfloor \frac{a}{c} \right\rfloor^2 \frac{n(n+1)(2n+1)}{6}+\left\lfloor \frac{b}{c} \right\rfloor^2 (n+1)+\left\lfloor \frac{a}{c} \right\rfloor \left\lfloor \frac{b}{c} \right\rfloor n(n+1) \end{align} \]

考虑 \(a<c,b<c\) 的情况，令 \(m=\left\lfloor \frac{ai+b}{c} \right\rfloor\)，\(t=\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor\)。

由于平方不好处理，可以按照下列方式拆分：

\[n^2=2·\frac{n(n+1)}{2}-n=\left(2\sum_{i=0}^{n}{i} \right)-n \]

这样做的意义在于，在添加变量 \(j\) 时只会变成一个求和算子，不会出现 \(\sum \times \sum\) 的形式。

\[\begin{align} h(a,b,c,n) &= \sum_{i=0}^{n}{\left\lfloor \frac{ai+b}{c} \right\rfloor^2}\\ &= \sum_{i=0}^{n}{\left[\left(2\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}{j} \right)-\left\lfloor \frac{ai+b}{c} \right\rfloor \right]}\\ &= \left(2\sum_{i=0}^{n}{\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}{j}} \right)-f(a,b,c,n)\\ \end{align} \]

接下来化简前一部分：

\[\begin{align} & \sum_{i=0}^{n}{\sum_{j=0}^{\left\lfloor \frac{ai+b}{c} \right\rfloor}{j}}\\ =& \sum_{i=0}^{n}{\sum_{j=0}^{\left\lfloor \frac{ai+b}{c}\right\rfloor-1}{(j+1)}}\\ =& \sum_{j=0}^{m-1}{(j+1)\sum_{i=0}^{n}{\left[j<\left\lfloor\frac{ai+b}{c}\right\rfloor \right]}}\\ =& \sum_{j=0}^{m-1}{(j+1)\sum_{i=0}^{n}{\left[i>t \right]}}\\ =& \sum_{j=0}^{m-1}{(j+1)(n-t)}\\ =& \frac{1}{2}nm(m+1)-\sum_{j=0}^{m-1}{(j+1)\left\lfloor \frac{jc+c-b-1}{a} \right\rfloor}\\ =& \frac{1}{2}nm(m+1)-g(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1) \end{align} \]

因此：

\[h(a,b,c,n)=nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n) \]

综上所述：

\[h(a,b,c,n)= \begin{cases} (\frac{b}{c})^2 &n=0\\ (n+1)\times(\frac{b}{c})^2 &a=0\\ h(a\%c,b\%c,c,n)+ 2 \left\lfloor \frac{b}{c} \right\rfloor f(a\%c,b\%c,c,n)+2 \left\lfloor \frac{a}{c} \right\rfloor g(a\%c,b\%c,c,n)+\\ \left\lfloor \frac{a}{c} \right\rfloor^2 \frac{n(n+1)(2n+1)}{6}+\left\lfloor \frac{b}{c} \right\rfloor^2 (n+1)+\left\lfloor \frac{a}{c} \right\rfloor \left\lfloor \frac{b}{c} \right\rfloor n(n+1) & a\geq c\ or\ b\geq c \\ nm(m+1)-2g(c,c-b-1,a,m-1)-2f(c,c-b-1,a,m-1)-f(a,b,c,n) & a<c\ and\ b<c \\ \end{cases} \]

在代码实现时，\(3\) 个函数会有交错递归，可以考虑 \(3\) 个一起整体递归，同步计算，否则会有很多项被多次计算，复杂度 \(O(\log n)\)。

模板题

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

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int mod=998244353;
int inv2=499122177,inv6=166374059;//mod下的2,6的逆元
struct data
{
    data()
    {
        f=g=h=0;
    }
    ll f,g,h;
};
data calc(ll a,ll b,ll c,ll n)
{
    ll ac=a/c,bc=b/c,m=(a*n+b)/c,nx=n+1,ny=2*n+1;
    data d;
    if(a==0)
    {
        d.f=bc*nx%mod;
        d.g=bc*n%mod*nx%mod*inv2%mod;
        d.h=bc*bc%mod*nx%mod;
        return d;
    }
    if(a>=c||b>=c)
    {
        d.f=n*nx%mod*inv2%mod*ac%mod+bc*nx%mod;
        d.g=ac*n%mod*nx%mod*ny%mod*inv6%mod+bc*n%mod*nx%mod*inv2%mod;
        d.h=ac*ac%mod*n%mod*nx%mod*ny%mod*inv6%mod+bc*bc%mod*nx%mod+ac*bc%mod*n%mod*nx%mod;
        d.f%=mod,d.g%=mod,d.h%=mod;
        data e=calc(a%c,b%c,c,n);// 迭代
        d.h+=e.h+2*bc%mod*e.f%mod+2*ac%mod*e.g%mod;
        d.g+=e.g,d.f+=e.f;
        d.f%=mod,d.g%=mod,d.h%=mod;
        return d;
    }
    data e=calc(c,c-b-1,a,m-1);
    d.f=n*m%mod-e.f,d.f=(d.f%mod+mod)%mod;
    d.g=m*n%mod*nx%mod-e.h-e.f,d.g=(d.g*inv2%mod+mod)%mod;
    d.h=n*m%mod*(m+1)%mod-2*e.g-2*e.f-d.f;
    d.h=(d.h%mod+mod)%mod;
    return d;
}
int main()
{
    int t;
    ll a,b,c,n;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
        data ans=calc(a,b,c,n);
        printf("%lld %lld %lld\n",ans.f,ans.h,ans.g);
    }
    return 0;
}