夏令营501-511NOIP训练17——蛇形矩阵

淩亂°似流年 2023-06-05 12:44 8阅读 0赞

传送门：[QAQQAQ][]

题意：话说小X在孩提时，都会做标准的蛇形矩阵了，发现很好玩。现在的小X很想对其进行改版，变为如下类型的一个无限大蛇形数阵：  
![1.jpg][]  
令S(x)表示以1为左上角，x为右下角的矩形内所有数之和。例如S(12)就是具有深色背景的数之和。  
给定n，对于“以1为左上角，n为右下角的矩形”内的每一个数i，计算所有S(i)之和。

思路：神仙数学题。。。

先看暴力：我们通过人工YY和简单推算可以得出一个公式：$ans=\\sum\_\{i=1\}^\{n\}\\sum\_\{j=1\}^\{m\}a\_\{i,j\}\*\\left ( n-i+1 \\right )\*\\left ( m-j+1 \\right )$

那么我们可以在$O(n)$的复杂度内解决问题

现在开始优化环节：我们可以把一个给定矩形拆分成若干条链，使每一条链都是一个连续的数列（莫名想到树链剖分。。。）

具体拆分方法为：若是一个长方形，则先把多余的行或列拆出来，把它变成一个正方形。然后对于每一圈，我们先竖着拆$rnd$个，再横着拆$rnd-1$个，拆完为止。

现在我们的任务是求一个连续序列的和。每一个的统计总和为$a\_\{i,j\}\*\\left ( n-i+1 \\right )\*\\left ( m-j+1 \\right )$，而且这个序列一定在同一行或同一列，所以我们可以把$\\left ( n-i+1 \\right )$或$\\left ( m-j+1 \\right )$中的一项提取出来，这样剩下的就是两项相乘求和了。而我们发现这两项分别在两个公差为$1$或$-1$的等差数列中。

现在开始数学环节——暴力推公式

令$f(s1,s2,d1,d2,n)$为首项为$s1$，公差为$d1$的等差数列和首项为$s2$，公差为$d2$的等差数列每一项分别相乘，前$n$项乘积的和

即：$f(s1,s2,d1,d2,n)$=$\\sum\_\{i=1\}^\{n\}\\left \[ s1+\\left ( i-1 \\right )\*d1 \\right \]\*\\left \[ s2+\\left ( i-1 \\right )\*d2 \\right \]$

=$\\sum\_\{i=1\}^\{n\}(s1\*s2)+\\sum\_\{i=1\}^\{n\}(s1\*d2\*(i-1))+\\sum\_\{i=1\}^\{n\}( d1\*s2\*(i-1))+\\sum\_\{i=1\}^\{n\}( d1\*d2\*(i-1)^\{2\} )$

=$n\*s1\*s2+s1\*d2\*\\sum\_\{i=1\}^\{n-1\}i+s2\*d1\*\\sum\_\{i=1\}^\{n-1\}i+d1\*d2\*\\sum\_\{i=1\}^\{n-1\}i^\{2\}$

=$n\*s1\*s2+s1\*d2\*\\frac\{(n-1)\*n\}\{2\}+d1\*s2\*\\frac\{(n-1)\*n\}\{2\}+d1\*d2\*\\frac\{(n-1)\*n\*(2\*n-1)\}\{6\}$

那么我们现在就可以用$O(1)$的时间算出每一条链，而一共拆成$\\sqrt\{n\}$条链，所以总复杂度$\\sqrt\{n\}$

代码：（注意代码实现时$x,y$很容易弄混）

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

#include<bits/stdc++.h>
    #define pii pair<ll,ll>
    #define mk make_pair
    using namespace std;
    typedef long long ll;
    const ll MOD=(ll)1e9+7;
    const ll inv2=500000004;
    const ll inv6=166666668;
     
    pii find(ll x)
    {
        ll rnd=(ll)sqrt(x);
        if(rnd*rnd!=x) rnd++;
        ll End=rnd*rnd;
        ll dif=End-x;
        ll posx,posy;
        if(rnd%2==0)
        {
            posx=rnd; posy=1;
            if(dif<=rnd-1) posy+=dif;
            else
            {
                posy=rnd; dif-=(rnd-1);
                posx-=dif;
            }
        }
        else
        {
            posx=1,posy=rnd;
            if(dif<=rnd-1) posx+=dif;
            else
            {
                posx=rnd; dif-=(rnd-1);
                posy-=dif;
            }
        }
        return mk(posx,posy);
    }
     
    void fn(ll &x)
    {
        x=(x%MOD+MOD)%MOD;
    }
     
    ll f(ll s1,ll s2,ll d1,ll d2,ll n)
    {
        if(n==0) return 0;
        ll ret=n*s1%MOD*s2%MOD; fn(ret);
        ret=(ret+n*(n-1)%MOD*inv2%MOD*(s1*d2%MOD+s2*d1%MOD))%MOD; fn(ret);
        ret=(ret+n*(n-1)%MOD*(2*n-1)%MOD*inv6%MOD*d1%MOD*d2%MOD); fn(ret);
        return ret;
    }
     
    int main()
    {
        ll x,ans=0;
        scanf("%lld",&x);
        pii pos=find(x);
        ll n=pos.first,m=pos.second;
        ll nowx=n,nowy=m;
        while(nowx>m)
        {
            if(nowx%2==1)
            {
                ll k=n-nowx+1;
                ll s1=(nowx-1)*(nowx-1)+1; ll d1=1;
                ll s2=m; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,m))%MOD;
            }
            else
            {
                ll k=n-nowx+1;
                ll s1=nowx*nowx; ll d1=-1;
                ll s2=m; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,m))%MOD;
            }
            nowx--;
        }
        while (nowy>n)
        {
            if(nowy%2==1)
            {
                ll k=m-nowy+1;
                ll s1=nowy*nowy; ll d1=-1;
                ll s2=n; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,n))%MOD;
            }
            else
            {
                ll k=m-nowy+1;
                ll s1=(nowy-1)*(nowy-1)+1; ll d1=1;
                ll s2=n; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,n))%MOD;
            }
            nowy--;
        }
        while(nowx>0&&nowy>0)
        {
            if(nowx%2==1)
            {
                ll k=m-nowy+1;
                ll s1=nowy*nowy; ll d1=-1;
                ll s2=n; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,nowx))%MOD;
                k=n-nowx+1;
                s1=(nowx-1)*(nowx-1)+1; d1=1;
                s2=m; d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,nowy-1))%MOD;
            }
            else
            {
                ll k=m-nowy+1;
                ll s1=(nowy-1)*(nowy-1)+1; ll d1=1;
                ll s2=n; ll d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,nowx))%MOD;
                k=n-nowx+1;
                s1=nowx*nowx; d1=-1;
                s2=m; d2=-1;
                ans=(ans+k*f(s1,s2,d1,d2,nowy-1))%MOD;
            }
            nowx--; nowy--;
        }
        cout<<ans<<endl;
        return 0;
    }

转载于:https://www.cnblogs.com/Forever-666/p/11286038.html

[QAQQAQ]: https://dev.xjoi.net/contest/1216/problem/2
[1.jpg]: https://dev.xjoi.net/images/problem/1734/1.jpg
[ContractedBlock.gif]: https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif
[ExpandedBlockStart.gif]: /images/20230601/fed27afa89ab4a848dfc9a37ac7994d0.png