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=;
const int N=;
int gcd[N][N];
inline int GCD(int x,int y){
return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==?y:GCD(y%x,x));
}
inline void Init(){
register int i,j;
for(i=;i<=;++i)
for(j=;j<=;++j)
gcd[i][j]=GCD(i,j);
}
inline int count_all(int n,int m){
int cnt=(n+)*(m+);
int ret=41666667LL*cnt%P*(cnt-)%P*(cnt-)%P*(cnt-)%P;
register int i,j,tmp,sum;
for(i=;i<=n;++i)
for(j=!i;j<=m;++j){
tmp=gcd[i][j]-;
if(!tmp)continue;
sum=(n-i+)*(m-j+)*(i&&j?:);
ret=(ret-(LL)tmp*sum*(cnt-)+P)%P;
if(tmp==)continue;
ret=(ret+(LL)tmp*(tmp-)/%P*sum%P*%P)%P;
}
return ret;
}
inline int count_rest(int n,int m){
int ret=,sum,s,tmp;
register int i,j,x,y,temp;
for(i=;i<=n;++i)
for(j=;j<=m;++j){
sum=(n-i+)*(m-j+);
s=i*j+;
tmp=(s-j-i-gcd[i][j])*;
for(x=;x<i;++x)
for(y=;y<j;++y){
temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]);
tmp=(tmp+temp*)%P;
}
x=i;
for(y=;y<j;++y){
temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]);
tmp=(tmp+temp*)%P;
}
y=j;
for(x=;x<i;++x){
temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]);
tmp=(tmp+temp*)%P;
}
x=;
for(y=;y<j;++y){
temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]);
tmp=(tmp+temp)%P;
}
y=;
for(x=;x<i;++x){
temp=s-x*y-(gcd[x][j]+gcd[i][y]+gcd[i-x][j-y]);
tmp=(tmp+temp)%P;
}
for(x=;x<i;++x)
for(y=;y<j;++y){
if(y*i-x*j<=)continue;
temp=y*i-x*j+-(gcd[x][y]+gcd[i-x][j-y]+gcd[i][j]);
tmp=(tmp+temp*)%P;
}
ret=(ret+(LL)tmp*sum)%P;
}
return ret;
}
inline int calc(int n,int m){
if((n+)*(m+)<)return ;
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 ;
}

Kod

调过之后,我观察我的程序,我想到了把我的程序优化到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=;
const int N=;
int gcd[N+][N+],gcd_s[N+][N+],gcd_t[N+][N+],f[N+][N+];
inline int GCD(int x,int y){
return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==?y:GCD(y%x,x));
}
inline void Init(){
register int i,j;
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[i][j]=gcd[i][j]=GCD(i,j);
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[i][j]+=gcd_s[i][j-];
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_t[i][j]=gcd_s[i][j];
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[j][i]+=gcd_s[j-][i];
}
inline int count_all(int n,int m){
int cnt=(n+)*(m+);
int ret=41666667LL*cnt%P*(cnt-)%P*(cnt-)%P*(cnt-)%P;
register int i,j,tmp,sum;
for(i=;i<=n;++i)
for(j=!i;j<=m;++j){
tmp=gcd[i][j]-;
if(!tmp)continue;
sum=(n-i+)*(m-j+)*(i&&j?:);
ret=(ret-(LL)tmp*sum*(cnt-)+P)%P;
if(tmp==)continue;
ret=(ret+(LL)tmp*(tmp-)/%P*sum%P*%P)%P;
}
ret=(ret+P)%P;
return ret;
}
inline int query(int a,int b,int x,int y){
if(a>x)return ;
if(b>y)return ;
int ret=gcd_s[x][y];
--a,--b;
if(a>=)ret-=gcd_s[a][y];
if(b>=)ret-=gcd_s[x][b];
if(a>=&&b>=)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=,sum,s;
register LL tmp;
register int i,j,x,y;
for(i=;i<=n;++i)
for(j=;j<=m;++j){
sum=(n-i+)*(m-j+);
if(i>j){
ret=(ret+(LL)f[i][j]*sum%P)%P;
continue;
}
s=i*j+;
tmp=s-j-i-gcd[i][j];
tmp+=s*(i-1LL)*(j-1LL)-(i)*(i-)/2LL*(j)*(j-)/2LL%P;
tmp-=query(,j,i-,j)*(j-1LL)%P+query(i,,i,j-)*(i-1LL)%P+query(,,i-,j-);
tmp+=s*(j-)-(i)*(j)*(j-)/;
tmp-=gcd[i][j]*(j-1LL)+query(i,,i,j-)+query(,,,j-);
tmp+=s*(i-)-(j)*(i)*(i-)/;
tmp-=gcd[i][j]*(i-1LL)+query(,j,i-,j)+query(,,i-,);
for(x=;x<i;++x){
y=x*j/i+;
if(y>=j)break;
tmp+=((j-+y)*(j-y)*i>>)-(x*j-+gcd[i][j])*(j-y)-query_(x,y,j-)-query_(i-x,,j-y);
}
tmp*=;
tmp+=s*(j-);
tmp-=gcd[][j]*(j-1LL)+query(i,,i,j-)*2LL;
tmp+=s*(i-);
tmp-=gcd[i][]*(i-1LL)+query(,j,i-,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+)*(m+)<)return ;
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 ;
}

Kod

