万能欧几里得

有一条定义域是 \((0,L]\) 的直线 \(l:y=\frac{ax+b}{c}\)，有两种操作是 \(P,Q\)，一旦 \(l\) 碰到整横线那么执行 \(P\)，碰到整竖线则执行 \(Q\)，碰到整点默认先执行 \(P\) 再 \(Q\)，问直线跑完的结果。

执行 \(k\) 次操作 \(A\) 写作 \(A^k\)。

用一个 \(f(a,b,c,L,P,Q)\) 表示执行这样一个过程，我们想把它化归成更小的问题来做。

先发现将 \(l\) 上下平移整数格是不影响的，于是我们可以把 \(b\) 拉到 \([0,c)\) 这个范围，即 \(f(a,b,c,L,P,Q)=f(a,b\bmod c,c,L,P,Q)\)。

当 \(a\geq c\) 时，发现在碰到一个竖线前必碰到 \(\lfloor\frac{a}{c}\rfloor\) 条横线，将这些个 \(P\) 缩成一个操作，于是有 \(f(a,b,c,L,P,Q)=f(a\bmod c,b,c,L,P,P^{\lfloor\frac{a}{c}\rfloor}Q)\)。

当 \(a<c\) 时，想化成前面的样子，想到交换一下坐标轴。现在考虑交换 \(x,y\) 轴，原坐标系里面定义域是 \((0,L]\)，值域是 \((\frac{b}{c},\frac{aL+b}{c}]\)，这两个交换，以及 \(P,Q\) 两种操作交换，显然有 \(\frac{b}{c}\in [0,1)\)。

首先若 \(\frac{aL+b}{c}<1\) 那么 \(l\) 不会经过竖线了，而经过 \(L\) 条横线，所以直接执行 \(Q^L\) 并结束。

否则，考虑到原问题的定义域是 \((0,k](k\in N^{+})\) 这样子，为了符合这个形式，再次注意到 \(\frac{b}{c}\in [0,1)\)，我们将 \((\frac{b}{c},\frac{aL+b}{c}]\) 从 \(1\) 处截断，再将 \((1,\frac{aL+b}{c}]\) 这部分平移到 \((0,\frac{aL+b}{c}-1]\)，最后把后面的分数值截取到整数。还有个小问题即是如果经过整点，那么 \(PQ\) 的顺序会反转，这里做个微扰，将 \(l\) 向下平移 \(\frac{1}{a}\) 即可。于是最终 \(l':y=\frac{cx+c-b-1}{a}\)后半部分化为以 \(f(c,c-b-1,a,\lfloor\frac{aL+b}{c}\rfloor-1,Q,P)\) 解决。

最后我们要处理 \((\frac{b}{c},1]\)，\((\lfloor\frac{aL+b}{c}\rfloor,\frac{aL+b}{c}]\) 这两个边界定义域。对于首段，交且仅交一个竖线 \(x=1\)，交 \(\lfloor\frac{c-b}{a}\rfloor\) 个横线，但是注意到中间那部分下移 \(\frac{1}{a}\) 的操作，假如 \(\frac{c-b}{a}\) 是个整数那么会被那边的多统计一次，所以我们不妨也直接把这个和中间对齐，就是 \(\lfloor\frac{c-b-1}{a}\rfloor\) 条横线。

再看尾段，这里显然不会再经过竖线了，令 \(t=\lfloor\frac{aL+b}{c}\rfloor\) 和上面同理，横线会经过 \(L-\lfloor\frac{ct-b-1}{a}\rfloor\) 条。

发现这个有点像三个数的欧几里得，如果做操作的代价是 \(O(t)\) 的，那么总复杂度 \(O(t\log V)\)。

如果题里面的定义域包括 \(0\)，那么要在 \(0\) 处向上跑截距条横线再向右跑一格再进行。

把一个减的直线换成增的扔到原点上且整点数不变：\(l:y=\frac{-ax+b}{c}\rightarrow l':y=\frac{ax-b\bmod a}{c}\)。

lg 模板的代码：

# include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef long double ldb;
const ll Mod = 998244353;


struct Element{
	ll s , u , r , ss , zs , si;
	Element operator + ( Element y ) { 
		Element ret = { 0 , 0 , 0 , 0 , 0 , 0 };
		ret . u = ( u + y . u ) % Mod , ret . r = ( r + y . r ) % Mod;
		ret . s = ( s + y . s + ( u * y . r % Mod ) ) % Mod;
		ret . ss = ( ss + y . ss + ( 2 * u * y . s % Mod ) + ( u * u % Mod * y . r % Mod ) ) % Mod;
		ret . si = ( si + y . si + r * y . r % Mod ) % Mod;
		ret . zs = ( zs + y . zs + ( u * y . si % Mod ) + ( u * r % Mod * y . r % Mod ) + ( r * y . s ) % Mod ) % Mod;
		return ret;
	}
};

Element Quick_Pow( Element u , ll k ){
	Element ret = { 0 , 0 , 0 , 0 , 0 , 0 };
	while( k ){
		if( k & 1 ) ret = ret + u;
		u = u + u , k >>= 1;
	}
	return ret;
}

Element Solve( ll a , ll b , ll c , ll L , Element P , Element Q ){
	if( ! L ) return ( Element ){ 0 , 0 , 0 , 0 , 0 , 0 };
	b %= c;
	if( a >= c ) return Solve( a % c , b , c , L , P , Quick_Pow( P , a / c ) + Q );
	ll t = ( a * L + b ) / c;
	if( ! t ) return Quick_Pow( Q , L );
	ll k = L - ( c * t - b - 1 ) / a , d = ( c - b - 1 ) / a;
	return Quick_Pow( Q , d ) + P + Solve( c , c - b - 1 , a , t - 1 , Q , P ) + Quick_Pow( Q , k );
}

int main(){
	ll T;
	scanf( "%lld" , & T );
	while( T -- ){
		ll L , a , b , c;
		scanf( "%lld%lld%lld%lld" , & L , & a , & b , & c );
		Element _P = { 0 , 1 , 0 , 0 , 0 , 0 } , _Q = { 0 , 0 , 1 , 0 , 0 , 0 };
		Element ans = Quick_Pow( _P , b / c ) + _Q + Solve( a , b , c , L , _P , _Q );
		printf( "%lld %lld %lld\n" , ans . s , ans . ss , ans . zs );
	}	
	return 0;
}