Skip to content

万能欧几里得 ​

记如果经过网格的一条横线,y=c(c是大于等于0的任意正整数)y=c(c是大于等于0的任意正整数),就乘上矩阵UU,如果经过x=c(c是大于等于0的任意正整数)x=c(c是大于等于0的任意正整数),乘上矩阵RR,如果同时经过则先UU再RR

例子:y=⌊ai+bc⌋y=\lfloor \frac{ai+b}{c} \rfloor,我们维护向量[y,∑y,1][y,\sum y,1]

则经过UU时,y+1y+1,经过RR,∑y+y\sum y+y

U=(100010101),R=(110010001)U=\begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 1 & 0 & 1 \end{pmatrix}, R=\begin{pmatrix} 1 & 1 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}

模板如下,复杂度为O(Tlog⁡max⁡(a,c))O(T \log \max(a,c))

cpp
template<typename T>
struct FloorSum{
    T qpow(T a,int b){
        T res;
        while(b){
            if(b&1)res=res*a;
            a=a*a;
            b>>=1;
        }
        return res;
    }
    T f(int a,int b,int c,int n,T U,T R){
        if(!n)return T();
        if(b>=c)return qpow(U,b/c)*f(a,b%c,c,n,U,R);
        if(a>=c)return f(a%c,b,c,n,U,qpow(U,a/c)*R);
        int m=((__int128)a*n+b)/c;
        if(!m)return qpow(R,n);
        return qpow(R,(c-b-1)/a)*U*f(c,(c-b-1)%a,a,m-1,R,U)*qpow(R,n-((__int128)c*m-b-1)/a);
    }
};

考虑模板题,要求的是∑y,∑y2,∑xy\sum y,\sum y^2 ,\sum xy,如何合并两个区间

定义Xl/r,Yl/rX_{l/r},Y_{l/r}表示左/右区间的cntr,cntucntr,cntu的值(经过多少个R,UR,U)

∑y=∑yl+∑(yr+Yl)=∑yl+∑yr+XrYl∑y2=∑yl2+∑(yr+Yl)2=∑yl2+∑yr2+2Yl∑yr+XrYl2∑xy=∑xlyl+∑(xr+Xl)(yr+Yl)=∑xlyl+∑xryr+Yl∑xr+Xl∑yr+XrXlYl∑x=∑xl+∑(xr+Xl)=∑xl+∑xr+XlXr\sum y=\sum y_l+\sum (y_r+Y_l)=\sum y_l+\sum y_r+X_rY_l\\ \sum y^2=\sum y_l^2+\sum(y_r+Y_l)^2=\sum y_l^2+\sum y_r^2+2Y_l\sum y_r+X_rY_l^2\\ \sum xy=\sum x_ly_l+\sum (x_r+X_l)(y_r+Y_l)=\sum x_ly_l+\sum x_ry_r+Y_l\sum x_r+X_l\sum y_r+X_rX_lY_l\\ \sum x=\sum x_l+\sum (x_r+X_l)=\sum x_l+\sum x_r+X_lX_r

总结:

  • 左边求和的数量为XlX_l,右边求和的数量为XrX_r
  • 合并时,右边的x,yx,y加上Xl,YlX_l,Y_l
  • U,RU,R初始化时,UU初始化yy,RR初始化xx以及其他和xx相关的答案

写出如下结构体

cpp
struct Info{
    int x,y,sumx,sumy,sumxy,sumy2;
    Info():x(0),y(0),sumx(0),sumy(0),sumxy(0),sumy2(0){}
    Info operator * (const Info &r){
        Info res;
        res.x=(x+r.x)%mod;
        res.y=(y+r.y)%mod;
        res.sumx=(sumx+r.sumx+x*r.x%mod)%mod;
        res.sumy=(sumy+r.sumy+y*r.x%mod)%mod;
        res.sumy2=(sumy2+r.sumy2+2*y%mod*r.sumy%mod+r.x*y%mod*y%mod)%mod;
        res.sumxy=(sumxy+r.sumxy+y*r.sumx%mod+x*r.sumy%mod+x*y%mod*r.x%mod)%mod;
        return res;
    }
};

初始化U,RU,R,求解的是[1,n][1,n]的范围,需要加上i=0i=0的答案

cpp
void solve(){
    int n,a,b,c;
    cin>>n>>a>>b>>c;
    Info U,R;
    U.y=R.x=R.sumx=1;
    FloorSum<Info> f;
    auto ans=f.f(a,b,c,n,U,R);
    ans.sumy=(ans.sumy+b/c)%mod;
    ans.sumy2=(ans.sumy2+(b/c)*(b/c))%mod;
    cout<<ans.sumy<<" "<<ans.sumy2<<" "<<ans.sumxy<<"\n";
}

求∑i=1nAiB⌊ai+bc⌋\sum_{i=1}^{n}A^iB^{\lfloor \frac{ai+b}{c} \rfloor} ​

