洛谷P5170 【模板】类欧几里得算法（数论）

传送门

此题剧毒，公式恐惧症患者请直接转去代码→_→

前置芝士

基本数论芝士

题解

本题就是要我们求三个函数的值

[f(a,b,c,n)=sum_{i=0}^n leftlfloorfrac{ai+b}{c} ight floor ]

[g(a,b,c,n)=sum_{i=0}^n leftlfloorfrac{ai+b}{c} ight floor^2 ]

[h(a,b,c,n)=sum_{i=0}^n ileftlfloorfrac{ai+b}{c} ight floor ]

先考虑(f)，如果(a)为(0)，那么显然$$f(a,b,c,n)=(n+1)frac{b}{c}$$

根据基本数论芝士，有$$leftlfloorfrac{ax}{b} ight floor=leftlfloorfrac{a(x% b)}{b} ight floor+aleftlfloorfrac{x}{b} ight floor$$

所以，当(ageq c)或(bgeq c)时，有$$f(a,b,c,n)=sum_{i=0}^n leftlfloorfrac{(amod c)i+(bmod c)}{c} ight floor+ileftlfloorfrac{a}{c} ight floor+leftlfloorfrac{b}{c} ight floor$$

[f(a,b,c,n)=f(amod c,bmod c,c,n)+frac{n(n+1)}{2}leftlfloorfrac{a}{c} ight floor+(n+1)leftlfloorfrac{b}{c} ight floor ]

然后考虑(a<c)且(b<c)的情况

[f(a,b,c,n)=sum_{i=0}^n leftlfloorfrac{ai+b}{c} ight floor ]

令(M=leftlfloorfrac{an+b}{c} ight floor)

[f(a,b,c,n)=sum_{i=0}^nsum_{j=1}^M left[leftlfloorfrac{ai+b}{c} ight floorgeq j ight] ]

[f(a,b,c,n)=sum_{i=0}^nsum_{j=1}^M left[ai+bgeq jc ight] ]

[f(a,b,c,n)=sum_{i=0}^nsum_{j=0}^{M-1} left[ai+b> jc+c-1 ight] ]

[f(a,b,c,n)=sum_{j=0}^{M-1}sum_{i=0}^nleft[i>leftlfloorfrac{jc+c-b-1}{a} ight floor ight] ]

[f(a,b,c,n)=sum_{j=0}^{M-1}n-leftlfloorfrac{jc+c-b-1}{a} ight floor ]

[f(a,b,c,d)=nM-f(c,c-b-1,a,M-1) ]

综上，(f)就可以递归计算了，

[f(a,b,c,n)=egin{cases}(n+1)lfloorfrac{b}{c} floor&,a=0\frac{n(n+1)}{2}lfloorfrac{a}{c} floor+(n+1)lfloorfrac{b}{c} floor+f(amod c,bmod c,c,n)&,age c or bge c\nM-f(c,c-b-1,a,M-1),M=lfloorfrac{an+b}{c} floor &,otherwiseend{cases} ]

发现这东西是((a,c)->(c,a)->(cmod a,a))，很像(gcd)的过程，所以这玩意儿的复杂度和(gcd)差不多，都是(O(log n))

接下来我们要计算(g,h)，这两个因为有关联所以我们可以放到一起算

当(a=0)时，有

[g(a,b,c,n)=(n+1)leftlfloorfrac{b}{c} ight floor^2 ]

[h(a,b,c,n)=frac{n(n+1)}{2}leftlfloorfrac{b}{c} ight floor ]

然后考虑(ageq c)或者(bgeq c)的情况，有

[g(a,b,c,n)=sum_{i=0}^nleft(leftlfloorfrac{(amod c)i+(bmod c)}{c} ight floor+ileftlfloorfrac{a}{c} ight floor+leftlfloorfrac{b}{c} ight floor ight)^2 ]

