类欧几里得算法
对于给定的元(a,b,c,n)
设(f(i)=lfloorfrac{ai+b}{c} floor)
求
[F(a,b,c,n)=sum_0^nf(i)
]
[G(a,b,c,n)=sum_0^nf(i)^2
]
[H(a,b,c,n)=sum_0^nicdot f(i)
]
Part1
(age c) 或 (bge c)
[lfloorfrac{ai+b}{c}
floor=lfloorfrac{(a mod c)i+(b mod c)}{c}
floor+ilfloor frac{a}{c}
floor +lfloor frac{b}{c}
floor
]
[F(a,b,c,n)=F(amod c,bmod c,c,n)+frac{n(n+1)lfloorfrac{a}{c}
floor}{2}+(n+1)lfloorfrac{b}{c}
floor
]
(a<c) 且 (b<c)
[F(a,b,c,n) =sum_0^n f(i)
]
[= sum_{i=0}^nsum_{j=0}^{f(i)-1} 1
]
[= sum_{i=0}^nsum_{j=0}^{infty}[j<f(i)]
]
[= sum_{j=0}^{f(n) -1}sum_{i=0}^{n}[j<f(i)]
]
其中
[[j<lfloorfrac{ai+b}{c}
floor] =[j<lceilfrac{ai+b-c+1}{c}
ceil]
]
[=[cj<ai+b-c+1]
]
[=[ai>cj-b+c-1]
]
[=[i>lfloorfrac{cj-b+c-1}{a}
floor]
]
[F(a,b,c,n)= sum_{j=0}^{f(n)-1 }sum_{i=0}^{n}[i>lfloorfrac{cj-b+c-1}{a}
floor]
]
[=sum_{j=0}^{f(n)-1} n-lfloorfrac{cj-b+c-1}{a}
floor
]
[F(a,b,c,n)=ncdot f(n)-F(c,-b+c-1,a,f(n)-1)
]
每次操作交换了(a,c)然后再次取模，复杂度是(O(log n))
递归边界是(a=0)
Part2
(G(a,b,c,n)=sum_0^nf(i)^2)
(age c) 或 (bge c)
设(lfloor frac{a}{c} floor=x,lfloor frac{b}{c} floor=y)
[lfloorfrac{ai+b}{c}
floor^2=(lfloorfrac{(a mod c)i+(b mod c)}{c}
floor+ix+y)^2
]
[=lfloorfrac{(a mod c)i+(b mod c)}{c}
floor^2+2lfloorfrac{(a mod c)i+(b mod c)}{c}
floorcdot (ix+y)+(ix+y)^2
]
[G(a,b,c,n)=
]
[G(amod c,bmod c,c,n)+2x H(amod c,bmod c,n)+2y F(amod c,bmod c,c,n)
]
[+frac{n(n+1)(2n+1)x^2}{6}+(n+1)y^2+n(n+1)xy
]
(a<c) 且 (b<c)
[G(a,b,c,n) =sum_0^nf(i)^2
]
[= 2 sum_{i=0}^nfrac{f(i)(f(i)-1)}{2}+frac{f(i)}{2}
]
[=2 sum_{i=0}^nsum_{j=0}^{f(i)-1}j+frac{1}{2}
]
[= 2sum_{i=0}^nsum_{j=0}^{infty}(j+frac{1}{2})[j<f(i)]
]
[= 2sum_{j=0}^{f(n)-1}sum_{i=0}^{n}(j+frac{1}{2})[j<f(i)]
]
[=2 sum_{j=0}^{f(n)-1}sum_{i=0}^{n}(j+frac{1}{2})[i>lfloorfrac{cj-b+c-1}{a}
floor]
]
[=2sum_{j=0}^{f(n)-1} (j+frac{1}{2})(n-lfloorfrac{cj-b+c-1}{a}
floor)
]
[G(a,b,c,n)=nf(n)^2-2H(c,-b+c-1,a,f(n)-1)-F(c,-b+c-1,a,f(n)-1)
]
Part3
(H(a,b,c,n)=sum_0^nicdot f(i))
(age c) 或 (bge c)
[icdot lfloorfrac{ai+b}{c}
floor=icdot (lfloorfrac{(a mod c)i+(b mod c)}{c}
floor+ilfloor frac{a}{c}
floor +lfloor frac{b}{c}
floor)
]
[icdot lfloorfrac{ai+b}{c}
floor=icdot lfloorfrac{(a mod c)i+(b mod c)}{c}
floor+i^2lfloor frac{a}{c}
floor +ilfloor frac{b}{c}
floor
]
[H(a,b,c,n)=H(amod c,bmod c,c,n)+frac{n(n+1)(2n+1)lfloorfrac{a}{c}
floor}{6}+frac{n(n+1)lfloorfrac{b}{c}
floor}{2}
]
(a<c) 且 (b<c)
[H(a,b,c,n)=sum_i^n icdot f(i)
]
[= sum_{i=0}^nsum_{j=0}^{f(i)-1} i
]
[= sum_{j=0}^{f(n) -1}sum_{i=0}^{n}icdot [j<lfloorfrac{ai+b}{c}
floor]
]
[=sum_0^{f(n)-1}sum_{i=0}^n icdot [i>lfloorfrac{cj-b+c-1}{a}
floor]
]
设(f'(i)=lfloorfrac{ci-b+c-1}{a} floor)
[H(a,b,c,n)= sum_{j=0}^{f(n)-1 }sum_{i=f'(j)+1}^{n}i
]
[= sum_{j=0}^{f(n)-1 }frac{(f'(j)+1+n)(n-f'(j))}{2}
]
[= sum_{j=0}^{f(n)-1 } frac{-f'(j)^2-f'(j)+n(n+1)}{2}
]
(H(a,b,c,n)=frac{f(n)n(n+1)-G(c,-b+c-1,a,f(n)-1)-F(c,-b+c-1,a,f(n)-1)}{2})
Tip:
这个多个函数的情况，如果直接递归写复杂度极高
所以需要把访问到的所有状态存下来递推，就能保证复杂度(O(log n))
#include<bits/stdc++.h>
using namespace std;
#define Mod1(x) ((x>=P)&&(x-=P))
#define Mod2(x) ((x<0)&&(x+=P))
#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)
#define pb push_back
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }
char IO;
int rd(){
int s=0;
int f=0;
while(!isdigit(IO=getchar())) if(IO=='-') f=1;
do s=(s<<1)+(s<<3)+(IO^'0');
while(isdigit(IO=getchar()));
return f?-s:s;
}
const int N=1010,P=998244353;
ll qpow(ll x,ll k=P-2) {
ll res=1;
for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
return res;
}
const ll Inv2=(P+1)/2,Inv6=qpow(6);
//ll F(ll a,ll b,ll c,ll n);
//ll G(ll a,ll b,ll c,ll n);
//ll H(ll a,ll b,ll c,ll n);
ll D2(ll n){
n%=P;
return n*(n+1)/2%P;
}
ll D3(ll n){
n%=P;
return n*(n+1)%P*(2*n+1)%P*Inv6%P;
}
ll sa[N],sb[N],sc[N],sn[N];
ll F[N],G[N],H[N];
int cnt;
void PreCalc(ll a,ll b,ll c,ll n){
sa[++cnt]=a; sb[cnt]=b;
sc[cnt]=c; sn[cnt]=n;
if(a==0) return;
if(a>=c || b>=c) PreCalc(a%c,b%c,c,n);
else {
ll t=(a*n+b)/c;
PreCalc(c,-b+c-1,a,t-1);
}
}
/*ll F(ll a,ll b,ll c,ll n){
if(a==0) return (b/c)%P*((n+1)%P)%P;
if(a>=c || b>=c) {
ll ans=(F(a%c,b%c,c,n)+D2(n)*(a/c)%P+(n+1)%P*(b/c))%P;
ans=(ans%P+P)%P;
return ans;
}
ll t=(a*n+b)/c;
ll ans=t%P*n%P-F(c,-b+c-1,a,t-1);
ans=(ans%P+P)%P;
return ans;
}
ll G(ll a,ll b,ll c,ll n){
if(a==0) return (b/c)*(b/c)%P*(n+1)%P;
if(a>=c || b>=c) {
ll x=a/c%P,y=b/c%P;
ll ans=(G(a%c,b%c,c,n)+2*x*H(a%c,b%c,c,n)+2*y*F(a%c,b%c,c,n))%P;
ans=(ans+D3(n)*x%P*x%P+(n+1)%P*y%P*y%P+2*D2(n)*x%P*y%P)%P;
ans=(ans%P+P)%P;
return ans;
}
ll t=(a*n+b)/c;
ll ans=(t%P)*(t%P)%P*(n%P)%P-2*H(c,-b+c-1,a,t-1)-F(c,-b+c-1,a,t-1)%P;
ans=(ans%P+P)%P;
return ans;
}
ll H(ll a,ll b,ll c,ll n){
if(a==0) return D2(n)%P*(b/c)%P;
if(a>=c || b>=c) {
ll ans=(H(a%c,b%c,c,n)+(a/c)*D3(n)+D2(n)*(b/c))%P;
ans=(ans%P+P)%P;
return ans;
}
ll t=(a*n+b)/c;
ll ans=(t%P*D2(n)%P*2%P-G(c,-b+c-1,a,t-1)-F(c,-b+c-1,a,t-1))%P;
ans=(ans*Inv2%P+P)%P;
return ans;
}*/
int main(){
rep(kase,1,rd()) {
int n=rd(),a=rd(),b=rd(),c=rd();
cnt=0;
PreCalc(a,b,c,n);
F[cnt+1]=G[cnt+1]=H[cnt+1]=0;
drep(i,cnt,1) {
ll a=sa[i],b=sb[i],c=sc[i],n=sn[i];
if(a==0) F[i]=(b/c)%P*((n+1)%P)%P;
else if(a>=c || b>=c) {
ll ans=(F[i+1]+D2(n)*(a/c)%P+(n+1)%P*(b/c))%P;
ans=(ans%P+P)%P;
F[i]=ans;
} else {
ll t=(a*n+b)/c;
ll ans=t%P*n%P-F[i+1];
ans=(ans%P+P)%P;
F[i]=ans;
}
if(a==0) G[i]=(b/c)*(b/c)%P*(n+1)%P;
else if(a>=c || b>=c) {
ll x=a/c%P,y=b/c%P;
ll ans=(G[i+1]+2*x*H[i+1]+2*y*F[i+1])%P;
ans=(ans+D3(n)*x%P*x%P+(n+1)%P*y%P*y%P+2*D2(n)*x%P*y%P)%P;
ans=(ans%P+P)%P;
G[i]=ans;
} else {
ll t=(a*n+b)/c;
ll ans=(t%P)*(t%P)%P*(n%P)%P-2*H[i+1]-F[i+1];
ans=(ans%P+P)%P;
G[i]=ans;
}
if(a==0) H[i]=D2(n)%P*(b/c)%P;
if(a>=c || b>=c) {
ll ans=(H[i+1]+(a/c)*D3(n)+D2(n)*(b/c))%P;
ans=(ans%P+P)%P;
H[i]=ans;
} else {
ll t=(a*n+b)/c;
ll ans=(t%P*D2(n)%P*2%P-G[i+1]-F[i+1])%P;
ans=(ans*Inv2%P+P)%P;
H[i]=ans;
}
}
printf("%lld %lld %lld
",F[1],G[1],H[1]);
}
}