Skip to content

类欧几里得求解问题有如下形式:f(a,b,c,n)=∑i=0n−1⌊ai+bc⌋f(a,b,c,n)=\sum_{i=0}^{n-1}\lfloor \frac{ai+b}{c} \rfloor

普通floor_sum ​

若a≥ca \geq c或者b≥cb \geq c,令a′=amod  c,b′=bmod  ca'=a \mod c,b'=b \mod c

f(a,b,c,n)=∑i=0n−1⌊ac⌋i+⌊bc⌋+⌊a′i+b′c⌋=n(n−1)2⌊ac⌋+n⌊bc⌋+f(a′,b′,c,n)f(a,b,c,n)=\sum_{i=0}^{n-1} \lfloor \frac{a}{c} \rfloor i+\lfloor \frac{b}{c} \rfloor+\lfloor \frac{a'i+b'}{c} \rfloor=\frac{n(n-1)}{2}\lfloor \frac{a}{c} \rfloor+n\lfloor \frac{b}{c} \rfloor+f(a',b',c,n)

就转化成了a<c,b<ca \lt c,b \lt c的情况

考虑a<c,b<ca \lt c,b \lt c的情况:

f(a,b,c,n)=∑i=0n−1∑j=0⌊ai+bc⌋−1[1]f(a,b,c,n)=\sum_{i=0}^{n-1}\sum_{j=0}^{\lfloor\frac{ai+b}{c}\rfloor-1}[1]

令m=⌊a(n−1)+bc⌋−1m=\lfloor\frac{a(n-1)+b}{c} \rfloor-1

f(a,b,c,n)=∑j=0m∑i=0n−1[j<⌊ai+bc⌋]=∑j=0m−1∑i=0n−1[c(j+1)≤ai+b]=∑j=0m−1∑i=0n−1[i>⌊cj+c−b−1a⌋]=∑j=0m−1[n−1−⌊cj+c−b−1a⌋]=(n−1)m−f(c,c−b−1,a,m)f(a,b,c,n)=\sum_{j=0}^m\sum_{i=0}^{n-1}[j \lt \lfloor \frac{ai+b}{c} \rfloor]\\ =\sum_{j=0}^{m-1}\sum_{i=0}^{n-1}[c(j+1) \leq ai+b]\\ =\sum_{j=0}^{m-1}\sum_{i=0}^{n-1}[i \gt \lfloor \frac{cj+c-b-1}{a}\rfloor]\\ =\sum_{j=0}^{m-1}[n-1-\lfloor \frac{cj+c-b-1}{a}\rfloor]\\ =(n-1)m-f(c,c-b-1,a,m)

复杂度为O(log⁡max⁡(a,c))O(\log \max(a,c))

从几何意义讲,表示的是y=ax+bcy=\frac{ax+b}{c}下的整数点

类似地,g(a,b,c,n)=∑i=0n−1i⌊ai+bc⌋,h(a,b,c,n)=∑i=0n−1⌊ai+bc⌋2g(a,b,c,n)=\sum_{i=0}^{n-1}i\lfloor \frac{ai+b}{c} \rfloor,h(a,b,c,n)=\sum_{i=0}^{n-1}\lfloor \frac{ai+b}{c} \rfloor^2

a≥c,b≥ca \geq c,b \geq c的情况:

g(a,b,c,n)=n(n−1)(2n−1)6⌊ac⌋+n(n−1)2⌊bc⌋+g(a′,b′,c,n)g(a,b,c,n)=\frac{n(n-1)(2n-1)}{6}\lfloor \frac{a}{c} \rfloor+\frac{n(n-1)}{2}\lfloor\frac{b}{c}\rfloor+g(a',b',c,n)

h(a,b,c,n)=h(a′,b′,c,n)+2⌊ac⌋g(a′,b′,c,n)+2⌊bc⌋f(a′,b′,c,n)+n(n−1)(2n−1)6⌊ac⌋2+n(n−1)⌊ac⌋⌊bc⌋+n⌊bc⌋2h(a,b,c,n)=h(a',b',c,n)+2\lfloor \frac{a}{c} \rfloor g(a',b',c,n)+2\lfloor\frac{b}{c}\rfloor f(a',b',c,n)+\frac{n(n-1)(2n-1)}{6}\lfloor \frac{a}{c} \rfloor^2+n(n-1)\lfloor \frac{a}{c} \rfloor \lfloor \frac{b}{c} \rfloor+n\lfloor\frac{b}{c}\rfloor^2