[g(a,b,c,n)=sum_{i=0}^nleftlfloorfrac{(amod c)i+(bmod c)}{c} ight floor^2+2left(ileftlfloorfrac{a}{c} ight floor+leftlfloorfrac{b}{c} ight floor ight)leftlfloorfrac{(amod c)i+(bmod c)}{c} ight floor+left(ileftlfloorfrac{a}{c} ight floor+leftlfloorfrac{b}{c} ight floor ight)^2 ]

[g(a,b,c,n)=g(amod c,bmod c,c,n)+2leftlfloorfrac{a}{c} ight floor h(amod c,bmod c,c,n)+2leftlfloorfrac{b}{c} ight floor f(amod c,bmod c,c,n)+sum_{i=0}^ni^2leftlfloorfrac{a}{c} ight floor+2ileftlfloorfrac{a}{c} ight floorleftlfloorfrac{b}{c} ight floor+leftlfloorfrac{b}{c} ight floor^2 ]

[g(a,b,c,n)=g(amod c,bmod c,c,n)+2leftlfloorfrac{a}{c} ight floor h(amod c,bmod c,c,n)+2leftlfloorfrac{b}{c} ight floor f(amod c,bmod c,c,n)+frac{n(n+1)(2n+1)}{6}leftlfloorfrac{a}{c} ight floor+n(n+1)leftlfloorfrac{a}{c} ight floorleftlfloorfrac{b}{c} ight floor+(n+1)leftlfloorfrac{b}{c} ight floor^2 ]

然后(h)的话就是

[h(a,b,c,n)=sum_{i=0}^n ileft(leftlfloorfrac{(amod c)i+(bmod c)}{c} ight floor+ileftlfloorfrac{a}{c} ight floor+leftlfloorfrac{b}{c} ight floor ight) ]

[h(a,b,c,n)=h(amod c,bmod c,c,n)+frac{n(n+1)(2n+1)}{6}leftlfloorfrac{a}{c} ight floor+frac{n(n+1)}{2}leftlfloorfrac{b}{c} ight floor ]

然后就是(a<c)且(b<c)的情况，依旧令(M=leftlfloorfrac{an+b}{c} ight floor)

发现(leftlfloorfrac{ai+b}{c} ight floor^2)不好搞，因为有(x^2=-x+2sum_{i=1}^x i)，那么就可以继续推倒了

[g(a,b,c,n)=sum_{i=0}^nleft(-leftlfloorfrac{ai+b}{c} ight floor+2sum_{j=1}^{leftlfloorfrac{ai+b}{c} ight floor} j ight) ]

[g(a,b,c,n)=-f(a,b,c,n)+2sum_{i=0}^nsum_{j=1}^M jleft[jleq leftlfloorfrac{ai+b}{c} ight floor ight] ]

后面那个东西就和算(f)的时候一样搞掉

[g(a,b,c,n)=-f(a,b,c,n)+2sum_{j=0}^{M-1}(j+1)(n-leftlfloorfrac{jc+c-b-1}{a} ight floor) ]

[g(a,b,c,n)=-f(a,b,c,n)+nm(m+1)-2h(c,c-b-1,a,M-1)-2f(c,c-b-1,a,M-1) ]

(g)就算好了，然后来算(h)

[h(a,b,c,n)=sum_{i=0}^nisum_{j=1}^M left[leftlfloorfrac{ai+b}{c} ight floorgeq j ight] ]

[h(a,b,c,n)=sum_{j=0}^{M-1}sum_{i=0}^nileft[i>leftlfloorfrac{jc+c-b-1}{a} ight floor ight] ]

[h(a,b,c,n)=sum_{j=0}^{M-1}frac{n(n+1)}{2}-sum_{i=0}^nileft[ileq leftlfloorfrac{jc+c-b-1}{a} ight floor ight] ]

[h(a,b,c,n)=sum_{j=0}^{M-1}frac{n(n+1)}{2}-frac{leftlfloorfrac{jc+c-b-1}{a} ight floor(leftlfloorfrac{jc+c-b-1}{a} ight floor+1)}{2} ]

[h(a,b,c,n)=frac{1}{2}left[nm(n+1)-g(c,c-b-1,a,M-1)-f(c,c-b-1,a,M-1) ight] ]

然后就可以了

