矩形内整点直角三角形有多少个

Question

小弟一个问题思考甚久，未果，忘各位高手指点：

给定一个坐标范围N(就是说0<=X,Y<=N)，怎样计算三个顶点都在坐标范围内并且顶点为整点的直角三角形个数？

如果实在没有数学公式的话，给出一个时间复杂度不超过（N*N*logN)的算法也可以，谢谢

急，在线等

样例：
N=10 总数：23596
N=20 总数：418716
N=30 总数：2288304
显然... 包括那种三边都不平行于坐标轴的，否则也不值得悬赏两百分求解啊~~~~~

只要顶点是整点且在坐标范围内就可以了

Answer 1

公式我推不出，算法倒是有一个：

设矩形X轴Y轴上长度分别是M，N。

两直角边平行于坐标轴时，确定两个X坐标M（M+1)/2，确定两个Y坐标
N（N+1）/2，可以得到4个直角三角形。于是共有MN（M+1)（N+1)个。

两直角边与坐标轴不平行时，先找两个点(x1,y1)(x2,y2)确定一条斜率为正的直角边（二重循环），然后要求找(x3,y3)，那么有|x2-x1|/|y2-y1|=（|y3-y1|/|x3-x1|或|y3-y2|/|x3-x2|），因为对于互补的角A和B有tgA=ctgB。

