类欧几里得算法推导
毒瘤中的毒瘤。
类欧几里得算法
给定 \(n,a,b,c\)，分别求
答案对 \(998244353\) 取模。多组数据。
\(1\le T\le10^5,\ 0\le n,a,b,c,\ c\neq0\)。
推导
分别令
推导 \(f\)
我们可以考虑分两类进行：一部分 \(a\ge c\vee b\ge c\)；另一部分 \(a<c\wedge b < c\)。对于前一部分，如果可以通过变换缩小范围，就可以化为第二部分。下面我们来分别推导两个部分的式子。
PART 1
化简为繁，我们将 \(a,b\) 拆成 \(\bmod c\) 的形式，尝试减小范围
若只留下 \(a\bmod c\) 和 \(b\bmod c\) 项，提出其余项
后面的和式可以再用 \(f\) 替换
现在新得到的 \(f\) 中的 \(a,b\) 都小于 \(c\) 了。接下来继续考虑这种情况的推导。
PART 2
化简为繁，我们把 \(\lfloor{ai+b\over c}\rfloor\) 用和式展开
交换和式进一步观察
尝试转化右边的条件：
根据下取整的性质
把 \(i\) 单独挪到一边
令 \(m=\lfloor{an+b\over c}\rfloor\)，再把上面的条件带回原式
对于每个 \(i\in[0,n]\)，要么大于 \(\lfloor{jc+c-b-1\over a}\rfloor\) 要么小于，总共的情况数就是 \(n-\lfloor{jc+c-b-1\over a}\rfloor\)。所以上式即为：
用 \(f\) 表示上式，得到递归式
递归调用上式即可求解 \(f\)。上式的 \(a\) 和 \(c\) 发生了交换，时间复杂度同欧几里得算法 \(O(\log n)\)。
推导 \(g,h\)
把 \(g,h\) 放一起的原因是因为它们的递归式之间存在一个依赖关系，相互包含。与 \(f\) 类似的，我们可以分成两部分来考虑。
PART 1
先考虑对于 \(a\ge c\vee b\ge c\) 的情况的 \(g()\)，类似的，有
对于 \(h()\) 在 \(a\ge c\vee b\ge c\)，先拆成 \(\bmod c\) 有
发现对于最外面的下取整，括号里面都是整数，于是可以改写成小括号，并且再展开平方，有
\(
\begin{aligned}
&\sum_{i=0}^n(\lfloor{a\over c}\rfloor^2i^2+\lfloor{b\over c}\rfloor^2+\lfloor{a\bmod c\times i+b\bmod c\over c}\rfloor^2+2\lfloor{a\over c}\rfloor\lfloor{b\over c}\rfloor i+2\lfloor{a\over c}\rfloor\lfloor{a\bmod c\times i+b\bmod c\over c}\rfloor+2\lfloor{b\over c}\rfloor\lfloor{a\bmod c\times i+b\bmod c\over c}\rfloor)\\
=&{n(n+1)(2n+1)\over 6}\lfloor{a\over c}\rfloor^2+(n+1)\lfloor{b\over c}\rfloor^2+h(a\bmod c,b\bmod c,c,n)+n(n+1)\lfloor{a\over c}\rfloor\lfloor{b\over c}\rfloor+2\lfloor{a\over c}\rfloor g(a\bmod c,b\bmod c,c,n)+2\lfloor{b\over c}\rfloor f(a\bmod c,b\bmod c,c,n)
\end{aligned}
\)
成功把 \(g()\) 和 \(h()\) 都化成小于 \(c\) 的形式了，下面我们来进一步推导小于 \(c\) 的情形。
PART 2
对于 \(g()\)，令 \(m=\lfloor{an+b\over c}\rfloor\)，与 \(f()\) 类似的，有
再令 \(t=\lfloor{jc+c-b-1\over a}\rfloor\)，有
最后用 \(h()\) 和 \(f()\) 替换掉对 \(t\) 求和的项，得到递推式
递归调用就能算了。
对于 \(h()\)，把 \(x^2\) 写成 \((2\sum_{i=1}^xi)-x\)，有
令 \(m=\lfloor{an+b\over c}\rfloor, t=\lfloor{cj+c-b-1\over a}\rfloor\)，就有
其中左式又可以用前文提到的方法处理一下，即有
所以 \(h()\) 的递归式就是
可以递归求解
方法总结
以上的推导蕴含了几种处理递归式和关于顶和底和式的思想方法：
- 对值的分布进行分类讨论，统一化成一种情形
- 把难以处理的顶和底的式子用艾弗森括号来表达，再考虑化简这个条件
- 惊人的计算准确性是必要的
代码实现
递归式
根据以上的推导，我们可以得到 \(f(),g(),h()\) 各两个递归式，即
\(
f(a,b,c,n)=\begin{cases}
{n(n+1)\over2}\lfloor{a\over c}\rfloor+(n+1){\lfloor{b\over c}\rfloor}+f(a\bmod c,b\bmod c,c,n) &(a\ge c\vee b\ge c)\\
mn-f(c,c-b-1,a,m-1),\ m=\lfloor{an+b\over c}\rfloor &(a< c\wedge b< c)
\end{cases}
\)
\( g(a,b,c,n)=\begin{cases} {n(n+1)(2n+1)\over 6}\lfloor{a\over c}\rfloor+{n(n+1)\over2}\lfloor{b\over c}\rfloor+g(a\bmod c,b\bmod c, c, n) &(a\ge c\vee b\ge c)\\ {1\over2}[mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)],\ m=\lfloor{an+b\over c}\rfloor &(a< c\wedge b< c) \end{cases} \)
\(
h(a,b,c,n)=\begin{cases}
{n(n+1)(2n+1)\over 6}\lfloor{a\over c}\rfloor^2+(n+1)\lfloor{b\over c}\rfloor^2+h(a\bmod c,b\bmod c,c,n)+n(n+1)\lfloor{a\over c}\rfloor\lfloor{b\over c}\rfloor+2\lfloor{a\over c}\rfloor g(a\bmod c,b\bmod c,c,n)+2\lfloor{b\over c}\rfloor f(a\bmod c,b\bmod c,c,n)
&(a\ge c\vee b\ge c)\\
m(m+1)n-2f(c,c-b-1,a,m-1)-2g(c,c-b-1,a,m-1)-f(a,b,c,n),\ m=\lfloor{an+b\over c}\rfloor &(a< c\wedge b< c)
\end{cases}
\)
递归调用，可以用一个结构体 node{ll f, g, h;}; 把三个函数同步求出来，免得写记忆化不好调还跑得慢。
复杂度与边界条件
我们来考虑边界条件是什么。如果 \(a>c\vee b>c\)，对于初始的 \((a,b,c,n)\) 递归一次得到 \((a\bmod c,b\bmod c,c,n)\)，递归第二次再得到 \((c, c-b\bmod c-1, a\bmod c,m-1)\)，发现对于 \((a,c)\) 经过两次递归得到 \((c,a\bmod c)\)，那说明类欧几里得算法的复杂度就是欧几里得算法的复杂度 \(O(\log(a+c))\)。自然而然的，边界条件就应该是 \(a=0\)。
代码
注意这里 f, g, h 的定义与原题不太一样，与博客上文保持一致。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mo = 998244353, inv2 = 499122177, inv6 = 166374059;
int T;
struct node {
ll f, g, h;
void out() {cout << f << " " << h << " " << g << endl;}
};
node solve(ll n, ll a, ll b, ll c) {
if(!a) return (node){(n + 1) * (b / c) % mo, n * (n + 1) % mo * inv2 % mo * (b / c) % mo, (n + 1) * (b / c) % mo * (b / c) % mo};
if(a >= c || b >= c) {
ll t1 = a / c, t2 = b / c, s1 = n * (n + 1) % mo * inv2 % mo, s2 = n * (n + 1) % mo * (2 * n + 1) % mo * inv6 % mo;
node tmp = solve(n, a % c, b % c, c), res;
res.f = (s1 * t1 % mo + (n + 1) * t2 % mo + tmp.f) % mo; (res.f += mo) %= mo;
res.g = (s2 * t1 % mo + s1 * t2 % mo + tmp.g) % mo; (res.g += mo) %= mo;
res.h = (s2 * t1 % mo * t1 % mo + (n + 1) * t2 % mo * t2 % mo + tmp.h + 2 * s1 % mo * t1 % mo * t2 % mo + 2 * t1 * tmp.g % mo + 2 * t2 * tmp.f % mo) % mo, (res.h += mo) %= mo;
return res;
}
ll m = (a * n + b) / c; node tmp = solve(m - 1, c, c - b - 1, a), res;
res.f = (m * n % mo - tmp.f) % mo, (res.f += mo) %= mo;
res.g = (m * n % mo * (n + 1) % mo - tmp.h - tmp.f) % mo * inv2 % mo, (res.g += mo) %= mo;
res.h = (m * (m + 1) % mo * n % mo - 2 * tmp.f - 2 * tmp.g - res.f) % mo, (res.h += mo) %= mo;
return res;
}
int main() {
ios :: sync_with_stdio(false); cin.tie(0); cout.tie(0);
cin >> T;
while(T--) {
ll n, a, b, c;
cin >> n >> a >> b >> c;
solve(n, a, b, c).out();
}
return 0;
}
取模挺多的，注意不要漏了。
浙公网安备 33010602011771号