Re: Eat an Ultimate Shit（[Ynoi2004] tars2 的题解）

Seauy · 2025-03-10 17:51:07 · 题解

死亡回放：15kb代码卡常不能

这次尝试点科学的办法。

我们切割修改区域，变成八个三角形区域的加法，并旋转或翻折成 X+Y>x+y-d,Y\leq y,X\leq x 的形式；查询常规地拆成四个左下角区域。每个三角形又能拆成两个向右无限延伸的三角形和两个右上角区域，设无限三角形左上角为 (x,y)，则限制为 X+Y\geq x+y,Y\leq y。

再设 (x,y) 处要加的值为 wd，查询为 (a,b)，则当 a+b\geq x+y,b\leq y 时，贡献为：

w\sum_{i=d-a+x}^{d-y+b} i(d-y+b+1-i)

不管 w，设 x'=d+x,y'=d-y，后面那坨化简为

\frac{1}{2}(x'+y'-a+b)(a+b-x'+y'+1)(y'+b+1)-\frac{1}{6}(y'+b)(y'+b+1)(2y'+2b+1)+\frac{1}{6}(x'-a-1)(x'-a)(2x'-2a-1)

我们需要以 a,b 主元，求出所有项系数关于 x',y' 的多项式（反过来也行），然后 CDQ 分治。如果不想暴力拆完，可以取若干点后高斯消元求出所有 a^ib^jx'^ky'^l\ (0\leq i+j+k+l\leq 3) 的系数，提个 \frac{1}{6} 的系数出来就全是整数。下面给个生成系数的代码，因为精度问题不一定每次都能算出正常的解，并且注意到实际上我求的是 (a+b)^i(b-2a)^jx'^ky'^l 的系数，因为会比较好看：

#include<bits/stdc++.h>
using namespace std;

const int MAXN=100;
const double eps=1e-5;

int n;
double A[MAXN+5][MAXN+5];

int f(int a,int b,int x,int y)
{
    int s1=(x+y-a+b)*(a+b-x+y+1)*(y+b+1)/2;
    int s2=(y+b)*(y+b+1)*(2*y+2*b+1)/6;
    int s3=(x-a-1)*(x-a)*(2*x-2*a-1)/6;
    return s1-s2+s3;
}

bool Solve()
{
    for(int i=1;i<=n;i++)
    {
        int minn=0;
        for(int j=i;j<=n;j++)
            if(fabs(A[j][i])>eps && (!minn || A[j][i]<A[minn][i]))
                minn=i;
        if(!minn) return 0;
        swap(A[minn],A[i]);
        for(int j=n+1;j>=i;j--) A[i][j]/=A[i][i];
        for(int j=i+1;j<=n;j++)
            for(int k=n+1;k>=i;k--)
                A[j][k]-=A[j][i]*A[i][k];
    }
    for(int i=n;i;i--)
        for(int j=i-1;j;j--)
            A[j][n+1]-=A[j][i]*A[i][n+1];
    return 1;
}

int main()
{
    //freopen("coe.txt","w",stdout);
    srand(time(0));
    while(1)
    {
        n=35;
        for(int T=1;T<=n;T++)
        {
            int cnt=0;
            int a=rand()%10,b=rand()%10;
            int x=rand()%10,y=rand()%10;
            A[T][n+1]=f(a,b,x,y);
            for(int i=0;i<=3;i++)
                for(int j=0;j<=3;j++)
                    for(int k=0;k<=3;k++)
                        for(int l=0;l<=3;l++)
                            if(i+j+k+l<=3)
                            {
                                ++cnt;
                                A[T][cnt]=pow(a+b,i)*pow(b-2*a,j)*pow(x,k)*pow(y,l);
                            }
        }
        if(!Solve()) continue;
        int cnt=0;
        printf("\\frac{1}{6}(");
        for(int i=0;i<=3;i++)
            for(int j=0;j<=3;j++)
            {
                if(i+j>3) continue;
                printf("+(");
                bool OK=0;
                for(int k=0;k<=3;k++)
                    for(int l=0;l<=3;l++)
                        if(i+j+k+l<=3)
                        {
                            ++cnt;
                            int v=round(A[cnt][n+1]*6);
                            if(!v) continue;
                            if(v>0 && OK) printf("+%d",v);
                            else printf("%d",v),OK=1;
                            if(k>1) printf("x\'^%d",k);
                            else if(k==1) printf("x\'");
                            if(l>1) printf("y\'^%d",l);
                            else if(l==1) printf("y\'");
                        }
                printf(")");
                if(i>1) printf("(a+b)^%d",i);
                else if(i==1) printf("(a+b)");
                if(j>1) printf("(b-2a)^%d",j);
                else if(j==1) printf("(b-2a)");
                printf("\\\\\n");
            }
        break;
    }
    return 0;
}

然后是结果：

\frac{1}{6}((2y'+3y'^2+y'^3+4x'+3x'y'-6x'^2-3x'^2y'+2x'^3)\\ +(2+3y'+y'^2-3x'-2x'y'+x'^2)\color{red}(b-2a)\\ \color{black}+(3y'+2y'^2+6x'+2x'y'-4x'^2)\color{red}(a+b)\\ \color{black}+(3+2y'-2x')\color{red}(a+b)(b-2a)\\ \color{black}+(y'+2x')\color{red}(a+b)^2\\ \color{black}+\color{red}(a+b)^2(b-2a)\color{black})

有没有可能项数更少？反正我试不出来，也不知道存不存在算法能高效解决这个问题。搞了半天才跟我三年前随手一拆的结果一样。

当 x\leq a,y+1\leq b 时：

w\sum_{i=x'-a}^d i(d+1-i)

化简后面一坨

\frac{1}{2}(x'-a+d)(d-x'+a+1)(d+1)-\frac{1}{6}d(d+1)(2d+1)+\frac{1}{6}(x'-a-1)(x'-a)(2x'-2a-1)

这次生成系数的代码：

#include<bits/stdc++.h>
using namespace std;

const int MAXN=100;
const double eps=1e-5;

int n;
double A[MAXN+5][MAXN+5];

int f(int a,int x,int d)
{
    int s1=(x-a+d)*(d-x+a+1)*(d+1)/2;
    int s2=d*(d+1)*(2*d+1)/6;
    int s3=(x-a-1)*(x-a)*(2*x-2*a-1)/6;
    return s1-s2+s3;
}

bool Solve()
{
    for(int i=1;i<=n;i++)
    {
        int minn=0;
        for(int j=i;j<=n;j++)
            if(fabs(A[j][i])>eps && (!minn || A[j][i]<A[minn][i]))
                minn=i;
        if(!minn) return 0;
        swap(A[minn],A[i]);
        for(int j=n+1;j>=i;j--) A[i][j]/=A[i][i];
        for(int j=i+1;j<=n;j++)
            for(int k=n+1;k>=i;k--)
                A[j][k]-=A[j][i]*A[i][k];
    }
    for(int i=n;i;i--)
        for(int j=i-1;j;j--)
            A[j][n+1]-=A[j][i]*A[i][n+1];
    return 1;
}

int main()
{
    //freopen("coe.txt","w",stdout);
    srand(time(0));
    while(1)
    {
        n=20;
        for(int T=1;T<=n;T++)
        {
            int cnt=0;
            int a=rand()%10,x=rand()%10,d=rand()%10;
            A[T][n+1]=f(a,x,d);
            for(int i=0;i<=3;i++)
                for(int j=0;j<=3;j++)
                    for(int k=0;k<=3;k++)
                        if(i+j+k<=3)
                        {
                            ++cnt;
                            A[T][cnt]=pow(a,i)*pow(x,j)*pow(d,k);
                        }
        }
        if(!Solve()) continue;
        int cnt=0;
        printf("\\frac{1}{6}(");
        for(int i=0;i<=3;i++)
        {
            printf("+(");
            bool OK=0;
            for(int j=0;j<=3;j++)
                for(int k=0;k<=3;k++)
                    if(i+j+k<=3)
                    {
                        ++cnt;
                        int v=round(A[cnt][n+1]*6);
                        if(!v) continue;
                        if(v>0 && OK) printf("+%d",v);
                        else printf("%d",v),OK=1;
                        if(j>1) printf("x\'^%d",j);
                        else if(j==1) printf("x\'");
                        if(k>1) printf("d^%d",k);
                        else if(k==1) printf("d");
                    }
            printf(")");
            if(i>1) printf("a^%d",i);
            else if(i==1) printf("a");
            printf("\\\\\n");
        }
        break;
    }
    return 0;
}

结果：

\frac{1}{6}((2d+3d^2+d^3+4x'+3x'd-6x'^2-3x'^2d+2x'^3)\\ +(-4-3d+12x'+6x'd-6x'^2)\color{red}a\color{black}\\ +(-6-3d+6x')\color{red}a^2\color{black}\\ -2\color{red}a^3\color{black})\\

最后是右上角加对查询的贡献，当 x\leq a,y\leq b 时（最左侧权值为 d）：

w\sum_{i=d-a+x}^d i(b-y+1)

化：

\frac{1}{2}(2d-a+x)(a-x+1)(b-y+1)

码：

#include<bits/stdc++.h>
using namespace std;

const int MAXN=100;
const double eps=1e-5;

int n;
double A[MAXN+5][MAXN+5];

int f(int a,int b,int x,int y,int d)
{return (2*d-a+x)*(a-x+1)*(b-y+1)/2;}

bool Solve()
{
    for(int i=1;i<=n;i++)
    {
        int minn=0;
        for(int j=i;j<=n;j++)
            if(fabs(A[j][i])>eps && (!minn || A[j][i]<A[minn][i]))
                minn=i;
        if(!minn) return 0;
        swap(A[minn],A[i]);
        for(int j=n+1;j>=i;j--) A[i][j]/=A[i][i];
        for(int j=i+1;j<=n;j++)
            for(int k=n+1;k>=i;k--)
                A[j][k]-=A[j][i]*A[i][k];
    }
    for(int i=n;i;i--)
        for(int j=i-1;j;j--)
            A[j][n+1]-=A[j][i]*A[i][n+1];
    return 1;
}

int main()
{
    //freopen("coe.txt","w",stdout);
    srand(time(0));
    while(1)
    {
        n=56;
        for(int T=1;T<=n;T++)
        {
            int cnt=0;
            int a=rand()%10,b=rand()%10;
            int x=rand()%10,y=rand()%10,d=rand()%10;
            A[T][n+1]=f(a,b,x,y,d);
            for(int i=0;i<=3;i++)
                for(int j=0;j<=3;j++)
                    for(int k=0;k<=3;k++)
                        for(int l=0;l<=3;l++)
                            for(int m=0;m<=3;m++)
                                if(i+j+k+l+m<=3)
                                {
                                    ++cnt;
                                    A[T][cnt]=pow(a,i)*pow(b,j)*pow(x,k)*pow(y,l)*pow(d,m);
                                }
        }
        if(!Solve()) continue;
        int cnt=0;
        printf("\\frac{1}{2}(");
        for(int i=0;i<=3;i++)
            for(int j=0;j<=3;j++)
            {
                if(i+j>3) continue;
                printf("+(");
                bool OK=0;
                for(int k=0;k<=3;k++)
                    for(int l=0;l<=3;l++)
                        for(int m=0;m<=3;m++)
                            if(i+j+k+l+m<=3)
                            {
                                ++cnt;
                                int v=round(A[cnt][n+1]*2);
                                if(!v) continue;
                                if(v>0 && OK) printf("+%d",v);
                                else printf("%d",v),OK=1;
                                if(k>1) printf("x^%d",k);
                                else if(k==1) printf("x");
                                if(l>1) printf("y^%d",l);
                                else if(l==1) printf("y");
                                if(m>1) printf("d^%d",m);
                                else if(m==1) printf("d");
                            }
                printf(")");
                if(i>1) printf("a^%d",i);
                else if(i==1) printf("a");
                if(j>1) printf("b^%d",j);
                else if(j==1) printf("b");
                printf("\\\\\n");
            }
        break;
    }
    return 0;
}

结果：

\frac{1}{2}((2d-2yd+x-2xd-xy+2xyd-x^2+x^2y)\\ +(2d+x-2xd-x^2)\color{red}b\color{black}\\ +(-1+2d+y-2yd+2x-2xy)\color{red}a\color{black}\\ +(-1+2d+2x)\color{red}ab\color{black}\\ +(-1+y)\color{red}a^2\color{black}\\ -\color{red}a^2b\color{black})\\

准备工作结束后开始写代码。麻烦点除了需要把这坨系数准确地打进程序里之外（或者可以另写个代码生成这段程序），还要获取八个方向所有三角的信息。一个个写坐标太折磨了，这里采用了旋转加翻折操作使得需要处理的三角都在修改中心的右方且角度与之前的分析一致，然后向 Solve 函数传入顶点坐标相对中心的偏移量，就能大大减少代码量。

时间 O(m \log^2 m) 空间 O(m)，至于常数，每个方向要做三次 4m 规模的 CDQ，需要维护的系数个数分别为 6,4,6，所以大概是 8\times 4\times (6+4+6)=512。感觉还是太愚昧了所以喜提最劣解，不如 FunnyCreatress 的差分做法。

关于答案对 2^{30} 取模，可以先 unsigned int 算出六倍答案模 2^{32} 的答案，然后去掉最高位再右移一位（写法上可以是先左移一位再右移两位）得到的是三倍答案模 2^{30} 的结果，此时 3 存在逆元，用 exgcd 算出来乘到答案里就好了。

参考代码：

#include<bits/stdc++.h>
using namespace std;

typedef unsigned int ui;

const int MAXN=1e5,MAXV=1e8,SIZE=4*MAXN;

int Read()
{
    int res=0;char c;
    while(!isdigit(c=getchar()));
    while(isdigit(c)) res=res*10+c-'0',c=getchar();
    return res;
}
void Print(ui x)
{
    if(x<10) {putchar('0'+x);return;}
    ui y=x/10; Print(y),putchar('0'+x-y*10);
}

int m,X1[MAXN+5],X2[MAXN+5],Y1[MAXN+5],Y2[MAXN+5];
bool opt[MAXN+5];
ui ans[MAXN+5];
int X[SIZE+5],Y[SIZE+5],ID[SIZE+5],Q[SIZE+5],q; ui val[SIZE+5][6];

void Reverse()//左右翻转 
{
    for(int i=1;i<=m;i++)
        if(opt[i]) X1[i]=MAXV-X1[i]+1,X2[i]=MAXV-X2[i]+1,swap(X1[i],X2[i]);
        else X1[i]=MAXV-X1[i]+1;
}
void Rotate()//顺时针旋转 
{
    for(int i=1;i<=m;i++)
        if(opt[i]) {int L=X1[i],R=X2[i],D=Y1[i],U=Y2[i]; X1[i]=D,X2[i]=U,Y1[i]=MAXV-R+1,Y2[i]=MAXV-L+1;}
        else swap(X1[i],X2[i]),X2[i]=MAXV-X2[i]+1;
}
int dsc[SIZE+5],dnum;
void Discrete()
{
    dnum=0;
    for(int i=1;i<=q;i++) dsc[++dnum]=Y[i];
    stable_sort(dsc+1,dsc+dnum+1),dnum=unique(dsc+1,dsc+dnum+1)-dsc-1;
    for(int i=1;i<=q;i++) Y[i]=lower_bound(dsc+1,dsc+dnum+1,Y[i])-dsc;
}

struct BIT
{
    ui sum[SIZE+5];
    void Plus(int x,ui v) {for(;x<=dnum;x+=x&-x) sum[x]+=v;}
    ui Ask(int x) {ui res=0; for(;x;x-=x&-x) res+=sum[x]; return res;}
}mapn[6];
void CDQ(int L,int R,int lim)
{
    if(L==R) return;
    int mid=(L+R)>>1;
    CDQ(L,mid,lim),CDQ(mid+1,R,lim);
    for(int i=mid+1,j=L;i<=R;i++)
        if(ID[Q[i]])
        {
            for(;j<=mid && X[Q[j]]<=X[Q[i]];j++)
                if(!ID[Q[j]]) for(int k=0;k<lim;k++) mapn[k].Plus(Y[Q[j]],val[Q[j]][k]);
            for(int k=0;k<lim;k++) ans[ID[Q[i]]]+=val[Q[i]][k]*mapn[k].Ask(Y[Q[i]]);
        }
    for(int i=mid+1,j=L;i<=R;i++)
        if(ID[Q[i]])
            for(;j<=mid && X[Q[j]]<=X[Q[i]];j++)
                if(!ID[Q[j]]) for(int k=0;k<lim;k++) mapn[k].Plus(Y[Q[j]],-val[Q[j]][k]);
    inplace_merge(Q+L,Q+mid+1,Q+R+1,[](int a,int b) {return X[a]<X[b];});
}

void Modify1(int x,int y,int d,ui w)
{
    ui x_=d+x,y_=d-y;
    ++q;
    X[q]=x+y,Y[q]=-y,ID[q]=0,Q[q]=q;
    val[q][0]=2*y_+3*y_*y_+y_*y_*y_+4*x_+3*x_*y_-6*x_*x_-3*x_*x_*y_+2*x_*x_*x_;
    val[q][1]=2+3*y_+y_*y_-3*x_-2*x_*y_+x_*x_;
    val[q][2]=3*y_+2*y_*y_+6*x_+2*x_*y_-4*x_*x_;
    val[q][3]=3+2*y_-2*x_;
    val[q][4]=y_+2*x_;
    val[q][5]=1;
    for(int i=0;i<6;i++) val[q][i]*=w;
}
void Query1(ui x,int y,int t,ui w)
{
    ++q;
    X[q]=x+y,Y[q]=-y,ID[q]=t,Q[q]=q;
    val[q][0]=1;
    val[q][1]=y-2*x;
    val[q][2]=x+y;
    val[q][3]=(x+y)*(y-2*x);
    val[q][4]=(x+y)*(x+y);
    val[q][5]=(x+y)*(x+y)*(y-2*x);
    for(int i=0;i<6;i++) val[q][i]*=w;
}
void Modify2(int x,int y,ui d,ui w)
{
    ui x_=d+x;
    ++q;
    X[q]=x,Y[q]=y+1,ID[q]=0,Q[q]=q;
    val[q][0]=2*d+3*d*d+d*d*d+4*x_+3*x_*d-6*x_*x_-3*x_*x_*d+2*x_*x_*x_;
    val[q][1]=-4-3*d+12*x_+6*x_*d-6*x_*x_;
    val[q][2]=-6-3*d+6*x_;
    val[q][3]=-2;
    for(int i=0;i<4;i++) val[q][i]*=w;
}
void Query2(ui x,int y,int t,ui w)
{
    ++q;
    X[q]=x,Y[q]=y,ID[q]=t,Q[q]=q;
    val[q][0]=1;
    val[q][1]=x;
    val[q][2]=x*x;
    val[q][3]=x*x*x;
    for(int i=0;i<4;i++) val[q][i]*=w;
}
void Modify3(ui x,ui y,ui d,ui w)
{
    ++q;
    X[q]=x,Y[q]=y,ID[q]=0,Q[q]=q;
    val[q][0]=2*d-2*y*d+x-2*x*d-x*y+2*x*y*d-x*x+x*x*y;
    val[q][1]=2*d+x-2*x*d-x*x;
    val[q][2]=-1+2*d+y-2*y*d+2*x-2*x*y;
    val[q][3]=-1+2*d+2*x;
    val[q][4]=-1+y;
    val[q][5]=-1;
    for(int i=0;i<6;i++) val[q][i]*=3*w;
}
void Query3(ui x,int y,int t,ui w)
{
    ++q;
    X[q]=x,Y[q]=y,ID[q]=t,Q[q]=q;
    val[q][0]=1;
    val[q][1]=y;
    val[q][2]=x;
    val[q][3]=x*y;
    val[q][4]=x*x;
    val[q][5]=x*x*y;
    for(int i=0;i<6;i++) val[q][i]*=w;
}
void Solve(int dx,int dy)
{
    q=0;
    for(int i=1;i<=m;i++)
        if(opt[i])
        {
            Query1(X2[i],Y2[i],i,1); 
            Query1(X1[i]-1,Y2[i],i,-1);
            Query1(X2[i],Y1[i]-1,i,-1);
            Query1(X1[i]-1,Y1[i]-1,i,1);
        }
        else if(dx<Y1[i])
        {
            Modify1(X1[i]+dx,X2[i]+dy,Y1[i]-dx,Y2[i]);
            Modify1(X1[i]+Y1[i],X2[i]-Y1[i]+dx+dy,0,-Y2[i]);
        }
    Discrete(),CDQ(1,q,6);
    q=0;
    for(int i=1;i<=m;i++)
        if(opt[i])
        {
            Query2(X2[i],Y2[i],i,1); 
            Query2(X1[i]-1,Y2[i],i,-1);
            Query2(X2[i],Y1[i]-1,i,-1);
            Query2(X1[i]-1,Y1[i]-1,i,1);
        }
        else if(dx<Y1[i])
        {
            Modify2(X1[i]+dx,X2[i]+dy,Y1[i]-dx,Y2[i]);
            Modify2(X1[i]+Y1[i],X2[i]-Y1[i]+dx+dy,0,-Y2[i]);
        }
    Discrete(),CDQ(1,q,4);
    q=0;
    for(int i=1;i<=m;i++)
        if(opt[i])
        {
            Query3(X2[i],Y2[i],i,1); 
            Query3(X1[i]-1,Y2[i],i,-1);
            Query3(X2[i],Y1[i]-1,i,-1);
            Query3(X1[i]-1,Y1[i]-1,i,1);
        }
        else if(dx<Y1[i])
        {
            Modify3(X1[i]+Y1[i],X2[i]-Y1[i]+dx+dy+1,0,-Y2[i]);
            Modify3(X1[i]+Y1[i],X2[i]+dy+1,0,Y2[i]);
        }
    Discrete(),CDQ(1,q,6);
}

int main()
{
    m=Read();
    for(int i=1;i<=m;i++) opt[i]=Read()-1,X1[i]=Read(),X2[i]=Read(),Y1[i]=Read(),Y2[i]=Read();
    Solve(1,0);
    Rotate();
    Solve(2,-1);
    Rotate();
    Solve(1,-1);
    Rotate();
    Solve(1,0);
    Reverse();
    Solve(1,0);
    Rotate();
    Solve(1,-1);
    Rotate();
    Solve(1,-1);
    Rotate();
    Solve(0,0);
    for(int i=1;i<=m;i++) if(opt[i]) Print((((ans[i]<<1)>>2)*715827883u)&1073741823u),putchar('\n');
    return 0;
}