a<c,b<ca \lt c,b \lt c的情况:

g(a,b,c,n)=12(mn(n−1)−h(c,c−b−1,a,m)−f(c,c−b−1,a,m))g(a,b,c,n)=\frac{1}{2}(mn(n-1)-h(c,c-b-1,a,m)-f(c,c-b-1,a,m))

h(a,b,c,n)=(n−1)m2−2g(c,c−b−1,a,m)−f(c,c−b−1,a,m)h(a,b,c,n)=(n-1)m^2-2g(c,c-b-1,a,m)-f(c,c-b-1,a,m)

洛谷P5170 【模板】类欧几里得算法 ​

cpp
struct FloorSum{
    static constexpr int MOD=998244353,inv2=499122177,inv6=166374059;
    struct Node{
        int f=0,g=0,h=0;
    };
    int mo(int x){return x>=MOD?x-MOD:x;}
    int sqr(int x){return 1ll*x*x%MOD;}
    Node f(int a,int b,int c,int n){
        Node res;
        if(!a){
            res.f=1ll*b/c*n%MOD;
            res.g=1ll*n*(n-1)/2%MOD*(b/c)%MOD;
            res.h=1ll*sqr(b/c)*n%MOD;
        }
        else if(a>=c||b>=c){
            Node o=f(a%c,b%c,c,n);
            res.f=(o.f+1ll*n*(n-1)/2%MOD*(a/c)+1ll*b/c*n)%MOD;
            res.g=(o.g+1ll*a/c*n%MOD*(n-1)%MOD*(n*2-1)%MOD*inv6+1ll*n*(n-1)/2%MOD*(b/c))%MOD;
            res.h=(o.h+2ll*(a/c)*o.g+2ll*(b/c)*o.f+1ll*n*(n-1)%MOD*(n*2-1)%MOD*inv6%MOD*sqr(a/c)+1ll*n*(n-1)%MOD*(a/c)%MOD*(b/c)+1ll*n*sqr(b/c))%MOD;
        }
        else{
            int m=(1ll*a*(n-1)+b)/c;
            Node o=f(c,c-b-1,a,m);
            res.f=(1ll*(n-1)*m-o.f+MOD)%MOD;
            res.g=((1ll*m*n%MOD*(n-1)-o.h-o.f)%MOD+MOD)*inv2%MOD;
            res.h=(1ll*(n-1)*m%MOD*m-mo(o.g*2)-o.f+MOD+MOD)%MOD;
        }
        return res;
    }
};
void solve(){
    FloorSum f;
    int a,b,c,n;
    cin>>n>>a>>b>>c;
    auto res=f.f(a,b,c,n+1);
    cout<<res.f<<" "<<res.h<<" "<<res.g<<"\n";
}

万能欧几里得算法 ​

具体内容参考这篇博客:https://www.cnblogs.com/apjifengc/p/17492207.html

下面只讲做法推导 记如果经过网格的一条横线,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';
}

其他例题 ​

CF1098E ​

给定长度为nn的序列aia_i,计算出每个gcd⁡(al,....ar)\gcd(a_l,....a_r)放进bb数组,从大到小排序;再计算∑b[l,...r]\sum b[l,...r]放进数组cc,然后求出数组cc的中位数,1≤n≤5∗104,1≤ai≤1051 \leq n \leq 5*10^4, 1\leq a_i \leq 10^5

思路: 首先计算bb数组,以ii为右端点,向左gcd最多迭代log⁡\log次,st表预处理区间gcd,然后二分跳log⁡\log次,处理出bb每个数字的出现次数,这部分为O(nlog⁡3n)O(n \log^3 n) 接下来二分答案mid,然后是推式子

  • 如果区间都是相同的数,计算出最多容纳的数字mx=min(cnt[i],mid/nums[i]),然后就是等差数列求和ans+=mx*(2*cnt[i]-mx+1)/2;
  • 然后是有不同的数字的情况,假设当前区间为[l,r][l,r],那么记sum=b[l+1,.....r−1]sum=b[l+1,.....r-1],记A=numsl,SA=cntl,B=numsr,SB=cntrA=nums_l,SA=cnt_l,B=nums_r,SB=cnt_r
  • 如果A∗SA+sum+B∗SB≤midA*SA+sum+B*SB \leq mid那么SASA个左端点和SBSB个右端点任选,答案为SA∗SBSA*SB
  • 否则答案为∑i=1SA∑j=1SBsum+Ai+Bj≤mid\sum_{i=1}^{SA}\sum_{j=1}^{SB}sum+Ai+Bj \leq mid