这样仍然是不够优美的,因为,即使加入了优化,在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=;
const int P=;
int gcd[N+][N+],gcd_s[N+][N+],f[N+][N+];
bool vis[N+][N+];
inline int GCD(int x,int y){
return gcd[x][y]?gcd[x][y]:gcd[x][y]=(x==?y:GCD(y%x,x));
}
inline void Init(){
register int i,j;
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[i][j]=gcd[i][j]=GCD(i,j);
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[i][j]+=gcd_s[i][j-];
for(i=;i<=N;++i)
for(j=;j<=N;++j)
gcd_s[j][i]+=gcd_s[j-][i];
}
inline int count_all(int n,int m){
//计算所有可能的四个点
int cnt=(n+)*(m+);
int ret=41666667LL*cnt%P*(cnt-)%P*(cnt-)%P*(cnt-)%P;
register int i,j,tmp,sum;
for(i=;i<=n;++i)
for(j=!i;j<=m;++j){
tmp=gcd[i][j]-;
if(!tmp)continue;
sum=(n-i+)*(m-j+)*(i&&j?:);
ret=(ret-(LL)tmp*sum*(cnt-)+P)%P;
if(tmp==)continue;
ret=(ret+(tmp*(tmp-1LL)>>)*sum*)%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 ;
int ret=gcd_s[x][y];
--a,--b;
if(a>=)ret-=gcd_s[a][y];
if(b>=)ret-=gcd_s[x][b];
if(a>=&&b>=)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=;
return ret;
}
if(a>=c||b>=c){
ret=lo(a%c,b%c,c,n);
ret.h+=n*(n+)*(*n+)/*(a/c)*(a/c)+(n+)*(b/c)*(b/c)+*(b/c)*ret.f+*(a/c)*ret.g+(a/c)*(b/c)*n*(n+);
ret.g+=n*(n+)*(*n+)/*(a/c)+n*(n+)/*(b/c);
ret.f+=n*(n+)/*(a/c)+(n+)*(b/c);
return ret;
}
LL m=(a*n+b)/c;
Data tmp=lo(c,c-b-,a,m-);
ret.f=n*m-tmp.f;
ret.g=(m*n*(n+)-tmp.f-tmp.h)/;
ret.h=m*(m+)*n-*tmp.g-*tmp.f-ret.f;
return ret;
}
inline int count_rest(int n,int m){
//计算内凹四边形做出的额外贡献
int ret=,sum,s,cnt;
Data get;
register LL tmp;
register int i,j;
for(i=;i<=n;++i)
for(j=;j<=m;++j){
sum=(n-i+)*(m-j+);
//这个框出现的次数
if(vis[i][j]){
ret=(ret+(LL)f[i][j]*sum%P)%P;
continue;
//用于剪枝
}
s=i*j+;
//初始化pick的部分内容
//以下所说的顶点均指矩形的顶点,所说的占有顶点的三角形均满足其所有顶点都在矩形框上
tmp=s-j-i-gcd[i][j];
//计算占有三个顶点的三角形的贡献
tmp+=s*(i-1LL)*(j-)-(i*(i-1LL)*j*(j-)>>);
//计算占有一个顶点的三角形的面积和
tmp-=query(,j,i-,j)*(j-1LL)+query(i,,i,j-)*(i-1LL)+query(,,i-,j-);
//计算占有一个顶点的三角形的沿边点数和
tmp+=s*(j-)-(i*j*(j-)>>);
//计算第一类占有两个顶点的三角形的面积和
tmp-=gcd[i][j]*(j-)+query(i,,i,j-)+query(,,,j-);
//计算第一类占有两个顶点的三角形的沿边点数和
tmp+=s*(i-)-(j*i*(i-)>>);
//计算第二类占有两个顶点的三角形的面积和
tmp-=gcd[i][j]*(i-)+query(,j,i-,j)+query(,,i-,);
//计算第二类占有两个顶点的三角形的沿边点数和
cnt=((i-)*(j-)-(gcd[i][j]-))>>;
//计算以对角线为分界线把矩形分为两部分后,其中一部分内部的点数
tmp-=(gcd[i][j]-)*cnt+query(,,i-,j-)-(gcd[i][j]*(gcd[i][j]-)>>);
//计算一边为对角线一点在矩形内部的三角形的沿边点数和
tmp*=;
//将由对称产生的贡献算入
get=lo(j,,i,i-);
tmp+=*j*get.g-i*get.h-i*get.f;
//计算一边为对角线一点在矩形内部的三角形的面积和
tmp+=s*(j-);
//计算第三类占有两个顶点的三角形的面积和
tmp-=j*(j-)+(query(i,,i,j-)<<);
//计算第三类占有两个顶点的三角形的沿边点数和
tmp+=s*(i-);
//计算第四类占有两个顶点的三角形的面积和
tmp-=i*(i-)+(query(,j,i-,j)<<);
//计算第四类占有两个顶点的三角形的沿边点数和
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+)*(m+)<)return ;
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 ;
}

Kod

在这道题中,有几点很值得学习.
I.由利用组合数计数引出思路.
II.在网格计数问题中利用矩形框和gcd.
III.利用容斥计数.
IV.对于图形,利用对称计数.
V.利用图形对称来计算贡献.
同时我感觉我在解决这道题的时候有许多不足的地方:
I.表现出了在解决计数问题方面的不足,比如调试计数的Kod,以及思维的不严谨,还有对于各种计数方法都不熟悉.
II.第一次学类欧,但是学得很草,证明看了一半,板子还是抄的,而且还不知道板子对不对……

上一篇:1-3年的Android开发工程师看过来,值得收藏!


下一篇:EF 通过时间戳实现自带 乐观锁