2. 程式人生 > >Project Euler 453 Lattice Quadrilaterals 困難的計數問題

Project Euler 453 Lattice Quadrilaterals 困難的計數問題

阿新 • • 發佈：2018-03-16
省選模擬賽 我想 dea tdi 初始 wap 存在 hide att

這是一道很綜合的計數問題,對於思維的全面性,解法的過渡性,代碼能力,細節處理,計數問題中的各種算法,像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

Project Euler 453 Lattice Quadrilaterals 困難的計數問題

相關推薦

Project Euler 453 Lattice Quadrilaterals 困難的計數問題

省選模擬賽 我想 dea tdi 初始 wap 存在 hide att 這是一道很綜合的計數問題,對於思維的全面性,解法的過渡性,代碼能力,細節處理,計數問題中的各種算法,像gcd、容斥、類歐幾裏德算法都有考察.在省選模擬賽中做到了這題,然而數據範圍是n,m小於等於1000

project euler 14 collatz

xiang wan urb col site asi pri hang dai Java%E5%AD%97%E8%8A%82%E3%80%81%E5%8D%81%E8%BF%9B%E5%88%B6%E3%80%81%E5%8D%81%E5%85%AD%E8%BF%9B%E5

Project Euler：problem 104

clu int hub printf div pro key urn 鏈接 題目鏈接：http://pe-cn.github.io/104/ # include <stdio.h> # include <stdlib.h> void sort(i

Project Euler：Problem 41 Pandigital prime

emp ont div track using string space number end We shall say that an n-digit number is pandigital if it makes use of all the digits

Project Euler：Problem 34 Digit factorials

class -a post ria urn back ble while [1] 145 is a curious number, as 1! + 4! + 5! = 1 + 24 + 120 = 145. Find the s

[Perl 6][Project Euler] Problem 9 - Special Pythagorean triplet

spc tid int exists post auto gsp none style [Perl 6][Project Euler] Problem 9 - Special Pythagorean triplet Description A Pythagorea

Project Euler 501 Eight Divisors (數論)

sum scan 就是 ont tinc distinct class 套路 pair 題目鏈接： https://projecteuler.net/problem=501 題意： \(f(n)\) be the count of numbers not exceeding

【Project Euler】530 GCD of Divisors 莫比烏斯反演

rac 困難 efi pan opened can tps blank 分塊 【題目】GCD of Divisors 【題意】給定f(n)=Σd|n gcd(d,n/d)的前綴和F(n)，n=10^15。 【算法】莫比烏斯反演 【題解】參考：任之洲數論函數.pdf 這個範圍

project euler 169

ret con href .net IT pac push_back tps tor project euler 169 題目鏈接：https://projecteuler.net/problem=169 參考題解：http://tieba.baidu.com/p/2738

Flea Circus（Project Euler 213）

main bsp make 單位 targe urn %d using oid original version hackerrank programming version 題目大意是N*N的格子，每個格子一開始有1個跳蚤，每過單位時間跳蚤會等概率向四周跳，問M

project euler之最大的素因子

數字 euler class tin urn break 如果 == 素因子 from math import floor def panduan(num): if num > 1: if num > 1: #

project euler之第10001個素數

range spa span for style num count math class import math count = 0 num = 1 while count<10001: num +=1 if num > 1:

project euler之Large sum

range readline file clas pan put NPU inpu txt import fileinput f=open("1.txt","r") m=0 for i in range(100): s=f.readline() m+=in

[Project Euler 429] Sum of squares of unitary divisors（數論）

blog ont namespace href 題解 str return void rap 題目鏈接：https://projecteuler.net/problem=429 題目： A unitary divisor dd of a number nn is a d

Project Euler Problem 46

Problem 46 : Goldbach’s other conjecture It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a

Project Euler Problem 45

Problem 45 : Triangular, pentagonal, and hexagonal Triangle, pentagonal, and hexagonal numbers are generated by the following formulae: Tria

Project Euler Problem 44

Problem 44 : Pentagon numbers Pentagonal numbers are generated by the formula, Pn=n(3n−1)/2. The first ten pentagonal numbers are: 1, 5, 12,

Project Euler Problem 43

Problem 43 : Sub-string divisibility The number, 1406357289, is a 0 to 9 pandigital number because it is made up of each of the digits 0 to 9

Project Euler Problem 42

Problem 42 : Coded triangle numbers The nth term of the sequence of triangle numbers is given by, tn =½n(n+1); so the first ten triangle numbe

Project Euler Problem 41

Problem 41 : Pandigital prime We shall say that an n-digit number is pandigital if it makes use of all the digits 1 to n exactly once. For exa