从这篇博客学的：link 。
解决这样的一类问题：
有一条直线 \(y=\frac{Px+B}{Q}\) ，其中 \(x\in(0,L],\mid B \mid<Q\) ，当直线穿过一条形如 \(y=h(h\in \mathbf{Z})\) 的直线的时候进行 \(U\) 操作，穿过一条形如 \(x=k(k\in \mathbf{Z})\) 的直线进行 \(R\) 操作，如果经过了一个整点就进行 \(UR\) 操作，问最后进行了所有操作之后的答案。
不妨把上面这个问题写作 \(f(P,Q,B,L,U,R)\) ，各个值和上面的定义相同。那么如果直接做的话，可以发现在进行第 \(a\) 个 \(R\) 操作前总共会进行 \(\lfloor \frac{Pa+B}{Q} \rfloor\) 个 \(U\) 操作，那么就是说这个操作序列可以用 \(\dots U^{\lfloor \frac{Pi+B}{Q}\rfloor-\lfloor\frac{P(i-1)+B}{Q}\rfloor}R\dots(i\in [1,L]\cap\mathbf{Z})\) 来表示，这里 \(U^a\) 表示 \(a\) 个 \(U\) 的拼接；同理一个长成这样的操作序列也可以用 \(f(P,Q,B,L,U,R)\) 来表示，即两者是等价的。
然后对 \(P,Q\) 的大小分类讨论来计算 \(f\) 。
如果 \(P\geq Q\) ，那么就意味着每次进行 \(R\) 操作前一定至少会有 \(\lfloor\frac{P}{Q}\rfloor\) 次 \(U\) 操作。
令 \(P=kQ+P'\) ，从上面的序列的角度考虑就是 \(U^{\lfloor \frac{Pi+B}{Q}\rfloor-\lfloor\frac{P(i-1)+B}{Q}\rfloor}R=U^{\lfloor \frac{P'i+B}{Q}\rfloor+ki-\lfloor\frac{P'(i-1)+B}{Q}\rfloor-k(i-1)}R=U^{\lfloor \frac{P'i+B}{Q}\rfloor-\lfloor\frac{P'(i-1)+B}{Q}\rfloor}U^kR\) 。
那么就可以得到 \(f(P,Q,B,L,U,R)=f(P',Q,B,L,U,U^kR)\) ，这样就变成了 \(P<Q\) 的问题。
如果 \(P<Q\) 那么会相对复杂一些。
注意到此时 \(U\) 前面会跟一些 \(R\) ，这启发我们去求第 \(b\) 个 \(U\) 前会有几个 \(R\) ，具体就是：\(b>\lfloor \frac{Pa+B}{Q} \rfloor \Rightarrow b>\frac{Pa+B}{Q}\Rightarrow a<\frac{Qb-B}{P}\Rightarrow a\leq \lfloor \frac{Qb-B-1}{P}\rfloor\) 。
可以注意到这个形式和上面是完全一致的，所以可以往子问题做。令 \(c=\lfloor \frac{PL+B}{Q} \rfloor\) 代表 \(U\) 的总出现次数，那么把前面 \(R^xU\) 的部分用 \(f(Q,P,-B-1,c,R,U)\) 表示，最后面的 \(L-\lfloor\frac{Qc-B-1}{P} \rfloor\) 个 \(R\) 单独计算，那么就是说 \(f(P,Q,B,L,U,R)=f(Q,P,-B-1,c,R,U)R^{L-\lfloor\frac{Qc-B-1}{P} \rfloor}\) 。
但是在具体实现的过程中会发现这样时错的，问题出在第一个 \(R^xU\) 上，因为 \(x\) 可能等于 \(0\) ，又由于定义中 \(x\in (0,L]\) ，所以会算错一部分内容，故把第一个 \(R^xU\) 单独拎出来计算，即 \(f(P,Q,B,L,U,R)=R^{\lfloor\frac{Q-B-1}{P} \rfloor}Uf(Q,P,Q-B-1,c-1,R,U)R^{L-\lfloor\frac{Qc-B-1} {P} \rfloor}\) 。在具体实现的时候还要让 \(Q-B-1\) 对 \(P\) 取模，因为在上面规定了 \(\mid B\mid\) ，并且因为直线的权值只和穿过横竖线有关，所以这样是对的。
计算复杂度可以观察到上面的过程和辗转相除法类似，所以复杂度为 \(\log V\) 。
最后就关于 \(U,R\) 的表示，一般可以写成矩阵的形式，这样比较方便。当然，写成矩阵就意味着更高的时间复杂度，所以尽量写成可合并多元组的形式，这样会更快。
例题的具体维护不想详细阐述，代码如下：
#include<bits/stdc++.h>
#define fo(i,a,b) for (int i=a;i<b;++i)
using namespace std;
typedef __int128_t i128;
typedef long long ll;
const int mod=998244353;
inline void add(int &x,int y){if ((x+=y)>=mod) x-=mod;}
inline int Mod(int x){return x>=mod?x-mod:x;}
struct node{
int cntu,cntr,sum,sqrsum,isum,mulsum;
node(int _cntu=0,int _cntr=0,int _sum=0,int _sqrsum=0,int _isum=0,int _mulsum=0){
cntu=_cntu,cntr=_cntr,sum=_sum,sqrsum=_sqrsum,isum=_isum,mulsum=_mulsum;
}
};
node U(1,0,0,0,0,0),R(0,1,0,0,1,0);
inline node operator + (node a,node b){
return node(
(a.cntu+b.cntu)%mod,(a.cntr+b.cntr)%mod,
(1ll*a.cntu*b.cntr%mod+a.sum+b.sum)%mod,
(1ll*a.cntu*a.cntu%mod*b.cntr%mod+2ll*a.cntu*b.sum%mod+a.sqrsum+b.sqrsum)%mod,
(1ll*a.cntr*b.cntr%mod+a.isum+b.isum)%mod,
(1ll*a.cntu*a.cntr%mod*b.cntr%mod+1ll*a.cntu*b.isum%mod+1ll*a.cntr*b.sum%mod+a.mulsum+b.mulsum)%mod
);
}
inline node ksm(node x,ll y){node c;for (;y;y>>=1,x=x+x)if (y&1) c=c+x;return c;}
inline node f(ll P,ll Q,ll B,ll L,node U,node R){
if (!L) return node();
if (P>=Q) return f(P%Q,Q,B,L,U,ksm(U,P/Q)+R);
ll c=static_cast<i128>((P*L+B)/Q);
if (!c) return ksm(R,L);
return ksm(R,(Q-B-1)/P)+U+f(Q,P,(Q-B-1)%P,c-1,R,U)+ksm(R,L-(Q*c-B-1)/P);
}
void work(){
ll n,a,b,c;
scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
node ans=node((b/c)%mod,0,(b/c)%mod,1ll*(b/c)*(b/c)%mod,0,0);
ans=ans+f(a,c,b%c,n,U,R);
printf("%d %d %d\n",ans.sum,ans.sqrsum,ans.mulsum);
}
int main(){
int T;scanf("%d",&T);
while (T--) work();
return 0;
}