UVa11139 Counting Quadrilaterals

题目链接

        UVa11139 Counting Quadrilaterals

题意

        在图下图 的10×10 网格中，你可以看到5 个四边形（注意，四边形的4 条边不能相交，而且没有三点共线）。

        

        当然，还可以连出其他网格四边形。编程统计出n×n（n≤120）网格里可以连出多少个网格四边形。如n=2 时有94 个，n=10 时有12046294 个。

分析

        四边形可按凹凸分类，先考虑凹四边形，它必然是一个三角形和其内一点形成的（见下图），并且三角形内任意一点与三角形的三个顶点能组成三个凹四边形，因此统计凹四边形可以枚举格点三角形和其内格点数（具体需要用到皮克定理 Pick's theorem）。

        

        枚举格点三角形又可以按照三角形和其外包矩形的关系来进行：

        

1. 三角形有两个顶点是外包矩形对角线的端点，第三个顶点为不在此对角线上的任意一点；

2. 三角形有两个顶点是外包矩形某条边的端点，第三个顶点为此边的平行边上不在端点的任意一点；

3. 三角形一个顶点是外包矩形的角顶点，另外两个顶点在与此角不相连的两边上任意一点且非端点；

        设c[m][n]表示以尺寸为为外包矩形的凹四边形计数（后面还会加上凸四边形计数），则对给定的m、n可以按照上面三种情况累加计数。

        接下来考虑凸四边形和其外包矩形的关系：

        

        和凹四边形的分类相似，但多出一种情况d：四个顶点各在一条边上。

        凸四边形的d类计数的结果为;

        b、c分类的计数复杂一点，首先需要把第四个顶点在外包矩形边上时单独考虑，在外包矩形内部时借助皮克定理可求；

        a分类的计算是最复杂的，对于另外两个端点跨对角线两边（图a所示黄色的四边形）好计算，但是在对角线同侧（图a所示紫色的四边形）时则复杂，下面详细展开：

        

        可以先枚举其中一个点，它与两对角点连线并延伸见上图所示，可见只有连线和延伸线所夹的格点才能与之构成凸四边形，可惜这些点的计数不能使用皮克定理求解了，因为延伸线与外包矩形的边界交点坐标可能不是整数，不构成格点多边形这一前提条件。

        这种可以借助动态规划单独计算出来：设f[w][h][s]表示宽高比为w:h（w、h互质）并且横向宽度为s的斜线下方的格点计数，那么f[w][h][s+1] = f[w][h][s] + ，这里当(s+1)*h % w == 0 时 否则 。

        实际的计算过程及一些细节请参见ac代码。

AC代码

#include <iostream>
using namespace std;

#define N 123
int f[N][N][N], g[N][N], n; long long c[N][N];

int gcd(int a, int b) {
    if (a > b) return gcd(b, a);
    if (a == 0) return b;
    if (a & 1) {
        if (b & 1) return gcd(a, (b-a)>>1);
        return gcd(a, b >> 1);
    } else if (b & 1) return gcd(a >> 1, b);
    return gcd(a >> 1, b >> 1) << 1;
}

void init() {
    for (int i=0; i<N; ++i) for (int j=i; j<N; ++j) g[i][j] = g[j][i] = gcd(i, j);
    for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) {
        f[i][j][0] = 0;
        for (int k=1; k<N; ++k) f[i][j][k] = f[i][j][k-1] + (i ? (k*j%i ? k*j/i : k*j/i-1) : 0);
    }
    for (int x=1; x<N; ++x) for (int y=x; y<N; ++y) {
        long long t = (x+1)*(y+1) - g[x][y] - 1;
        long long cc = t*t/2-1 + (x-1)*(x-2) + (y-1)*(y-2) + (x-1)*(y-1)*(6*(x+y) - 8) + (x-1)*(x-1)*(y-1)*(y-1);
        for (int x1=1; x1<=x; ++x1) {
            int y1 = y*x1; y1 = y1%x ? y1/x : y1/x-1;
            while (y1 >= 0) {
                cc += 2*(f[y1][x1][y1] - f[y-y1][x-x1][y1] - g[y-y1][x-x1]*y1/(y-y1));
                cc += 2*(f[x-x1][y-y1][x-x1] - f[x1][y1][x-x1] - g[x1][y1]*(x-x1)/x1);
                cc += 6*(x1*y-x*y1 + 2 - g[x1][y1] - g[x-x1][y-y1] - g[x][y]);
                --y1;
            }
        }
        for (int x1=1; x1<x; ++x1) {
            cc += (x-2)*(y-1) + 2 - g[x1][y] - g[x-x1][y];
            cc += 3*(x*y + 2 - x - g[x1][y] - g[x-x1][y]);
        }
        for (int y1=1; y1<y; ++y1) {
            cc += (x-1)*(y-2) + 2 - g[x][y1] - g[x][y-y1];
            cc += 3*(x*y + 2 - y - g[x][y1] - g[x][y-y1]);
        }
        for (int x1=1; x1<x; ++x1) for (int y1=1; y1<y; ++y1) {
            cc += 2*(x*y+x1*y1 + 6 - 2*x - 2*y - g[x1][y] - g[x][y1] - g[x-x1][y-y1]);
            cc += 6*(x*y-x1*y1 + 2 - g[x1][y] - g[x][y1] - g[x-x1][y-y1]);
        }
        c[x][y] = c[y][x] = cc;
    }
}

void solve() {
    long long ans = 0;
    for (int x=1; x<=n; ++x) for (int y=1; y<=n; ++y) ans += (n+1-x)*(n+1-y) * c[x][y];
    cout << n << ' ' << ans << endl;
}

int main() {
    init();
    while (cin>>n && n) solve();
    return 0;
}