j≤mid−sum−AiB∑i=1SA∑j=1SBsum+Ai+Bj≤mid=∑i=1SAmin⁡(max⁡(0,mid−sum−AiB),SB)因为min⁡(a,b)=a−max(0,a−b)=∑i=1SAmax⁡(0,mid−sum−AiB)−∑i=1SAmax⁡(0,max⁡(0,mid−sum−AiB)−SB)=∑i=1SAmax⁡(0,mid−sum−AiB)−∑i=1SAmax⁡(0,max⁡(−SB,mid−sum−Ai−B∗SBB))=∑i=1SAmax⁡(0,mid−sum−AiB)−∑i=1SAmax⁡(0,mid−sum−Ai−B∗SBB)j \leq \frac{mid-sum-Ai}{B}\\ \sum_{i=1}^{SA}\sum_{j=1}^{SB}sum+Ai+Bj \leq mid=\sum_{i=1}^{SA}\min(\max(0, \frac{mid-sum-Ai}{B}),SB)\\ \text{因为}\min(a,b)=a-max(0,a-b)\\ =\sum_{i=1}^{SA}\max(0, \frac{mid-sum-Ai}{B})-\sum_{i=1}^{SA}\max(0,\max(0, \frac{mid-sum-Ai}{B})-SB)=\\ \sum_{i=1}^{SA}\max(0, \frac{mid-sum-Ai}{B})-\sum_{i=1}^{SA}\max(0,\max(-SB, \frac{mid-sum-Ai-B*SB}{B}))\\ =\sum_{i=1}^{SA}\max(0, \frac{mid-sum-Ai}{B})-\sum_{i=1}^{SA}\max(0, \frac{mid-sum-Ai-B*SB}{B})

然后就可以套板子了 可以双指针维护左右端点

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(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);
}
template<typename T,typename F>
struct SpraseTable{
    int n;
    vector<vector<T>> st;
    F op;
    void init(int n,const vector<T>& a){
        this->n=n;
        st.assign(n+1,vector<T>(__lg(n)+1));
        for(int i=1;i<=n;i++)st[i][0]=a[i];
        for(int j=1;1<<j<=n;j++){
            for(int i=1;i+(1<<j)-1<=n;i++){
                st[i][j]=op(st[i][j-1],st[i+(1<<j-1)][j-1]);
            }
        }
    }
    T query(int l,int r){
        int k=__lg(r-l+1);
        return op(st[l][k],st[r-(1<<k)+1][k]);
    }
};
SpraseTable<int,decltype([](int a,int b){return __gcd(a,b);})> st;
void solve(){
	int n;
	cin>>n;
	vector<int> a(n+1);
	for(int i=1;i<=n;i++)cin>>a[i];
	st.init(n,a);
	map<int,int> mp;
	for(int i=1;i<=n;i++){
		int x=i;
		while(x>=1){
			int l=1,r=x;
			int g=st.query(x,i);
			while(l<r){
				int mid=l+r>>1;
				if(st.query(mid,i)==g)r=mid;
				else l=mid+1;
			}
			mp[g]+=x-l+1;
			x=l-1;
		}
	}
	vector<int> nums,cnt;
	for(auto &[x,y]:mp){
		nums.push_back(x);
		cnt.push_back(y);
	}
	int C=n*(n+1)/2;
	C=C*(C+1)/2;
	if(C&1)C=(C+1)/2;
	else C=C/2;
	auto check=[&](int mid){
		int ans=0;
		for(int i=0;i<nums.size();i++){
			int mx=min(cnt[i],mid/nums[i]);
			ans+=mx*(2*cnt[i]-mx+1)/2;
		}
		if(nums.size()<=1)return ans>=C;
		auto get=[&](int l,int r,int sum){
			if(r>=nums.size())return 0ll;
			int a=nums[l],b=nums[r],sa=cnt[l],sb=cnt[r];
			if(a*sa+sum+b*sb<=mid)return sa*sb;
			return FloorSum(-a,mid-sum,b,sa)-max(0ll,(mid-sum)/b)-FloorSum(-a,mid-sum-b*sb,b,sa-1)+max(0ll,(mid-sum-b*sb)/b);
		};
		for(int i=0,j=1,k=0,sum=0,tot=0;i<nums.size();i++){
			if(nums[i]>mid)break;
			if(i>=j)j=i+1,sum=tot=0;
			ans+=tot*cnt[i];
			while(j<nums.size()&&sum+cnt[j]*nums[j]<=mid){
				ans+=get(i,j,sum);
				sum+=cnt[j]*nums[j];
				tot+=cnt[j];
				j++;
			}
			if(j<nums.size())ans+=get(i,j,sum);
			if(i+1<j){
				sum-=cnt[i+1]*nums[i+1];
				tot-=cnt[i+1];
			}
		}
		return ans>=C;
	};
	int l=1,r=1e15;
	while(l<r){
		int mid=l+r>>1;
		if(check(mid))r=mid;
		else l=mid+1;
	}
	cout<<l<<"\n";
}

