Project Euler 453 Lattice Quadrilaterals 困难的计数问题
这是一道很综合的计数问题,对于思维的全面性,解法的过渡性,代码能力,细节处理,计数问题中的各种算法,像gcd、容斥、类欧几里德算法都有考察.
在省选模拟赛中做到了这题,然而数据范围是n,m小于等于1000.
首先有一个O(n^4m^4)的暴力.
然后开始计数,思路是:答案等于任取4个点的方案数+2*取4个点不为凸的方案.
前一部分相对来说容易统计,先用组合数算所有的,再把存在3点、4点共线的矩形的贡献减掉就好了.
这里用到了矩形框的思路,利用了容斥,而且在计数的时候用gcd作为工具,这个思路下面还会用到,这一部分的具体实现见代码.
后一部分,对我来说,我觉得是十分困难的,我经历了从O(n^2m^2)到O(n^2m),再到O(nmlog)的历程.
先从最基本的计数方向考虑,我们怎么统计取4个点不为凸的方案.
答案就是找到所有三角形,算出其内部点的和.
那么我们找所有三角形,继续沿用上面矩形框的思路,我们枚举最小的能把三角形框起来的矩形框.
然后我分了几种情况:
(以下所说的占有的顶点均指矩形的顶点,所说的占有顶点的三角形均满足其所有顶点都在矩形框上,设矩形的某一对平行边为长,另一对平行边为宽)
I.占有三个顶点的三角形
II.占有一个顶点的三角形
III.四种占有两个顶点的三角形
1.一边为对角线,另一顶点在宽上的三角形
2.一边为对角线,另一顶点在高上的三角形
3.以宽为底,另一顶点在对面边上的三角形
4.以长为底,另一顶点在对面边上的三角形
IV.一边为对角线一点在矩形内部的三角形(我用一晚上查错,才发现我漏掉了这种情况)
在计算其内部点的个数的时候,需要用到pick定理.
对于以上情况的计数,我采用的方法都是,先固定一个或两个顶点,算其贡献,然后再算由于对称性而产生的贡献.
这样我就想到了一种O(n^2m^2)的做法:
先O(nm)枚举矩形框,然后再O(1)计算I,O(nm)计算II(假定其一顶点为(0,0),另外两顶点在与(0,0)相对的边上滑动),O(n+m)计算III里的四种(与计算II的思路相似),O(nm)计算IV(枚举内部点作为顶点).
面积直接计算,沿边点数用gcd计算.
具体实现见代码:
#include <cstdio> #include <cstring> #include <algorithm> typedef long long LL; const int P=1000000007; const int N=1010; int gcd[N][N]; inline int GCD(int x,int y){ return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==0?y:GCD(y%x,x)); } inline void Init(){ register int i,j; for(i=0;i<=1000;++i) for(j=0;j<=1000;++j) gcd[i][j]=GCD(i,j); } inline int count_all(int n,int m){ int cnt=(n+1)*(m+1); int ret=41666667LL*cnt%P*(cnt-1)%P*(cnt-2)%P*(cnt-3)%P; register int i,j,tmp,sum; for(i=0;i<=n;++i) for(j=!i;j<=m;++j){ tmp=gcd[i][j]-1; if(!tmp)continue; sum=(n-i+1)*(m-j+1)*(i&&j?2:1); ret=(ret-(LL)tmp*sum*(cnt-3)+P)%P; if(tmp==1)continue; ret=(ret+(LL)tmp*(tmp-1)/2%P*sum%P*3%P)%P; } return ret; } inline int count_rest(int n,int m){ int ret=0,sum,s,tmp; register int i,j,x,y,temp; for(i=1;i<=n;++i) for(j=1;j<=m;++j){ sum=(n-i+1)*(m-j+1); s=i*j+2; tmp=(s-j-i-gcd[i][j])*2; for(x=1;x<i;++x) for(y=1;y<j;++y){ temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]); tmp=(tmp+temp*2)%P; } x=i; for(y=1;y<j;++y){ temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]); tmp=(tmp+temp*2)%P; } y=j; for(x=1;x<i;++x){ temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]); tmp=(tmp+temp*2)%P; } x=0; for(y=1;y<j;++y){ temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]); tmp=(tmp+temp)%P; } y=0; for(x=1;x<i;++x){ temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]); tmp=(tmp+temp)%P; } for(x=1;x<i;++x) for(y=1;y<j;++y){ if(y*i-x*j<=0)continue; temp=y*i-x*j+2-(gcd[x][y]+gcd[i-x][j-y]+gcd[i][j]); tmp=(tmp+temp*2)%P; } ret=(ret+(LL)tmp*sum)%P; } return ret; } inline int calc(int n,int m){ if((n+1)*(m+1)<4)return 0; return (count_all(n,m)+2LL*count_rest(n,m))%P; } int main(){ Init(); int n,m; scanf("%d%d",&n,&m); printf("%d\n",calc(n,m)); return 0; }
调过之后,我观察我的程序,我想到了把我的程序优化到O(n^2m)的做法:
仍然先O(nm)枚举矩形框,仍然O(1)计算I,对于II、III的情况,我们在一开始先预处理gcd二维前缀和,用矩形求和O(1)解决沿边点数,面积的话也可以O(1)计算,但是在计算IV的时候,我们虽然可以沿用刚刚的思路,但是我们还是要先O(n)枚举一维坐标,另一维O(1)解决.
具体实现见代码:
#include <cstdio> #include <cstring> #include <algorithm> typedef long long LL; const int P=1000000007; const int N=1000; int gcd[N+10][N+10],gcd_s[N+10][N+10],gcd_t[N+10][N+10],f[N+10][N+10]; inline int GCD(int x,int y){ return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==0?y:GCD(y%x,x)); } inline void Init(){ register int i,j; for(i=0;i<=N;++i) for(j=0;j<=N;++j) gcd_s[i][j]=gcd[i][j]=GCD(i,j); for(i=0;i<=N;++i) for(j=1;j<=N;++j) gcd_s[i][j]+=gcd_s[i][j-1]; for(i=0;i<=N;++i) for(j=0;j<=N;++j) gcd_t[i][j]=gcd_s[i][j]; for(i=0;i<=N;++i) for(j=1;j<=N;++j) gcd_s[j][i]+=gcd_s[j-1][i]; } inline int count_all(int n,int m){ int cnt=(n+1)*(m+1); int ret=41666667LL*cnt%P*(cnt-1)%P*(cnt-2)%P*(cnt-3)%P; register int i,j,tmp,sum; for(i=0;i<=n;++i) for(j=!i;j<=m;++j){ tmp=gcd[i][j]-1; if(!tmp)continue; sum=(n-i+1)*(m-j+1)*(i&&j?2:1); ret=(ret-(LL)tmp*sum*(cnt-3)+P)%P; if(tmp==1)continue; ret=(ret+(LL)tmp*(tmp-1)/2%P*sum%P*3%P)%P; } ret=(ret+P)%P; return ret; } inline int query(int a,int b,int x,int y){ if(a>x)return 0; if(b>y)return 0; int ret=gcd_s[x][y]; --a,--b; if(a>=0)ret-=gcd_s[a][y]; if(b>=0)ret-=gcd_s[x][b]; if(a>=0&&b>=0)ret+=gcd_s[a][b]; return ret; } #define query_(a,b,c) (gcd_t[a][c]-gcd_t[a][b-1]) inline int count_rest(int n,int m){ if(n>m)std::swap(n,m); int ret=0,sum,s; register LL tmp; register int i,j,x,y; for(i=1;i<=n;++i) for(j=1;j<=m;++j){ sum=(n-i+1)*(m-j+1); if(i>j){ ret=(ret+(LL)f[i][j]*sum%P)%P; continue; } s=i*j+2; tmp=s-j-i-gcd[i][j]; tmp+=s*(i-1LL)*(j-1LL)-(i)*(i-1)/2LL*(j)*(j-1)/2LL%P; tmp-=query(1,j,i-1,j)*(j-1LL)%P+query(i,1,i,j-1)*(i-1LL)%P+query(1,1,i-1,j-1); tmp+=s*(j-1)-(i)*(j)*(j-1)/2; tmp-=gcd[i][j]*(j-1LL)+query(i,1,i,j-1)+query(0,1,0,j-1); tmp+=s*(i-1)-(j)*(i)*(i-1)/2; tmp-=gcd[i][j]*(i-1LL)+query(1,j,i-1,j)+query(1,0,i-1,0); for(x=1;x<i;++x){ y=x*j/i+1; if(y>=j)break; tmp+=((j-1+y)*(j-y)*i>>1)-(x*j-2+gcd[i][j])*(j-y)-query_(x,y,j-1)-query_(i-x,1,j-y); } tmp*=2; tmp+=s*(j-1); tmp-=gcd[0][j]*(j-1LL)+query(i,1,i,j-1)*2LL; tmp+=s*(i-1); tmp-=gcd[i][0]*(i-1LL)+query(1,j,i-1,j)*2LL; tmp=(tmp%P+P)%P; f[j][i]=f[i][j]=tmp; ret=(ret+tmp*sum)%P; } return ret; } inline int calc(int n,int m){ if((n+1)*(m+1)<4)return 0; return (count_all(n,m)+2LL*count_rest(n,m))%P; } int main(){ Init(); int n,m; scanf("%d%d",&n,&m); printf("%d\n",calc(n,m)); return 0; }
这样仍然是不够优美的,因为,即使加入了优化,在n、m均为1000的数据下,仍然会跑两三秒,所以我们瞄准我们的瓶颈IV,对他进行优化:
I、II、III的计算方法同上,对于IV,我们只要可以在O(1)或者O(log)的时间复杂度下解决就可以了,我们发现,他的沿边点数实际上是,对角线的贡献加上我们枚举的矩形框内部点的gcd之和,再减去枚举点在对角线上的贡献,这个可以用矩形求和O(1)解决,而面积呢,我们可以推出叉积的式子,然后利用类欧几里德算法解决.
具体实现见代码:
#include <cstdio> #include <cstring> #include <algorithm> typedef long long LL; const int N=1000; const int P=1000000007; int gcd[N+10][N+10],gcd_s[N+10][N+10],f[N+10][N+10]; bool vis[N+10][N+10]; inline int GCD(int x,int y){ return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==0?y:GCD(y%x,x)); } inline void Init(){ register int i,j; for(i=0;i<=N;++i) for(j=0;j<=N;++j) gcd_s[i][j]=gcd[i][j]=GCD(i,j); for(i=0;i<=N;++i) for(j=1;j<=N;++j) gcd_s[i][j]+=gcd_s[i][j-1]; for(i=0;i<=N;++i) for(j=1;j<=N;++j) gcd_s[j][i]+=gcd_s[j-1][i]; } inline int count_all(int n,int m){ //计算所有可能的四个点 int cnt=(n+1)*(m+1); int ret=41666667LL*cnt%P*(cnt-1)%P*(cnt-2)%P*(cnt-3)%P; register int i,j,tmp,sum; for(i=0;i<=n;++i) for(j=!i;j<=m;++j){ tmp=gcd[i][j]-1; if(!tmp)continue; sum=(n-i+1)*(m-j+1)*(i&&j?2:1); ret=(ret-(LL)tmp*sum*(cnt-3)+P)%P; if(tmp==1)continue; ret=(ret+(tmp*(tmp-1LL)>>1)*sum*3)%P; } ret=(ret+P)%P; return ret; } inline int query(int a,int b,int x,int y){ //查询一个矩形内部的gcd和 if(a>x||b>y)return 0; int ret=gcd_s[x][y]; --a,--b; if(a>=0)ret-=gcd_s[a][y]; if(b>=0)ret-=gcd_s[x][b]; if(a>=0&&b>=0)ret+=gcd_s[a][b]; return ret; } struct Data{ LL f,g,h; }; inline Data lo(LL a,LL b,LL c,LL n){ //类欧几里德算法 Data ret; if(!a){ ret.f=ret.g=ret.h=0; return ret; } if(a>=c||b>=c){ ret=lo(a%c,b%c,c,n); ret.h+=n*(n+1)*(2*n+1)/6*(a/c)*(a/c)+(n+1)*(b/c)*(b/c)+2*(b/c)*ret.f+2*(a/c)*ret.g+(a/c)*(b/c)*n*(n+1); ret.g+=n*(n+1)*(2*n+1)/6*(a/c)+n*(n+1)/2*(b/c); ret.f+=n*(n+1)/2*(a/c)+(n+1)*(b/c); return ret; } LL m=(a*n+b)/c; Data tmp=lo(c,c-b-1,a,m-1); ret.f=n*m-tmp.f; ret.g=(m*n*(n+1)-tmp.f-tmp.h)/2; ret.h=m*(m+1)*n-2*tmp.g-2*tmp.f-ret.f; return ret; } inline int count_rest(int n,int m){ //计算内凹四边形做出的额外贡献 int ret=0,sum,s,cnt; Data get; register LL tmp; register int i,j; for(i=1;i<=n;++i) for(j=1;j<=m;++j){ sum=(n-i+1)*(m-j+1); //这个框出现的次数 if(vis[i][j]){ ret=(ret+(LL)f[i][j]*sum%P)%P; continue; //用于剪枝 } s=i*j+2; //初始化pick的部分内容 //以下所说的顶点均指矩形的顶点,所说的占有顶点的三角形均满足其所有顶点都在矩形框上 tmp=s-j-i-gcd[i][j]; //计算占有三个顶点的三角形的贡献 tmp+=s*(i-1LL)*(j-1)-(i*(i-1LL)*j*(j-1)>>2); //计算占有一个顶点的三角形的面积和 tmp-=query(1,j,i-1,j)*(j-1LL)+query(i,1,i,j-1)*(i-1LL)+query(1,1,i-1,j-1); //计算占有一个顶点的三角形的沿边点数和 tmp+=s*(j-1)-(i*j*(j-1)>>1); //计算第一类占有两个顶点的三角形的面积和 tmp-=gcd[i][j]*(j-1)+query(i,1,i,j-1)+query(0,1,0,j-1); //计算第一类占有两个顶点的三角形的沿边点数和 tmp+=s*(i-1)-(j*i*(i-1)>>1); //计算第二类占有两个顶点的三角形的面积和 tmp-=gcd[i][j]*(i-1)+query(1,j,i-1,j)+query(1,0,i-1,0); //计算第二类占有两个顶点的三角形的沿边点数和 cnt=((i-1)*(j-1)-(gcd[i][j]-1))>>1; //计算以对角线为分界线把矩形分为两部分后,其中一部分内部的点数 tmp-=(gcd[i][j]-2)*cnt+query(1,1,i-1,j-1)-(gcd[i][j]*(gcd[i][j]-1)>>1); //计算一边为对角线一点在矩形内部的三角形的沿边点数和 tmp*=2; //将由对称产生的贡献算入 get=lo(j,0,i,i-1); tmp+=2*j*get.g-i*get.h-i*get.f; //计算一边为对角线一点在矩形内部的三角形的面积和 tmp+=s*(j-1); //计算第三类占有两个顶点的三角形的面积和 tmp-=j*(j-1)+(query(i,1,i,j-1)<<1); //计算第三类占有两个顶点的三角形的沿边点数和 tmp+=s*(i-1); //计算第四类占有两个顶点的三角形的面积和 tmp-=i*(i-1)+(query(1,j,i-1,j)<<1); //计算第四类占有两个顶点的三角形的沿边点数和 tmp=(tmp%P+P)%P; vis[i][j]=vis[j][i]=true; f[i][j]=f[j][i]=tmp; //用于剪枝 ret=(ret+tmp*sum)%P; //将本次贡献算入返回值 } return ret; } inline int calc(int n,int m){ if((n+1)*(m+1)<4)return 0; return (count_all(n,m)+2LL*count_rest(n,m))%P; } int main(){ Init(); int n,m; scanf("%d%d",&n,&m); printf("%d\n",calc(n,m)); return 0; }
在这道题中,有几点很值得学习.
I.由利用组合数计数引出思路.
II.在网格计数问题中利用矩形框和gcd.
III.利用容斥计数.
IV.对于图形,利用对称计数.
V.利用图形对称来计算贡献.
同时我感觉我在解决这道题的时候有许多不足的地方:
I.表现出了在解决计数问题方面的不足,比如调试计数的Kod,以及思维的不严谨,还有对于各种计数方法都不熟悉.
II.第一次学类欧,但是学得很草,证明看了一半,板子还是抄的,而且还不知道板子对不对……
苟利国家生死以， 岂因祸福避趋之。
