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

传送门：QAQQAQ

题意：话说小X在孩提时，都会做标准的蛇形矩阵了，发现很好玩。现在的小X很想对其进行改版，变为如下类型的一个无限大蛇形数阵：

令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 ight )*left ( m-j+1 ight )$

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

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

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

现在我们的任务是求一个连续序列的和。每一个的统计总和为$a_{i,j}*left ( n-i+1 ight )*left ( m-j+1 ight )$，而且这个序列一定在同一行或同一列，所以我们可以把$left ( n-i+1 ight )$或$left ( m-j+1 ight )$中的一项提取出来，这样剩下的就是两项相乘求和了。而我们发现这两项分别在两个公差为$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 ight )*d1 ight ]*left [ s2+left ( i-1 ight )*d2 ight ]$

　　　　　　　　　　=$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$很容易弄混）

#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;
}