维护∑AxBy,AX,BY\sum A^xB^y,A^X,B^Y

∑AxBy=∑AxlByl+∑AXl+xrBYl+yr=∑AxlByl+AXl∑AxrByrBYl\sum A^xB^y=\sum A^{x_l}B^{y_l}+\sum A^{X_l+x_r}B^{Y_l+y_r}=\sum A^{x_l}B^{y_l}+A^{X_l}\sum A^{x_r}B^{y_r}B^{Y_l}

cpp
struct Matrix{
    static constexpr int maxn=22;
    int m[maxn][maxn];
    Matrix(){
        memset(m,0,sizeof m);
    }
    Matrix(int n){
        init(n);
    }
    void init(int n){
        memset(m,0,sizeof m);
        for(int i=0;i<n;i++)m[i][i]=1;
    }
    Matrix operator*(const Matrix &a){
        Matrix res;
        for(int k=0;k<maxn;k++){
            for(int i=0;i<maxn;i++){
                for(int j=0;j<maxn;j++){
                    res.m[i][j]=(res.m[i][j]+m[i][k]*a.m[k][j]%mod)%mod;
                }
            }
        }
        return res;
    }
    Matrix operator+(const Matrix &a){
        Matrix res;
        for(int i=0;i<maxn;i++){
            for(int j=0;j<maxn;j++){
                res.m[i][j]=(m[i][j]+a.m[i][j])%mod;
            }
        }
        return res;
    }
};
struct Info{
    Matrix sum,x,y;
    Info():sum(),x(22),y(22){}
    Info operator * (const Info &r){
        Info res;
        res.x=x*r.x;
        res.y=y*r.y;
        res.sum=sum+x*r.sum*y;
        return res;
    }
};
template<typename T>
struct FloorSum{
    T qpow(T a,int b){
        T res;
        while(b){
            if(b&1)res=res*a;
            a=a*a;
            b>>=1;
        }
        return res;
    }
    T f(int a,int b,int c,int n,T U,T R){
        if(!n)return T();
        if(b>=c)return qpow(U,b/c)*f(a,b%c,c,n,U,R);
        if(a>=c)return f(a%c,b,c,n,U,qpow(U,a/c)*R);
        int m=((__int128)a*n+b)/c;
        if(!m)return qpow(R,n);
        return qpow(R,(c-b-1)/a)*U*f(c,(c-b-1)%a,a,m-1,R,U)*qpow(R,n-((__int128)c*m-b-1)/a);
    }
};
void solve(){
    int p,q,r,l,n;
    cin>>p>>q>>r>>l>>n;
    Info R,U;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)cin>>R.x.m[i][j],R.sum.m[i][j]=R.x.m[i][j];
    }
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)cin>>U.y.m[i][j];
    }
    FloorSum<Info> f;
    auto ans=f.f(p,r,q,l,U,R);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)cout<<ans.sum.m[i][j]<<" \n"[j==n-1];
    }
}

求∑i=0n⌊ai+bc⌋\sum_{i=0}^n \lfloor\frac{ai+b}{c}\rfloor(简单版) ​

cpp
int FloorSum(int a,int b,int c,int n){// \sum_{i=0}^{n} (ai+b)/c
	if(!a)return b/c*(n+1);
	if(a<0){
		b+=n*a;
		a=-a;
	}
	if(b<0){
		int d=(a-b-1)/a;
		n-=d;
		b+=a*d;
        if(n<0)return 0;
    }
    //======== 负数计入答案 ======================
    // if(b<0){
    //     int k=(c-b-1)/c;
    //     return FloorSum(a,b+k*c,c,n)-k*(n+1);
	// }
    //===========================================
	if(a>=c||b>=c)return n*(n+1)/2*(a/c)+(b/c)*(n+1)+FloorSum(a%c,b%c,c,n);
	int m=(a*n+b)/c;
	return m*n-FloorSum(c,c-b-1,a,m-1);
}

求1~n除M模R的数的二进制1的数量 ​

结论:对于一个数AA二进制第kk位为11,有⌊A+2k2k+1⌋−⌊A2k+1⌋=1\lfloor \frac{A+2^k}{2^{k+1}}\rfloor-\lfloor \frac{A}{2^{k+1}} \rfloor=1,另一种写法是⌊A2k⌋−2⌊A2k+1⌋=1\lfloor \frac{A}{2^k} \rfloor-2\lfloor \frac{A}{2^{k+1}} \rfloor=1

可以用类欧求解

求∑x=0nxk1⌊ax+bc⌋k2\sum_{x=0}^n x^{k_1}\lfloor \frac{ax+b}{c} \rfloor^{k_2} ​