CF1182F ​

题意: 求[a,b][a,b]内f(x)=∣sin⁡pqπx∣f(x)=|\sin \frac{p}{q}\pi x|的最大值对应的最小整数xx,0≤a,b≤109,1≤p,q≤1090 \leq a,b \leq 10^9,1 \leq p,q \leq 10^9

思路:

  • pxmod  qq\frac{px \mod q}{q}接近12\frac{1}{2}更优,二分这个距离,即2pxmod  2q∈[q−mid,q+mid]2px \mod 2q \in [q-mid,q+mid]
  • 结论:[pxmod  q∈[l,r]]=[⌊px−lq⌋−⌊px−r−1q⌋>0][px \mod q \in [l,r]]=[\lfloor \frac{px-l}{q} \rfloor-\lfloor \frac{px-r-1}{q} \rfloor \gt 0]
  • 因此只需求和判断是否存在即可
  • 求出最小距离后,用exgcd还原求解
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);
}
int exgcd(int a,int b,int &x,int &y){
	if(!b){
		x=1,y=0;
		return a;
	}
	int g=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return g;
}
void solve(){
	int a,b,p,q;
	cin>>a>>b>>p>>q;
	int P=p<<1,Q=q<<1;
	auto check=[&](int mid){
		int l=q-mid,r=q+mid;
		return FloorSum(P,-l,Q,b)-FloorSum(P,-l,Q,a-1)-FloorSum(P,-r-1,Q,b)+FloorSum(P,-r-1,Q,a-1)>0;
	};
	int l=0,r=1e9;
	while(l<r){
		int mid=l+r>>1;
		if(check(mid))r=mid;
		else l=mid+1;
	}
	int x,y;
	int g=exgcd(P,Q,x,y);
	int ans=1e9;
	int dx=Q/g,dy=P/g;
	if((q-l)%g==0){
		int val=(q-l)/g*x;
		val+=(a-val)/dx*dx;
		while(val>=a)val-=dx;
		while(val<a)val+=dx;
		ans=min(ans,val);
	}
	if((q+l)%g==0){
		int val=(q+l)/g*x;
		val+=(a-val)/dx*dx;
		while(val>=a)val-=dx;
		while(val<a)val+=dx;
		ans=min(ans,val);
	}
	cout<<ans<<"\n";
}

abc313G ​

abc402G ​

2020ICPC沈阳I ​

在这里插入图片描述

∑i=0H−1∑j=0M−1[∣2π(iM+j)HM−2πjM∣≤2πAHM]∑i=0H−1∑j=0M−1[∣iM+j−Hj∣≤A]HM−(∑i=0H−1∑j=0M−1[iM+j−Hj>A]+∑i=0H−1∑j=0M−1[iM+j−Hj<−A])\sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[|\frac{2 \pi(iM+j)}{HM}-\frac{2 \pi j}{M}| \leq \frac{2 \pi A}{HM}]\\ \sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[|iM+j-Hj| \leq A]\\ HM-(\sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[iM+j-Hj\gt A]+\sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[iM+j-Hj\lt -A])

对于前半部分:

