类欧几里得算法

难得有空写学习笔记。

写这篇笔记的时候同学在机房互吐冰红茶。

问题：给定 $n,a,b,c$，请你分别求出 $f(a,b,c,n)=\sum _{i=0}^n\lfloor \frac{ai+b}{c}\rfloor$，$g(a,b,c,n)=\sum _{i=0}^{n}i\lfloor \frac{ai+b}{c}\rfloor$，$h(a,b,c,n)=\sum _{i=0}^n{\lfloor \frac{ai+b}{c}\rfloor }^2$。

---

Part 1

考虑如何求 $f(a,b,c,n)$，分成两部分求解：

当 $a,b>c$ 时：

则 $f(a,b,c,n)=\sum _{i=0}^n\lfloor \frac{amodc\times i+bmodc+\lfloor \frac{a}{c}\rfloor \times c\times i+\lfloor \frac{b}{c}\rfloor \times c}{c}\rfloor$。

把常数项提出来即为 $f(a,b,c,n)=\lfloor \frac{a}{c}\rfloor \times \frac{n(n+1)}{2}+\lfloor \frac{b}{c}\rfloor \times (n+1)+\sum _{i=0}^n\lfloor \frac{amodc\times i+bmodc}{c}\rfloor$。

换一下求和部分：$f(a,b,c,n)=\lfloor \frac{a}{c}\rfloor \times \frac{n(n+1)}{2}+\lfloor \frac{b}{c}\rfloor \times (n+1)+f(amodc,bmodc,c,n)$，发现此时 $amodc,bmodc$ 都小于 $c$。

当 $a,b<c$ 时：

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

把 $j$ 提到外面去变为 $f(a,b,c,n)=\sum _{j=0}^{\lfloor \frac{an+b}{c}\rfloor -1}\sum _{i=0}^n[j<\lfloor \frac{ai+b}{c}\rfloor ]$。

考虑转换条件式：

$j<\lfloor \frac{ai+b}{c}\rfloor \Rightarrow j\le \frac{ai+b}{c}-1\Rightarrow i\ge \frac{jc+c-b}{a}\Rightarrow i>\lfloor \frac{jc+c-b-1}{a}\rfloor$。

令 $m=\lfloor \frac{an+b}{c}\rfloor$，则 $f(a,b,c,n)=\sum _{j=0}^{m-1}\sum _{i=0}^n[i>\lfloor \frac{jc+c-b-1}{a}\rfloor ]$。

最后 $f(a,b,c,n)=\sum _{j=0}^{m-1}(n-\lfloor \frac{jc+c-b-1}{a}\rfloor )$，即为 $f(a,b,c,n)=nm-f(c,c-b-1,a,m-1)$。

---

Part 2

考虑求 $g(a,b,c,n)$。

同样的思路可以算出当 $a,b>c$ 时： $g(a,b,c,n)=\frac{n(n+1)(2n+1)}{6}\times \lfloor \frac{a}{c}\rfloor +\frac{n(n+1)}{2}\times \lfloor \frac{b}{c}\rfloor +g(amodc,bmodc,c,n)$。

当 $a,b<c$ 时，令 $m=\lfloor \frac{an+b}{c}\rfloor$，$t=\lfloor \frac{jc+c-b-1}{a}\rfloor$。

则 $g(a,b,c,n)=\sum _{j=0}^{m-1}\sum _{i=0}^{n}i[i>t]$。

不难得到 $g(a,b,c,n)=\sum _{j=0}^{m-1}\frac{1}{2}(n(n+1)-t(t+1))$。

提出常数项得到 $g(a,b,c,n)=\frac{1}{2}(nm(n+1)-\sum _{j=0}^{m-1}t-\sum _{j=0}^{m-1}{t}^2)$。

最后把 $t$ 带进去得到 $g(a,b,c,n)=\frac{1}{2}(nm(n+1)-f(c,c-b-1,a,m-1)-h(c,c-b-1,a,m-1))$。

---

Part 3

考虑求 $h(a,b,c,n)$。

当 $a,b>c$ 时：

考虑取模变形：$h(a,b,c,n)=\sum _{i=0}^n(i\lfloor \frac{a}{c}\rfloor +\lfloor \frac{b}{c}\rfloor +\frac{amodc\times i+bmodc}{c}{)}^2$。

展开可得：

$h(a,b,c,n)=\sum _{i=0}^n(i^2{\lfloor \frac{a}{c}\rfloor }^2+{\lfloor \frac{b}{c}\rfloor }^2+{\lfloor \frac{amodc\times i+bmodc}{c}\rfloor }^2+2i\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor +2i\lfloor \frac{a}{c}\rfloor \lfloor \frac{amodc\times i+bmodc}{c}\rfloor +2\lfloor \frac{b}{c}\rfloor \lfloor \frac{amodc\times i+bmodc}{c}\rfloor )$。