∑xk1yk2=∑xlk1ylk2+∑(xr+Xl)k1(yr+Yl)k2∑(xr+Xl)k1(yr+Yl)k2=∑∑i=0k1∑j=0k2(k1i)Xlixrk1−i(k2j)Yljyrk2−j=∑i=0k1∑j=0k2(k1i)(k2j)XliYlj∑xrk1−iyrk2−j\sum x^{k_1}y^{k_2}=\sum x_l^{k_1}y_l^{k_2}+\sum(x_r+X_l)^{k_1}(y_r+Y_l)^{k_2}\\ \sum(x_r+X_l)^{k_1}(y_r+Y_l)^{k_2}=\sum \sum_{i=0}^{k_1}\sum_{j=0}^{k_2}\binom{k_1}{i}X_l^ix_r^{k_1-i}\binom{k_2}{j}Y_l^jy_r^{k_2-j}=\sum_{i=0}^{k_1}\sum_{j=0}^{k_2}\binom{k_1}{i}\binom{k_2}{j}X_l^iY_l^j \sum x_r^{k_1-i}y_r^{k_2-j}

因此启发我们维护x,y,∑xiyjx,y,\sum x^iy^j

cpp
struct Info{
    int x,y;
    vector<vector<int>> sum;
    Info():sum(vector (k1+1,vector<int>(k2+1))),x(0),y(0){}
    Info operator * (const Info &r){
        Info res;
        res.x=(x+r.x)%mod;
        res.y=(y+r.y)%mod;
        int cx=x%mod,cy=y%mod;
        for(int u=0;u<=k1;u++){
            for(int v=0;v<=k2;v++){
                res.sum[u][v]=sum[u][v];
                for(int i=0,X=1;i<=u;i++,X=X*cx%mod){
                    for(int j=0,Y=1;j<=v;j++,Y=Y*cy%mod){
                        res.sum[u][v]=(res.sum[u][v]+comb.binom(u,i)*comb.binom(v,j)%mod*X%mod*Y%mod*r.sum[u-i][v-j]%mod)%mod;
                    }
                }
            }
        }
        return res;
    }
};
void solve(){
    int n,a,b,c;
    cin>>n>>a>>b>>c>>k1>>k2;
    FloorSum<Info> f;
    Info U,R;
    U.y=R.x=1;
    for(int i=0;i<=k1;i++)R.sum[i][0]=1;
    comb.init(15);
    auto ans=f.f(a,b,c,n,U,R);
    int res=ans.sum[k1][k2];
    if(k1==0)res=(res+comb.qpow(b/c,k2))%mod;//k1=0,0^0=1,答案记得加上
    cout<<res<<"\n";
}

求ax+by≤nax+by \leq n的对数(x,y)(x,y) ​

∑i=0na⌊n−ia+bb⌋\sum_{i=0}^{\frac{n}{a}}\lfloor \frac{n-ia+b}{b} \rfloor

cpp
void solve(){
	int n,a,b;
	cin>>n>>a>>b;
	cout<<FloorSum(-a,n+b,b,n/a)<<'\n';
}

求最简分数ab<pq<cd\frac{a}{b} \lt \frac{p}{q} \lt \frac{c}{d} ​

cpp
// P5179
void work(int a,int b,int c,int d,int &p,int &q){
    if(a<b&&c>d){
        p=q=1;
        return;
    }
    work(d%c,c,b-(d/c)*a,a,q,p);
    q+=d/c*p;
}
void solve(){
    int a,b,c,d;
    while(cin>>a>>b>>c>>d){
        int p,q;
        work(a,b,c,d,p,q);
        cout<<p<<"/"<<q<<"\n";
    }
}

求(nmod  1)⊕(nmod  2)⊕...⊕(nmod  n)(n \mod 1) \oplus (n \mod 2)\oplus...\oplus(n \mod n) ​

按位考虑一个数第kk位为11,等价于⌊x2k⌋−2⌊x2k+1⌋\lfloor \frac{x}{2^k} \rfloor -2\lfloor \frac{x}{2^{k+1}} \rfloor

定义ansians_i表示∑i=1n⌊n−i⌊ni⌋2k⌋\sum_{i=1}^n \lfloor \frac{n-i\lfloor \frac{n}{i} \rfloor}{2^k}\rfloor

对于小于等于n\sqrt n的数暴力做,这些块比较密集;大于的用整除分块+类欧

cpp
void solve(){
    int n;
    cin>>n;
    int c=0;
    for(int i=1;i<=min(n,(int)1e6);i++)c^=n%i;
    vector<int> ans(50);
    for(int l=1e6+1,r;l<=n;l=r+1){
        int v=n/l;
        r=n/(n/l);
        for(int i=0;i<40;i++){
            ans[i]+=FloorSum(-v,n-l*v,1ll<<i,r-l);
        }
    }
    int res=0;
    for(int i=0;i<39;i++){
        res+=((ans[i]-2*ans[i+1])&1)*(1ll<<i);
    }
    cout<<(c^res)<<"\n";
}