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;
}

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=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;
}

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=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;
}

Kod

在這道題中,有幾點很值得學習.I.由利用組合數計數引出思路.II.在網格計數問題中利用矩形框和gcd.III.利用容斥計數.IV.對於圖形,利用對稱計數.V.利用圖形對稱來計算貢獻.同時我感受我在解決這道題的時候有許多不足的地方:I.表現出了在解決計數問題方面的不足,好比調試計數的Kod,以及思惟的不嚴謹,還有對於各類計數方法都不熟悉.II.第一次學類歐,可是學得很草,證實看了一半,板子仍是抄的,並且還不知道板子對不對……