题目翻译
题解
首先很容易想到,我们只需要算“二维前缀和”,即以某点为右下角的子矩阵内的元素和。然后要求的那个值用容斥原理减一减就能算出来。
下面看怎么算二维前缀和(当然是模意义下的),记为S(i,j):
以6*6矩阵为例:
假设我们想算S(2,5)。它由两部分组成:第一部分是左上角的2*2矩阵,也就是1+2+3+4,这是一个自然数列求和。第二部分是三列:5,6;10,11;和17,18。特点是什么呢?我们来看每一列的和,也就是5+6,10+11,17+18这三项。第一项5+6是等差数列求和。然后相邻两项之差,第一次是(10+11)-(5+6)=2*5,第二次是(17+18)-(10+11)=2*(5+2),如果继续算的话,还有(26+27)-(17+18)=2*(5+2+2)……也就是说,相邻两项的差每次增加2*2(第一个2是x,因为x行的规律相同)。
可以发现,我们只需要实现两个函数:①等差数列求和(一阶级数求和),②二阶级数求和。二阶级数的意义就是说,其差分序列是一个等差数列(相邻元素的差匀速增大)。
先不考虑我们怎么实现一阶/二阶级数求和。有一个现实的问题:题目要求输出后十位,那么我们可以对10^10取模。但怎么判断答案是否小于10^10呢?
方法很tricky:先对10^10取模算一下答案,再对10^10+9取模算一下答案。如果答案相同,就认为答案小于10^10.
那么,具体怎么实现两个级数求和呢?
我先是看着数据组数T<=10^5不大,所以写了个二分快速乘+矩阵乘法(二分快速乘是因为10^10的平方会爆long long),然后光荣TLE……
然后想着既然矩阵乘法会T那就推公式。等差数列的公式谁都会吧……二阶级数的公式可以展开后
利用公式:1^2 + 2^2 + ... + n^2 = n(n+1)(2n+1) / 6.
注意,不能用逆元做除法(10^10非质数)。除法怎么办呢?二阶级数的结果中有类似(n-2)(n-1)n/6这样的式子,显然n-2,n-1,n三项都不会爆long long,那我们把这三个数写出来,看谁是2或3的倍数,直接除掉,最后乘起来……许多乘法都需要用二分快速乘,尤其是等差数列中项数最大可能是(10^9)^2=10^18,不注意这一点会WA。
然后光荣地又T了。没错,对于T=10^5,O(TlogN)的算法TLE了……
然后把二分快速乘改为常数快速乘:
return ((x*y-(LL)(((long double)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;
就过了。
(doge脸)人类为何要互相伤害……卡常数的意义是什么……
代码
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define sqr(x) ((x)*(x))
typedef long long LL;
LL MOD;
inline LL realmod(LL x){
x%=MOD;
if(x<0) x+=MOD;
return x;
}
inline LL quickmul(LL x,LL y){
x=realmod(x),y=realmod(y);
return ((x*y-(LL)(((long double)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;
}
inline LL series_1(LL a,LL b,LL n){//一阶级数
//首项为a,公差为b,共n项
a=realmod(a),b=realmod(b);
LL ans=quickmul(a,n);
static LL c[2];
c[0]=n-1,c[1]=n;
for(int i=0;i<2;i++){if(c[i]%2==0){c[i]/=2;break;}}
ans+=quickmul(b,quickmul(c[0],c[1]));
return realmod(ans);
}
inline LL series_2(LL a2,LL a1,LL b,LL n){//二阶级数
//首项为a2,第一次加a1,然后每次多加b,共n项
a2=realmod(a2),a1=realmod(a1),b=realmod(b);
LL ans=series_1(a2,a1,n);
static LL c[3];
c[0]=n-2,c[1]=n-1,c[2]=n;
for(int i=0;i<3;i++){if(c[i]%2==0){c[i]/=2;break;}}
for(int i=0;i<3;i++){if(c[i]%3==0){c[i]/=3;break;}}
ans+=quickmul(b,quickmul(c[0],quickmul(c[1],c[2])));
ans=realmod(ans);
return realmod(ans);
}
inline LL calc(LL x,LL y){
if(x<1||y<1) return 0;
LL ans;
ans=series_1(1,1,sqr(min(x,y)));
if(x>y){
LL a2=series_1((y+1)*(y+1),-1,y);
LL a1=realmod((2*y+3)*y);
LL b=realmod(2*y);
ans+=series_2(a2,a1,b,x-y);
}
else if(y>x){
LL a2=series_1(x*x+1,1,x);
LL a1=realmod((2*x+1)*x);
LL b=realmod(2*x);
ans+=series_2(a2,a1,b,y-x);
}
return ans;
}
inline LL calc(LL x1,LL y1,LL x2,LL y2){
return realmod(calc(x2,y2)-calc(x2,y1-1)-calc(x1-1,y2)+calc(x1-1,y1-1));
}
LL mod1=1e10,mod2=1e10+9;
inline void work(void){
int x1,y1,x2,y2;
scanf("%d%d%d%d",&x1,&y1,&x2,&y2);
MOD=mod1;
LL res1=calc(x1,y1,x2,y2);
MOD=mod2;
LL res2=calc(x1,y1,x2,y2);
if(res1!=res2){
printf("...");
printf("%.10lld\n",res1);
}
else printf("%lld\n",res1);
}
int main(){
freopen("endlessmatrix.in","r",stdin);
freopen("endlessmatrix.out","w",stdout);
int T;
scanf("%d",&T);
while(T--) work();
return 0;
}