【数学】[Baltic2016][BZOJ5184] Spiral

分析：

恶心至极。。。
首先，把图拆成一个个的环，
第一层是：1
第二层是：2，3，4，5，6，7，8，9
……
然后，就可以发现，单独的一个环的贡献，可以在O(1)复杂度内算出来。（拆成四条边）

然而这要T

问题就在于如何快速求多个环的贡献。

首先，我们可以把原矩阵也拆分了，拆分成数个部分环的组合：

这样一来，每个部分环中，被剖分的矩形只有两种可能：


可以通过推倒证明，在这两种情况下，每一个环的贡献一定满足：
$Si=A?i^3+B?i^2+C?i^1+D?i^0$
（其中i表示是从上一层开始的第i个环）

所以，要求K个环的贡献和，就只需要算：
$A?(0^3+1^3+2^3+\dots \dots +（k?1{)}^3)+B?(0^2+1^2+\dots \dots +(j?1{)}^2)+C?(0^1+1^1+\dots \dots +（k?1{)}^1)+D?(0^0+1^0+2^0+\dots \dots +(k?1{)}^0)$
这几个都是有公式的，自己去套就好了。

然后，从证明过程中可以看出，推出A,B,C,D是几乎不可能的。但有一个更好的方法：我们可以暴力算出前4个环的贡献，然后将A,B,C,D反解出来。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#define SF scanf
#define PF printf
#define MAXN 100010
#define MOD 1000000007
#define y0 warspitewarspitewarspitewarspitewarspitewarspite
using namespace std;
typedef long long ll;
vector<ll> st;
ll ans;
ll inv6=(MOD+1ll)/6ll,inv2=(MOD+1ll)/2ll,v[5];
ll pres(ll x){
    return (1ll+x)*x/2ll%MOD;
}
ll calc(ll l1,ll r1,ll l2,ll r2,ll val){
    ll l=max(l1,l2);ll r=min(r1,r2);if(l<=r)return (pres(r)-pres(l-1)+(r-l+1)*(val%MOD)%MOD)%MOD;return 0;
}
ll count(ll L,ll l,ll r,ll d,ll u){
    if(L==0)return 	l<=0&&r>=0&&u>=0&&d<=0;ll val=(2ll*L+1ll)*(2ll*L+1ll);ll res=0;val-=L;if(u>=-L&&d<=-L)res=res+calc(l,r,-L,L,val);val-=2ll*L;if(l<=-L&&r>=-L)res=res+calc(d,u,-L+1,L-1,val);val-=2ll*L;if(u>=L&&d<=L)res=res+calc(l,r,-L,L,val);val-=2ll*L;if(l<=L&&r>=L)res=res+calc(d,u,-L+1,L-1,val);res=(res%MOD+MOD)%MOD;return res;
}	
int main(){
    
// freopen("t.in.4","r",stdin);
// freopen("wa.out","w",stdout);int n,q;SF("%d%d",&n,&q);for(int i=1;i<=q;i++){
    ans=0;ll x0,x1,y0,y1;SF("%lld%lld%lld%lld",&x0,&y0,&x1,&y1);st.clear();st.push_back(0);st.push_back(1);st.push_back(2);st.push_back(abs(x0)-1);st.push_back(abs(x0));st.push_back(abs(x0)+1);st.push_back(abs(x1)-1);st.push_back(abs(x1));st.push_back(abs(x1)+1);st.push_back(abs(y0)-1);st.push_back(abs(y0));st.push_back(abs(y0)+1);st.push_back(abs(y1)-1);st.push_back(abs(y1));st.push_back(abs(y1)+1);sort(st.begin(),st.end());for(int i=1;i<int(st.size());i++){
    int las=st[i-1];int now=st[i];if(las<0||las==now)continue;int k=now-las;if(k<4){
    for(int i=las;i<now;i++)(ans+=count(i,x0,x1,y0,y1))%=MOD;}else{
    for(int i=0;i<4;i++)v[i]=count(i+las,x0,x1,y0,y1);ll A=(((v[3]-2ll*v[2]+v[1])-(v[2]-2ll*v[1]+v[0]))%MOD+MOD)%MOD;A=A*inv6%MOD;v[1]=(v[1]-A+MOD)%MOD;v[2]=(v[2]-8ll*A%MOD+MOD)%MOD;v[3]=(v[3]-27ll*A%MOD+MOD)%MOD;ll B=((v[2]-2ll*v[1]+v[0])%MOD+MOD)%MOD;B=B*inv2%MOD;v[1]=(v[1]-B+MOD)%MOD;v[2]=(v[2]-4ll*B%MOD+MOD)%MOD;v[3]=(v[3]-9ll*B%MOD+MOD)%MOD;ll C=v[1]-v[0];ll D=v[0];ll K=k;ll SUM=K*(K-1ll)%MOD*inv2%MOD;
// PF("{%lld %lld %lld %lld(%lld,%lld,%lld,%lld) [%d %d]}\n",A,B,C,D,v[0],v[1],v[2],v[3],las,now);ans=(ans+K*D%MOD)%MOD;ans=(ans+SUM*C%MOD)%MOD;ans=(ans+K*(K-1ll)%MOD*(2ll*k-1ll)%MOD*inv6%MOD*B%MOD)%MOD;ans=(ans+SUM*SUM%MOD*A%MOD)%MOD;}}ans=(ans+MOD)%MOD;PF("%lld\n",ans);}
}