经过整理可得：

$h(a,b,c,n)=h(amodc,bmodc,c,n)+S_2(n){\lfloor \frac{a}{c}\rfloor }^2+(n+1){\lfloor \frac{b}{c}\rfloor }^2+2S_1(n)\lfloor \frac{a}{c}\rfloor \lfloor \frac{b}{c}\rfloor +2\lfloor \frac{a}{c}\rfloor g(amodc,bmodc,c,n)+2\lfloor \frac{b}{c}\rfloor f(amodc,bmodc,c,n)$。

其中 $S_1(n)=\frac{n(n+1)}{2},S_2(n)=\frac{n(n+1)(2n+1)}{6}$。

当 $a,b<c$ 时：

一个简单有用的 trick：$x^2=2\frac{x(x+1)}{2}-x=2\sum _{i=1}^{x}i-x$。

然后带入原式：$h(a,b,c,n)=\sum _{i=0}^{n}2(\sum _{j=0}^{\lfloor \frac{ai+b}{c}\rfloor }j)-\lfloor \frac{ai+b}{c}\rfloor$。

同样地，令 $m=\lfloor \frac{an+b}{c}\rfloor$，$t=\lfloor \frac{jc+c-b-1}{a}\rfloor$。

做最后的变形：

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

$=2\sum _{j=0}^{m-1}(j+1)(n-t)-f(a,b,c,n)$

$=2\sum _{j=0}^{m-1}(j+1)(n-\lfloor \frac{jc+c-b-1}{a}\rfloor )-f(a,b,c,n)$

$=nm(m+1)-2f(c,c-b-1,a,m-1)-2g(c,c-b-1,a,m-1)-f(a,b,c,n)$

---

最后一起考虑 $a=0$ 时的情况，也就是把含 $a$ 的项去掉即可。

三个公式一起递归求解即可，时间复杂度为 $O(\log n)$。

代码回家再写。

全机房在摸鱼。我也要颓颓颓

---

Code

#include<bits/stdc++.h> 
#define int long long 
#define x first 
#define y second 
#define il inline 
#define debug() puts("-----") 
using namespace std; 
typedef pair<int,int> pii; 
il int read(){ 
	int x=0,f=1; char ch=getchar(); 
	while(ch<'0'||ch>'9'){ if(ch=='-') f=-1; ch=getchar(); } 
	while(ch>='0'&&ch<='9') x=(x<<1)+(x<<3)+(ch^48),ch=getchar(); 
	return x*f; 
} 
const int Mod=998244353; 
const int inv2=499122177;  
const int inv6=166374059; 
struct Node{ int f,g,h; };  
il int s1(int x){ return x*(x+1)%Mod*inv2%Mod; } 
il int s2(int x){ return x*(x+1)%Mod*(2*x+1)%Mod*inv6%Mod; }  
il Node solve(int a,int b,int c,int n){ 
	Node ans,res; 
	int Ac=a/c,Bc=b/c; 
	int Ac2=Ac*Ac%Mod,Bc2=Bc*Bc%Mod; 
	if(a==0){ 
		ans.f=Bc*(n+1)%Mod; 
		ans.g=Bc*s1(n)%Mod; 
		ans.h=Bc2*(n+1)%Mod; 
		return ans; 
	} if(a>=c||b>=c){ 
		res=solve(a%c,b%c,c,n); 
		ans.f=(res.f+Ac*s1(n)%Mod+Bc*(n+1)%Mod)%Mod; 
		ans.g=(res.g+Ac*s2(n)%Mod+Bc*s1(n)%Mod)%Mod; 
		ans.h=(res.h+Ac2*s2(n)%Mod+Bc2*(n+1)%Mod)%Mod; 
		ans.h=(ans.h+2*(s1(n)*Ac%Mod*Bc%Mod+Ac*res.g%Mod+Bc*res.f%Mod)%Mod)%Mod; 
		return ans;  
	} int m=(a*n+b)/c; 
	res=solve(c,c-b-1,a,m-1); 
	ans.f=(n*m%Mod-res.f)%Mod; 
	ans.g=(n*m%Mod*(n+1)%Mod-res.f-res.h)%Mod*inv2%Mod; 
	ans.h=(n*m%Mod*(m+1)%Mod-2*(res.f+res.g)%Mod-ans.f)%Mod; 
	ans.f=(ans.f+Mod)%Mod,ans.g=(ans.g+Mod)%Mod,ans.h=(ans.h+Mod)%Mod;  
	return ans; 
} 
il void work(){ 
	int n=read(),a=read(),b=read(),c=read(); 
	Node ans=solve(a,b,c,n); printf("%lld %lld %lld\n",ans.f,ans.h,ans.g); 
} 
signed main(){ 
	int T=read(); 
	while(T--) work(); 
	return 0; 
}