HDU 4483 HDOJ Lattice triangle

http://acm.hdu.edu.cn/showproblem.php?pid=4483

题目意思：

　　给定一个n*n的广场，问选取三个整点构成三角形的方案数。

今天做hust的练习赛的时候，碰到了UVALive3295，是求n*m的矩阵上的方案数。不过那题n,m没这题大。然后就讲到这题，就做了一下，还是很神奇的一题。

整体思路：总方案数 - 三点成线方案数。

n*n的矩阵，整点数为(n+1)*(n+1)，以下讨论，直接另n = n+1

总方案数C(n*n, 3)

三点构成的线与坐标轴平行，C(n, 3)*n*2

本题重点就在求三点成斜线的方案数。

我们来这么考虑三个点，让三个点的横坐标依次增大把旁边两个点的连线看做矩形的对角线，中间的点就在对角线上。

易知矩形对角线有两条，我们不妨只计算对角线方向为右的方案数，*2就是总方案数。

长为i，宽为j的矩形，对角线上的整点数为gcd(i, j)-1，在正方形中有(n-i)*(n-j)个

于是，让tmp表示三点成右斜线的方案数，有以下代码

tmp=0;
for(int i=2;i<n;i++)
for(int j=2;j<n;j++)
    tmp+=(n-i)*(n-j)*(__gcd(i,j)-1);

这样复杂度为n^2，还不能满足题目要求。不过我们可以在此基础上思考。

我们反过来枚举gcd的值呢？

假设i, j的gcd为k

那么gcd(i/k, j/k)==1，且 i, j<n

令a=i/k, b=j/k

则a,b的pair数 = 1 + 2*(ph(2) + ph(3) + ph(4) + ... + ph((n-1)/k))         ph()为欧拉函数

最前面的1表示 a=1,b=1

其他情况下a==b，a,b是不互质的。

令a<b，枚举a，b的取值就是ph(a)，当然，a, b顺序可以换，所以要最前面要*2

而这么多的a,b计算

   sum((n-i)*(n-j)*(__gcd(i, j)-1))

= sum((n - k*a) * (n - k*b) * (k-1))

= (k-1)*sum(n^2 - n*k*(a+b) + k^2*a*b)

n^2的系数就是a,b的pair数

n*k的系数就是这么多对a,b的a+b的总和

k^2的系数就是这么多对a,b的a*b的总和

插播一段数论

　　对于n>2, 若存在一个m，使得gcd(n,m)==1，则gcd(n, n-m)==1。反证法可证。所以可以理解为和n互素的数总是成对出现，和为n。

　　所以[1, n]中与n互质的数的总和为ph(n)*n/2

　　而对于n>2的偶数，gcd(n/2, n)!=1，不用讨论。n==2时，总和为1，恰好也满足这个式子。

于是上面三个系数，都可以通过欧拉函数求出。

本题解决。

代码如下：

#include<cstdio>
using namespace std;
typedef long long ll;
const int N=100005;
const ll md=1000000007,md2=md*6;
int n,ph[N],sc[N],ss[N],sm[N];
void gph(){
    ph[1]=sc[1]=sm[1]=1;
    ss[1]=2;
    for(int i=2;i<N;i++)if(!ph[i])
        for(int j=i;j<N;j+=i){
            if (!ph[j]) ph[j]=j;
            ph[j]-=ph[j]/i;
        }
    for(int i=2;i<N;i++){
        sc[i]=(sc[i-1]+(ll)ph[i]*2)%md;
        ss[i]=(ss[i-1]+(ll)i*ph[i]*3)%md;
        sm[i]=(sm[i-1]+(ll)i*ph[i]%md*i)%md;
    }
}
int main()
{
    gph();
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%d",&n);n++;
        ll tmp=(ll)n*n%md,ans;
        ans=tmp*(tmp-1)%md2*(tmp-2)%md2/6;
        ans-=tmp*(n-1)%md2*(n-2)*2%md2/6;
        tmp=0;
        for(int k=2;k<n;k++){
            int m=(n-1)/k;
            tmp=(tmp+(k-1)*((ll)sc[m]*n%md*n%md - (ll)n*k%md*ss[m]%md + (ll)k*k%md*sm[m]%md))%md;
        }
        ans-=tmp*2;
        printf("%I64d\n",(ans%md+md)%md);
    }
    return 0;
}