类欧几里得算法
类欧几里得算法
对于给定的元\(a,b,c,n\)
设\(f(i)=\lfloor\frac{ai+b}{c}\rfloor\)
求
\[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^ni\cdot f(i)
\]
Part1
\(a\ge c\) 或 \(b\ge c\)
\[\lfloor\frac{ai+b}{c}\rfloor=\lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor+i\lfloor \frac{a}{c}\rfloor +\lfloor \frac{b}{c}\rfloor
\]
\[F(a,b,c,n)=F(a\mod c,b\mod c,c,n)+\frac{n(n+1)\lfloor\frac{a}{c}\rfloor}{2}+(n+1)\lfloor\frac{b}{c}\rfloor
\]
\(a<c\) 且 \(b<c\)
\[F(a,b,c,n) =\sum_0^n f(i)
\]
\[= \sum_{i=0}^n\sum_{j=0}^{f(i)-1} 1
\]
\[= \sum_{i=0}^n\sum_{j=0}^{\infty}[j<f(i)]
\]
\[= \sum_{j=0}^{f(n) -1}\sum_{i=0}^{n}[j<f(i)]
\]
其中
\[[j<\lfloor\frac{ai+b}{c}\rfloor] =[j<\lceil\frac{ai+b-c+1}{c}\rceil]
\]
\[=[cj<ai+b-c+1]
\]
\[=[ai>cj-b+c-1]
\]
\[=[i>\lfloor\frac{cj-b+c-1}{a}\rfloor]
\]
\[F(a,b,c,n)= \sum_{j=0}^{f(n)-1 }\sum_{i=0}^{n}[i>\lfloor\frac{cj-b+c-1}{a}\rfloor]
\]
\[=\sum_{j=0}^{f(n)-1} n-\lfloor\frac{cj-b+c-1}{a}\rfloor
\]
\[F(a,b,c,n)=n\cdot 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\)
\(a\ge c\) 或 \(b\ge c\)
设\(\lfloor \frac{a}{c}\rfloor=x,\lfloor \frac{b}{c}\rfloor=y\)
\[\lfloor\frac{ai+b}{c}\rfloor^2=(\lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor+ix+y)^2
\]
\[=\lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor^2+2\lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor\cdot (ix+y)+(ix+y)^2
\]
\[G(a,b,c,n)=
\]
\[G(a\mod c,b\mod c,c,n)+2x H(a\mod c,b\mod c,n)+2y F(a\mod c,b\mod 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}^n\frac{f(i)(f(i)-1)}{2}+\frac{f(i)}{2}
\]
\[=2 \sum_{i=0}^n\sum_{j=0}^{f(i)-1}j+\frac{1}{2}
\]
\[= 2\sum_{i=0}^n\sum_{j=0}^{\infty}(j+\frac{1}{2})[j<f(i)]
\]
\[= 2\sum_{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>\lfloor\frac{cj-b+c-1}{a}\rfloor]
\]
\[=2\sum_{j=0}^{f(n)-1} (j+\frac{1}{2})(n-\lfloor\frac{cj-b+c-1}{a}\rfloor)
\]
\[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^ni\cdot f(i)\)
\(a\ge c\) 或 \(b\ge c\)
\[i\cdot \lfloor\frac{ai+b}{c}\rfloor=i\cdot (\lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor+i\lfloor \frac{a}{c}\rfloor +\lfloor \frac{b}{c}\rfloor)
\]
\[i\cdot \lfloor\frac{ai+b}{c}\rfloor=i\cdot \lfloor\frac{(a \mod c)i+(b \mod c)}{c}\rfloor+i^2\lfloor \frac{a}{c}\rfloor +i\lfloor \frac{b}{c}\rfloor
\]
\[H(a,b,c,n)=H(a\mod c,b\mod c,c,n)+\frac{n(n+1)(2n+1)\lfloor\frac{a}{c}\rfloor}{6}+\frac{n(n+1)\lfloor\frac{b}{c}\rfloor}{2}
\]
\(a<c\) 且 \(b<c\)
\[H(a,b,c,n)=\sum_i^n i\cdot f(i)
\]
\[= \sum_{i=0}^n\sum_{j=0}^{f(i)-1} i
\]
\[= \sum_{j=0}^{f(n) -1}\sum_{i=0}^{n}i\cdot [j<\lfloor\frac{ai+b}{c}\rfloor]
\]
\[=\sum_0^{f(n)-1}\sum_{i=0}^n i\cdot [i>\lfloor\frac{cj-b+c-1}{a}\rfloor]
\]
设\(f'(i)=\lfloor\frac{ci-b+c-1}{a}\rfloor\)
\[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\n",F[1],G[1],H[1]);
}
}