∑i=0H−1∑j=0M−1[iM+j−Hj>A]=∑i=0H−1∑j=0M−1[i>⌊A+(H−1)jM⌋]=∑j=0M−1H−1−⌊A+(H−1)jM⌋=M(H−1)−∑j=0M−1⌊A+(H−1)jM⌋\sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[iM+j-Hj\gt A]=\\ \sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[i \gt \lfloor\frac{A+(H-1)j}{M}\rfloor]=\\ \sum_{j=0}^{M-1}H-1-\lfloor\frac{A+(H-1)j}{M}\rfloor=\\ M(H-1)-\sum_{j=0}^{M-1}\lfloor\frac{A+(H-1)j}{M}\rfloor

对于后半部分:

∑i=0H−1∑j=0M−1[iM+j−Hj<−A]=∑i=0H−1∑j=0M−1[i<⌊(H−1)j−AM⌋]=∑i=0H−1∑j=0M−1[i≤⌊(H−1)j−A−1M⌋]=∑j=0M−1⌊(H−1)j−A−1M⌋+1=∑j=0M−1⌊(H−1)j−A−1+MM⌋\sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[iM+j-Hj\lt -A]=\\ \sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[i \lt \lfloor\frac{(H-1)j-A}{M}\rfloor]=\\ \sum_{i=0}^{H-1}\sum_{j=0}^{M-1}[i \leq \lfloor\frac{(H-1)j-A-1}{M}\rfloor]=\\ \sum_{j=0}^{M-1} \lfloor\frac{(H-1)j-A-1}{M}\rfloor+1=\\ \sum_{j=0}^{M-1} \lfloor\frac{(H-1)j-A-1+M}{M}\rfloor

似乎要计算负数部分和特判HM

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);
}
void solve(){
    int H,M,A;
    cin>>H>>M>>A;
    int ans=H*M;
    ans-=M*(H-1)-FloorSum(H-1,A,M,M-1);
    ans-=FloorSum(H-1,M-A-1,M,M-1);
    cout<<min(H*M,ans)<<"\n";
}

补一发易理解的正解,设时间为TT,时针角速度ω1=2πHM\omega_1=\frac{2\pi}{HM},分针角速度ω2=2πM\omega_2=\frac{2\pi}{M}

因此θ=T2π(H−1)HMmod  2π\theta=T\frac{2 \pi (H-1)}{HM} \mod 2\pi,min⁡(θ,2π−θ)≤2πAHM\min(\theta,2\pi-\theta) \leq \frac{2 \pi A}{HM}

答案就是[T(H−1)mod  HM≤A]+[T(H−1)mod  HM≥HM−A][T(H-1) \mod HM \leq A] + [T(H-1) \mod HM \geq HM-A]

当A=HM2A=\frac{HM}{2}时会有重合,此时相当于角度小于等于π\pi,则任意时刻都满足,答案为HMHM

否则来看一下式子T(H−1)mod  HM≤AT(H-1) \mod HM \leq A,令G=gcd⁡(H−1,HM)G=\gcd(H-1,HM)

如果G=1G=1,因为T∈[0,HM−1]T \in [0,HM-1]是HMHM的完全剩余系,所以T(H−1)T(H-1)也是HMHM的完全剩余系,答案就是A+1A+1

否则,考虑TH−1Gmod  HMG≤⌊AG⌋T\frac{H-1}{G} \mod \frac{HM}{G} \leq \lfloor \frac{A}{G} \rfloor,变成上面的情况,答案就是⌊AG⌋+1\lfloor \frac{A}{G} \rfloor+1,相当于将[0,HM−1][0,HM-1]切分成GG段,因此总答案为G(⌊AG⌋+1)G(\lfloor \frac{A}{G} \rfloor+1)

考虑T(H−1)mod  HM≥HM−AT(H-1) \mod HM \geq HM-A

TH−1Gmod  HMG≥⌈HM−AG⌉T \frac{H-1}{G} \mod \frac{HM}{G} \geq \lceil \frac{HM-A}{G} \rceil,答案就是G(HM−1G−⌊HM−A+G−1G⌋+1)G(\frac{HM-1}{G}-\lfloor \frac{HM-A+G-1}{G} \rfloor+1)

cpp
void solve(){
    int H,M,A;
    cin>>H>>M>>A;
    if(2*A==H*M){
        cout<<H*M<<"\n";
        return;
    }
    int G=__gcd(H-1,H*M);
    cout<<G*(A/G+1+(H*M-1)/G-(H*M-A+G-1)/G+1)<<"\n";
}

LOJ6344 ​

在这里插入图片描述 按位考虑一个数第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
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);
}
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";
}