[g(a,b,c,n)=egin{cases}(n+1)lfloorfrac{b}{c} floor^2&,a=1\g(amod c,bmod c,c,n)+2lfloorfrac{a}{c} floor h(amod c,bmod c,c,n)+\2lfloorfrac{b}{c} floor f(amod c,bmod c,c,n)+frac{n(n+1)(2n+1)}{6}lfloorfrac{a}{c} floor^2+\n(n+1)lfloorfrac{a}{c} floorlfloorfrac{b}{c} floor+(n+1)lfloorfrac{b}{c} floor^2&,age c or bge c\nM(M+1)-f(a,b,c,n)-2h(c,c-b-1,a,M-1)-\2f(c,c-b-1,a,M-1)&,otherwiseend{cases} ]

[h(a,b,c,n)=egin{cases}frac{n(n+1)}{2}lfloorfrac{b}{c} floor&,a=0\h(amod c,bmod c,c,n)+frac{n(n+1)(2n+1)}{6}lfloorfrac{a}{c} floor+frac{n(n+1)}{2}lfloorfrac{b}{c} floor&,age c or bge c\frac{1}{2}[Mn(n+1)-g(c,c-b-1,a,M-1)-f(c,c-b-1,a,M-1)]&,otherwiseend{cases} ]

递归计算就行了，顺便记得把(f,g,h)同步计算

总算写完了→_→

//minamoto
#include<bits/stdc++.h>
#define R register
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res=1,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int K=-1,Z=0;
inline void Ot(){fwrite(sr,1,K+1,stdout),K=-1;}
void print(R int x){
    if(K>1<<20)Ot();if(x<0)sr[++K]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++K]=z[Z],--Z);sr[++K]=' ';
}
const int P=998244353,inv2=499122177,inv6=166374059;
inline int add(R int x,R int y){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y){return x-y<0?x-y+P:x-y;}
inline int mul(R int x,R int y){return 1ll*x*y-1ll*x*y/P*P;}
inline int pow(R int x){return mul(x,x);}
inline int s(R int x){return 1ll*x*(x+1)%P*inv2%P;}
inline int ss(R int x){return 1ll*x*(x+1)%P*((x<<1)+1)%P*inv6%P;}
int ksm(R int x,R int y){
    R int res=1;
    for(;y;y>>=1,x=mul(x,x))if(y&1)res=mul(res,x);
    return res;
}
struct node{int f,g,h;}res;
int n;
node get(int a,int b,int c,int n){
    node res;
    int x=a/c,y=b/c;
    if(!a){
        res.f=1ll*y*(n+1)%P;
        res.g=1ll*pow(y)*(n+1)%P;
        res.h=1ll*y*s(n)%P;
        return res;
    }
    if(a>=c||b>=c){
        res=get(a%c,b%c,c,n);
        res.g=add(res.g,add(1ll*(x<<1)*res.h%P,add(1ll*(y<<1)*res.f%P,add(1ll*ss(n)*pow(x)%P,add(1ll*n*(n+1)%P*x%P*y%P,1ll*(n+1)*pow(y)%P)))));
        res.h=add(res.h,add(1ll*ss(n)*x%P,1ll*s(n)*y%P));
        res.f=add(res.f,add(1ll*s(n)*x%P,1ll*(n+1)*y%P));
        return res;
    }
    int M=(1ll*a*n+b)/c;
    res=get(c,c-b-1,a,M-1);
    int h=res.h,g=res.g,f=res.f;
    res.f=dec(1ll*n*M%P,res.f);
    res.g=dec(dec(dec(1ll*n*M%P*(M+1)%P,res.f),mul(h,2)),mul(f,2));
    res.h=1ll*inv2*dec(dec(1ll*M*n%P*(n+1)%P,g),f)%P;
    return res;
}
int main(){
//	freopen("testdata.in","r",stdin);
    int T=read();
    while(T--){
        int n=read(),a=read(),b=read(),c=read();
        res=get(a,b,c,n);
        print(res.f),print(res.g),print(res.h),sr[K]='
';
    }
    return Ot(),0;
}