类欧几里得算法

基础

设

[f(a,b,c,n)=sum_{i=0}^{n}{leftlfloor frac{ai+b}{c} ight floor} ]

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

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

[egin{align} f(a,b,c,n) &= sum_{i=0}^{n}{leftlfloor frac{ai+b}{c} ight floor}\ &=sum_{i=0}^{n}{leftlfloor frac{(lfloor frac{a}{c} floor c+a mod c)i+(lfloor frac{b}{c} floor c+b mod c)}{c} ight floor}\ &=frac{(1+n)n}{2}·leftlfloor frac{a}{c} ight floor+(n+1)·leftlfloor frac{b}{c} ight floor+sum_{i=0}^{n}{ leftlfloor frac{(a mod c)i+(b mod c)}{c} ight floor}\ &=frac{(1+n)n}{2}·leftlfloor frac{a}{c} ight floor+(n+1)·leftlfloor frac{b}{c} ight floor+f(a mod c ,b mod c,c,n)\ end{align} ]

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

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

[f(a,b,c,n) = sum_{i=0}^{n}{leftlfloor frac{ai+b}{c} ight floor}=sum_{i=0}^{n}{sum_{j=0}^{leftlfloor frac{ai+b}{c} ight floor-1}{1}} ]

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

[Rightarrow sum_{j=0}^{leftlfloor frac{an+b}{c} ight floor-1}{sum_{i=0}^{n}{left[j<leftlfloor frac{ai+b}{c} ight floor ight]}} ]

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

先将取整符号拿掉：

[j< leftlfloor frac{ai+b}{c} ight floor Leftrightarrow j+1 leq leftlfloor frac{ai+b}{c} ight floor 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 leftlfloor frac{jc+c-b-1}{a} ight floor <i ]

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

[egin{align} f(a,b,c,n) &= sum_{j=0}^{m-1}{sum_{i=0}^{n}{left[ i>leftlfloor frac{jc+c-b-1}{a} ight floor ight]}}\ &= sum_{j=0}^{m-1}{left(n-leftlfloor frac{jc+c-b-1}{a} ight floor ight)}\ &= nm-sum_{j=0}^{m-1}{leftlfloor frac{jc+c-b-1}{a} ight floor}\ &= nm-f(c,c-b-1,a,m-1)\ &= nleftlfloor frac{an+b}{c} ight floor-f(c,c-b-1,a,leftlfloor frac{an+b}{c} ight floor-1) end{align} ]

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

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

综上所述：

[f(a,b,c,n)= egin{cases} lfloor frac{b}{c} floor & n=0\ (n+1) imes lfloor frac{b}{c} floor & a=0\ frac{(1+n)·n}{2} imes leftlfloor frac{a}{c} ight floor+(n+1) imes leftlfloor frac{b}{c} ight floor+f(a \% c ,b\%c,c,n) & ageq c or bgeq c\ n imes leftlfloor frac{an+b}{c} ight floor-f(c,c-b-1,a,leftlfloor frac{an+b}{c} ight floor-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}{ileftlfloor frac{ai+b}{c} ight floor}\ h(a,b,c,n)=sum_{i=0}^{n}{leftlfloor frac{ai+b}{c} ight floor^2} ]

g 函数求解

首先，考虑取模：

[g(a,b,c,n)=g(a\%c,b\%c,c,n)+frac{n(n+1)(2n+1)}{6} leftlfloor frac{a}{c} ight floor+frac{n(n+1)}{2}leftlfloor frac{b}{c} ight floor ]

接下来，考虑 (a<c,b<c) 的情况，令 (m=leftlfloor frac{an+b}{c} ight floor)，有：

[g(a,b,c,n)=sum_{i=0}^{n}{ileftlfloor frac{ai+b}{c} ight floor}=sum_{j=0}^{m-1}{sum_{i=0}^{n}{i·left[j<leftlfloor frac{ai+b}{c} ight floor ight]}} ]

令 (t=leftlfloor frac{jc+c-b-1}{a} ight floor)，同理有：

[egin{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} ight)\ &= frac{1}{2} left(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1) ight) end{align} ]

综上所述：

[g(a,b,c,n)=egin{cases} 0 & n=0\ frac{n(n+1)}{2} imes leftlfloor frac{b}{c} ight floor & a=0\ g(a\%c,b\%c,c,n)+frac{n(n+1)(2n+1)}{6} leftlfloor frac{a}{c} ight floor+frac{n(n+1)}{2}leftlfloor frac{b}{c} ight floor & ageq c or bgeq c \ frac{1}{2} left(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1) ight) & a<c and b<c\ end{cases} ]

h 函数求解

首先取模：

[egin{align} h(a,b,c,n) &= h(a\%c,b\%c,c,n)\ &+ 2 leftlfloor frac{b}{c} ight floor f(a\%c,b\%c,c,n)+2 leftlfloor frac{a}{c} ight floor g(a\%c,b\%c,c,n)\ &+ leftlfloor frac{a}{c} ight floor^2 frac{n(n+1)(2n+1)}{6}+leftlfloor frac{b}{c} ight floor^2 (n+1)+leftlfloor frac{a}{c} ight floor leftlfloor frac{b}{c} ight floor n(n+1) end{align} ]

考虑 (a<c,b<c) 的情况，令 (m=leftlfloor frac{ai+b}{c} ight floor)，(t=leftlfloor frac{jc+c-b-1}{a} ight floor)。

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

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

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

[egin{align} h(a,b,c,n) &= sum_{i=0}^{n}{leftlfloor frac{ai+b}{c} ight floor^2}\ &= sum_{i=0}^{n}{left[left(2sum_{j=0}^{leftlfloor frac{ai+b}{c} ight floor}{j} ight)-leftlfloor frac{ai+b}{c} ight floor ight]}\ &= left(2sum_{i=0}^{n}{sum_{j=0}^{leftlfloor frac{ai+b}{c} ight floor}{j}} ight)-f(a,b,c,n)\ end{align} ]

接下来化简前一部分：

[egin{align} & sum_{i=0}^{n}{sum_{j=0}^{leftlfloor frac{ai+b}{c} ight floor}{j}}\ =& sum_{i=0}^{n}{sum_{j=0}^{leftlfloor frac{ai+b}{c} ight floor-1}{(j+1)}}\ =& sum_{j=0}^{m-1}{(j+1)sum_{i=0}^{n}{left[j<leftlfloorfrac{ai+b}{c} ight floor ight]}}\ =& sum_{j=0}^{m-1}{(j+1)sum_{i=0}^{n}{left[i>t ight]}}\ =& sum_{j=0}^{m-1}{(j+1)(n-t)}\ =& frac{1}{2}nm(m+1)-sum_{j=0}^{m-1}{(j+1)leftlfloor frac{jc+c-b-1}{a} ight floor}\ =& 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)= egin{cases} (frac{b}{c})^2 &n=0\ (n+1) imes(frac{b}{c})^2 &a=0\ h(a\%c,b\%c,c,n)+ 2 leftlfloor frac{b}{c} ight floor f(a\%c,b\%c,c,n)+2 leftlfloor frac{a}{c} ight floor g(a\%c,b\%c,c,n)+\ leftlfloor frac{a}{c} ight floor^2 frac{n(n+1)(2n+1)}{6}+leftlfloor frac{b}{c} ight floor^2 (n+1)+leftlfloor frac{a}{c} ight floor leftlfloor frac{b}{c} ight floor n(n+1) & ageq c or bgeq 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
",ans.f,ans.h,ans.g);
    }
    return 0;
}