HDU4483 Lattice triangle 数论

题意：给定一个N*N个网格（每条边上共有N+1个点），从这个网格中取出三个点构成三角形，问一共可以构成多少种三角形。

解法：首先令N' = N+1，那么不考虑直线相交的情况下从N'*N'个点中选出三个点的方案数为C(3, N'*N')；然后考虑到每条平行于水平和垂直线的线段上共有2*N'*C(3, N')种情况需要减去，最后还要减去斜线直线上的情况。

斜线上则只考虑从左上到右下的情况，乘以2表示从右到左是对称的。对于每一种斜线的情况都可以将左上端点和右下端点固定，然后就可以看作是一个矩形对角线上除去两个端点外，中间还有多少个点的问题了，对于这个问题，在求凸包上有多少个整数点时应该知道，假设边为(a, b)，线段上点的个数为gcd(a, b) + 1，除去两个端点后为gcd(a, b) - 1，因此枚举出所有的矩形就可以了。

对于边长为(a, b)的矩形，N*N的矩形中共有(N' - a) * (N' - b)个。因此边长为(a, b)的情况下，共有(N' - a) * (N' - b) * (gcd(a, b) - 1)中不合法的情况需要减去。但是该题由于N的范围过大，因此在枚举边长是达到了O(N^2)的时间复杂度，肯定会超时的。因此需要优化这个求解的过程：

具体做法是枚举gcd值，由于gcd值至少为2才会对结果有影响，因此gcd值就可以从2开始枚举起，且这个gcd值最大为N。若枚举了gcd值为k，那么有gcd(a, b) = k，则gcd(a/k, b/k) = 1，因此只要知道在某一范围内互质数的对数，然后再将这些互质的数均乘上一个x即为所求，约定a <= b，由于a、b的取值均为1-N，因此b的取值最大为N / x，那么互质数的对数就是2 * (phi[1] + phi[2] + ... + phi[N / x])，乘以2是因为a，b可以交换位置，后面的求和是因为这么些互质的数均满足要求。那么对于每一组互质的数m, n，有m = a/k, b = b/k，即a = k*m, b = k*n，因此代入到(N' - a) * (N' - b) * (gcd(a, b) - 1)有(N'*N' - (m + n)*N'*k + m*n*k*k) * (k - 1)。这是对于一组来说，由于枚举k值，因此k值相同的数对应一起计算，因此可知所有k值相同的数对，N'*N'的系数为前面所讨论的总对数，N'*k以及k*k的系数由互质对的成对情况分别能够求出递推式。成对是指若gcd(y, N) = 1，那么gcd(N-y, N) = 1 (y < N)。

维护好上面提到的三个系数便可以计算出最后的结果了，计算过程中应防止溢出，代码如下：

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int MaxN = 100000;
const int MOD = int(1e9+7);
int phi[MaxN+5];
int spair[MaxN+5];
int sadd[MaxN+5];
int smul[MaxN+5];

void Euler() {
    phi[1] = 1;
    spair[1] = 1, sadd[1] = 2, smul[1] = 1;
    for (int i = 2; i <= MaxN; ++i) {
        if (!phi[i]) {
            for (int j = i; j <= MaxN; j += i) {
                if (!phi[j]) phi[j] = j;
                phi[j] -= phi[j] / i;
            }
        }
        // 由于于i互质的数一定是互补的即gcd(x, n) = 1，则gcd(n-x, n) = 1 
        spair[i] = (1LL * spair[i-1] + phi[i] * 2) % MOD;
        sadd[i]  = (1LL * sadd[i-1] + 1LL * phi[i] * 3 % MOD * i) % MOD;
        smul[i]  = (1LL * smul[i-1] + 1LL * phi[i] * i % MOD * i) % MOD;
    }
}

LL cal(LL x[]) {
    // 因为有除法，因此要将除数先进行处理
    for (int i = 0; i < 3; ++i) {
        if (x[i] % 2 == 0) {
            x[i] /= 2;
            break;
        }
    }
    for (int i = 0; i < 3; ++i) {
        if (x[i] % 3 == 0) {
            x[i] /= 3;
            break;
        }
    }
    for (int i = 0; i < 3; ++i) x[i] %= MOD;
    return x[0]%MOD*x[1]%MOD*x[2]%MOD;
}

int main() {
    Euler();
    int T, N;
    scanf("%d", &T);
    while (T--) {
        LL ret, exp = 0;
        scanf("%d", &N); // 坐标的取值范围[0,N]，共有N+1种选择 
        ++N;
        LL a[3] = {1LL*N*N, 1LL*N*N-1, 1LL*N*N-2};
        LL b[3] = {N, N-1, N-2};
        ret = cal(a);
        ret = ((ret - cal(b)*2%MOD*N%MOD) % MOD + MOD) % MOD;
        // 也可将mod乘以6，先做乘法然后mod这个新的数，相当于把余数扩大了6倍，然后再除以6 
        for (int k = 2; k < N; ++k) {
            int t = (N-1) / k;
            exp = (exp + 1LL*(k-1)*(1LL*spair[t]*N%MOD*N%MOD))%MOD;
            exp = ((exp - 1LL*(k-1)*(1LL*sadd[t]*k%MOD*N%MOD)%MOD)%MOD + MOD) % MOD;
            exp = (exp + 1LL*(k-1)*k%MOD*k%MOD*smul[t]%MOD)%MOD;
        }
        ret = ((ret-2*exp)%MOD+MOD)%MOD;
        printf("%I64d\n", ret);
    }
    return 0;
}