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
",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
",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
",calc(n,m));
  return 0;
}

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