那么如何求第三个点呢？我想了一种方法，举例如下
设矩形左下为原点（0，0），右上为（100，200）。已知确定一条直角边两端点（6，7）和（10，13），现在要找第三个点，可能为（6-m，7+n）、(6+m，7-n）、（10-m，13+n）、(10+m，13-n）。那么就要找满足n/m=(10-6)/(13-7)=4/6的整数对（m,n）。先用辗转相除求出4和6的最大公因子，然后相约得到互质的整数对(2,3)，则(m，n）必然是(3,2)、（6，4）、（9，3）...这样的形式。具体有多少个这样的第三个点，对于(10+m,13-n)来说，就是10+m<100,13-n>0，因此为(100-10)/3取整和（13-0）/2取整中较小的那个，其他三种形式同理。

复杂度计算：显然是（M*N*辗转相除法复杂度）

辗转相除法复杂度忘了，我大概推了一下，考虑边长m和n的矩形（m>n），相除一次后变成(n,p)，显然有p<1/2(n+p)=1/2m，新矩形面积减少到原来一半以下，最坏情况一直相除直到2和1的情形，这时候的面积为2，所以相除次数不超过log(mn)=logm+logn。

显然算法中有m<=M，n<=N

因此算法复杂度为（M*N*（logM+logN）），满足楼主要求

修改：

我果然很粗心。。。

我刚刚又思考了一下，可以换个思路，先考虑斜边而不是两直角边，借用你的方法，首先不考虑平移，将斜边中心定为原点，由于原来中心点坐标可能为m/2，m为奇数的情况，因此将坐标均乘以2，确定了一条斜边后，要找到直角顶点（知道直角顶点的情况下确定斜边其实是一样的）。问题本质上就成为以该斜边直径的圆经过多少整点，本来要遍历的点是O（m+n）级别的(对所有不大于直径的x坐标，求y坐标是否整数)，这样整个复杂度就变成了n^3，因此求整点的算法必须改进。

现在的问题就是对于(m,n)，要求出所有满足m^2+n^2=a^2+b^2的(a,b)，最后(n,m)(a,b)(b,a)都是可能满足条件的点(x,y)，“可能”是因为还需要x-m和y-n都是偶数（因为坐标之前被乘以2）。所以事先可以对任意m和n，将(m,n)保存在S数组D的D[m^2+N^2]指向的地址中，这个地址指向的是一个保存了如上的(m,n)(n,m)(a,b)(b,a)...的链表（只统计个数是不行的，因为对于不同的(m,n)，平移得到的三角形个数是不同的，必须保存(m,n)，其实就是一个散列结构），这样对(m,n)求所有(a,b)可以直接查询，不过这样做空间开销很大，效率上不一定比之前求最大公约数高，虽然求公因子理论上为logN复杂度，实际运行平均效率远高于这个最坏值。

最后三角形大小形状角度确定后进行平移计算确实是O(1)复杂度，如楼下所说，因为三角形三顶点最大横坐标之差，最大纵坐标之差已经确定了，也就是说覆盖这个三角形的最小子矩形（边与坐标轴平行）已经确定了，把两个方向上平移的可能相乘就可以了。

感觉还是数学不过关，m^2+n^2=a^2+b^2应该有更好的计算方法的，就是想不出来。

Answer 2

楼上的算法的基本框架差不多了，但是做的不够精细。
首先LS的复杂度计算错误，楼上的算法的复杂度是O(N^4*loglogN*K)，K为一个较大的且依赖N的常数。
枚举直角边的时候不用枚举两个点，只要固定一个点在原点，枚举另外一个点即可，增量的思路也同LS。但是计数的时候直接根据当前三角形的外框直接用乘法原理就可以了。这样的算法复杂度就

变成了O(N^2loglogN*K)，代码如下：
#include <cstdio>
#include <iostream>
using namespace std;

int n;

int gcd(int a,int b){
if (b==0) return a;
return gcd(b,a%b);
}

int main(){
while (scanf("%d",&n)!=EOF){
if (n==0) return 0;
long long ans=0LL;
for (int i=1;i<=n;++i)
for (int j=1;j<=n;++j){
ans+=4LL*(n-i+1)*(n-j+1);
int dx=j/gcd(i,j),dy=i/gcd(i,j);
long long cur=0LL;
for (int k=1;;++k){
int w=i+k*dx;
int h=j>?k*dy;
if (w>n) break;
if (h>n) break;
cur+=(long long)(n-w+1)*(n-h+1);
}
if (i==2&&j==3) cout << cur << endl;
for (int k=1;;++k){
int w=i>?k*dx;
int h=j+k*dy;
if (w>n) break;
if (h>n) break;
cur+=(long long)(n-w+1)*(n-h+1);
}
if (i==2&&j==3) cout << cur << endl;
ans+=cur*2LL;
}
cout << ans << endl;
}
return 0;
}
更进一步，最后计数的时候还可以直接用sigma(i)和sigma(i^2)的求和公式来优化求和的过程，算法复杂度就变成了O(N^2loglogN)，此外题目中的矩形为正方形，还可以根据对称性进行常数优化

，以下为代码：
#include <cstdio>
#include <iostream>
using namespace std;

int gcd(int a,int b);
long long sum2[2010][2010],num[2010][2010];
int f[2010][2010][2],sum[2010][2010],n;

int main(){
for (int i=1;i<=2000;++i)
for (int j=1;j<=2000;++j){
int temp=gcd(i,j);
f[i][j][0]=i/temp;
f[i][j][1]=j/temp;
}
for (int i=1;i<=2000;++i)
for (int j=i;j<=2000;++j){
sum[i][j]=(i+j)*(j-i+1)/2;
sum2[i][j]=((long long)j*(j+1)*(j+j+1)-(long long)i*(i-1)*(i+i-1))/6;
}
while (scanf("%d",&n)!=EOF){
if (n==0) return 0;
long long ans=0LL;
for (int i=1;i<=n;++i){
int di=n-i+1;
for (int j=1;j<=n;++j){
int dj=n-j+1;
ans+=(long long)(4*di*dj);
if (i>j){
ans+=num[j][i]*2LL;
continue;
}
long long cur=0LL;
int dx=f[i][j][1],dy=f[i][j][0];
int range1=(n-i)/dx;
if (n/dy<range1) range1=n/dy;
int range2=j/dy;
int range=range2;
if (range1<range) range=range1;
cur+=(long long)range*di*dj-(long long)dx*dj*sum[1][range];
if (range<range1){
cur+=(long long)(range1-range)*di*(n+1);
cur-=(long long)((n+1)*dx+di*dy)*sum[range+1][range1];
cur+=(long long)dx*dy*sum2[range+1][range1];
}
range1=(n-j)/dy;
if (n/dx<range1) range1=n/dx;
range=range2=i/dx;
if (range1<range) range=range1;
cur+=(long long)range*di*dj-(long long)dy*di*sum[1][range];
if (range<range1){
cur+=(long long)(range1-range)*dj*(n+1);
cur-=(long long)((n+1)*dy+dj*dx)*sum[range+1][range1];
cur+=(long long)dx*dy*sum2[range+1][range1];
}
num[i][j]=cur;
ans+=cur*2LL;
}
}
cout << ans << endl;
}
return 0;
}

int gcd(int a,int b){
if (b==0) return a;
return gcd(b,a%b);
}
这个代码跑n=2000我本地需要0.5s，但其中预处理时间较长，跑15个2000就只要3s不到的时间了。

Answer 3

直接对前两个点循环，然后设第三个点坐标为(x,y),
前两个点坐标为(x1,y1)和(x2,y2)，然后用勾股定理解一个二元二次方程，解得后在判断是否在矩形内，这样的话，第一步循环的复杂度是O(N*N)，第二步解方程可以给出公式，但是要求平方根，不过估计也不会超过O(log2n),这样算法总时间复杂度就会小于O(N*N*log2n)，其实应该就是略大于O(N*N)一点点。你去试试吧，代码肯定很长，我不想写了

这里说明一下，可能会以为前两个点循环要O(N^4)的复杂度，但是其实不用。用高数的方法可以说明如果有一个点固定了，那么我那里要枚举的第二个点的可能位置不超过log2(n)个。（因为只有当第一点、第二点与任一定点的角度β的正弦和正切均为一个有理数域的扩域中的数才行）这样算法复杂度应该是O(n*n*log2n)左右，而log2n和logn差不多（log2n≈3logn），所以我的算法复杂度等价于O(n*n*3logn)，虽然稍微大了一点，但是作一点剪枝应该还是符合要求的。

Answer 4

我来给你一个答案吧，这应该属于一个“复杂网络”问题

没有统一的数学公式！

我是这样证明的：编程计算出：
N 直角三角形
1 4
2 44
3 200
4 596
5 1444
6 2960
7 5520
8 9496
9 15332
10 23596
20 418716
30 2288304
你可以轻易的算出当直角三角形的直角边平行于坐标轴时，共有(N+1)*(N+1)*N*N ～ N^4,然后当直角边是于轴成45度角时应该有～N^3量级个。所以总数应该是一个N^4的多项式，选取任意5个点拟合，但是选取不同的点拟合出来不同的结果，即这个假定的通式对不同N值是不同的，所以没有固定的通式。
证毕

Answer 5

不太好算，就算是用点模型的方法，四个点有四个，六个点有十四个，到十二个点的时候就不太好数了，因为牵扯到这十二个点怎么排的问题了，是2*6，还是3*4。