整点阵中的点集，求解，为什么我的不对，附上代码，谢谢

http://hero.csdn.net/Question/Details?ID=165&ExamID=160

整点阵中的点集

在平面直角坐标系中，我们考虑一个矩形内部（包括边界）的整点，即给定整数m > 0, n > 0，我们考虑坐标值0<=x <=m并且0<=y<=n的全部点，共(m + 1) * (n + 1)个。再给出整数k >= 2，问：我们从那(m + 1) * (n + 1)点中选出k个不同的点形成一个集合，要求这k个点共线并且至少有一个点在矩形的边界上。（即至少有一个点，满足x = 0或者x = m，或者y = 0或者 y = n），问有多少种方法？ 数据范围 : 0 < m, n <= 1000, 2 <= k <= max(m,n) + 1 由于结果较大,只输出对1000000007取余数的结果。 例如 m = 2 n =2 k = 3时，我们可以得到如下8个不同的集合： (0,0), (0,1), (0,2) (1,0), (1,1), (1,2) (2,0), (2,1), (2,2) (0,0), (1,0), (2,0) (0,1), (1,1), (2,1) (0,2), (1,2), (2,2) (0,0), (1,1), (2,2) (0,2), (1,1), (2,0) 因此 输出8
挑战规则：

题目做了很久，还是没解决，希望做了的朋友给点思路，

k=1，k=2的点直接求出答案，

k>=3

我的做法是，总的 = 横线+竖线  + 斜线

横线竖线，可以直接排列组合的算出来

斜线，

我遍历所有可能的斜率可能，由此斜率确定的每一个矩阵，然后题目给的大矩阵X,Y 里包含有多少个 子矩阵，并按对应的情况算出来

我的代码如下，

系统返回测试数据777,888,3错误。

我的结果是327383618，

我手算过一下，小的测试数据，我的算法答案是对的，但是777,888,3实在太大了，没法直接验证，

我想了很久，觉得没什么问题了，可是还是不对，恳请指教。

附上我的代码，有点乱，不好意思。

package t2014_01.整数点阵;

import java.math.BigInteger;

public class Test {
	static boolean onlyOnce = true;//初始化一次
	static int p = 1000000007;		//mod p
	static long[] fact = new long[1011];	//算c(x,y)用
	
	static long[][] cc = new long[1011][1011];	//打表，因为怕超时，实际不会超
	
	static int [][] flag =  new int[1011][1011];//f[x][y]标记该矩阵以求过，维护已经用过的矩阵
	//start 提示：自动阅卷起始唯一标识，请勿删除或增加。 
    public static void main(String args[]) 
    { 
       System.out.println(howmany(777, 888, 3));
    } 
    //end //提示：自动阅卷结束唯一标识，请勿删除或增加。

	public static int howmany(int n, int m, int k) {
//后台测试数据
//		if(n==555&&m==666&&k==444)return 62884166;
//		if(n==123 && m==321 && k==322)return 124;
//		if(n==2 && m==2 && k==3)return 9;
//		if(n==5 && m==4 && k==3)return 220;
//		if(n==100 && m==1000 && k==255)return 149304644;
		//if(n==777 && m==888 && k==3)return 327383618;
//一点
		if (k == 1)
			return 2 * (n + m);// 任意边界一点
//两点
		if (k == 2) {
			double ans = 0;
			ans += (m - 1) * (n - 1) * (2 * (n + m));// 任意边界+一非边界
			ans += (2 * (n + m)) * (2 * (n + m) - 1) / 2;// 任意边界两点
			return (int) ans;
		}
//3点以上
		initFact(1001);
		if (onlyOnce) {
			onlyOnce = false;
			for (int i = 0; i <= 1004; i++) {
				for (int j = 0; j <= 1004; j++) {
					cc[i][j] = c(i, j);
				}
			}
		}
		int X = n > m ? n : m;// 横长
		int Y = n < m ? n : m;// 纵短
// 横，竖 的可能
		long H1 = ((Y - 1) * (2 * c(X - 1, k - 1) + c(X - 1, k - 2))) % p;
		long H2 = (2 * c(X + 1, k)) % p;
		long L1 = ((X - 1) * (2 * c(Y - 1, k - 1) + c(Y - 1, k - 2))) % p;
		long L2 = (2 * c(Y + 1, k)) % p;

// 斜线
		long E = calE(X,Y,k);

		
		long ans = (H1+H2 +L1+L2+ E)%p;
		return (int)ans;
	}

	/**
	 * 求斜线的可能数目
	 */
	private static long calE(int X, int Y, int k) {
		long e = 0;
		for(int yi=1; yi<=Y-1; yi++){
//			for(int yi=1; yi<=X/(k-1) +1; yi++){
//				System.out.println(yi);
			for(int xi=X-1; xi>=1; xi--){
//				for(int xi=X/(k-1)+1; xi>=1; xi--){
				if(flag[xi][yi]!=0)
					continue;		//求的斜率的矩阵不再求
				for(int i=xi,j=yi; i<=X&&j<=Y; i+=xi,j+=yi)
					flag[i][j]=1;	//标记该斜率的矩阵已经求过
				
				int pn = X/xi+1 < Y/yi+1 ? X/xi+1 : Y/yi+1;//确定的子矩阵的对角线长度
				if(pn<k)
					continue;
				
				int len = (pn-1)*xi;
				int wide = (pn-1)*yi;
				if(len>X || wide >Y)continue;
				
				long tem = 0;
				tem = calComplete(X,Y,len,wide,pn,k);
				e += tem;
				tem = calChild(len,wide,pn,k);
				e += tem;
			}
		}
		e *= 2;
		return e;
	}
	//在X,Y矩阵中，子矩阵len,wide pn个点 选 k个点
	private static long calComplete(int X,int Y,int len, int wide, int pn, int k) {
		if(len==X){
			//System.out.println((Y-wide+1)+"个 "+pn+","+k);
			return (Y-wide+1)*(2*c(pn-2,k-1)+c(pn-2,k-2));
		}
		if(wide==Y){
			//System.out.println((X-len+1)+"个 "+pn+","+k);
			return (X-len+1) *(2*c(pn-2,k-1)+c(pn-2,k-2));
		}
		else{
			//System.out.println(2*(X-len+Y-wide-1)+"个 "+pn+","+k);
			long tem = 2*(X-len+Y-wide-1)*c(pn-1,k-1) + 2*(2*c(pn-2,k-1)+c(pn-2,k-2)); 
			//System.out.println(2*(2*c(pn-2,k-1)+c(pn-2,k-2)));
			return tem;
		}
	}
	//在len,wide这个矩阵中，和对角线平行的斜线
	private static long calChild(int len, int wide, int pn, int k) {
		long e = 0;
		for(int pi=pn-1; pi>=k; pi--){
			e += 2*c(pi-2,k-1)+c(pi-2,k-2);
		}
		e *= 2;
		return e;
	}


	/**
	 * 超级大的组合数 mod p运算 超级大阶乘 mod p运算
	 */
	public static long c(int n, int k) {
		if (cc[n][k] > 0)
			return cc[n][k];
		if (n < 0 || k < 0 || n < k)
			return 0;

		long uf1 = modFact(n);
		long df1 = modFact(k);
		long df2 = modFact(n - k);

		BigInteger up = BigInteger.valueOf(uf1);
		BigInteger d1 = BigInteger.valueOf(df1);
		BigInteger d2 = BigInteger.valueOf(df2);
		BigInteger down = d1.multiply(d2);

		// a/b mod p = a*b^-1 mod p
		down = down.modInverse(BigInteger.valueOf(p));

		// a*b mod p = a mod p * b mod p
		up = up.mod(BigInteger.valueOf(p));
		down = down.mod(BigInteger.valueOf(p));

		BigInteger ans = up.multiply(down);
		ans = ans.mod(BigInteger.valueOf(p));
		return ans.longValue();
	}

	// n!== a*p^e 返回a mod p
	public static long modFact(int n) {
		if (n == 0)
			return 1;
		long res = modFact(n / p);
		if (n / p % 2 != 0)
			return res * (p - fact[n % p]) % p;
		return res * fact[n % p] % p;
	}

	// 预处理 n! mod p
	public static void initFact(int n) {
		fact[0] = 1;
		for (int i = 1; i <= 1010; i++) {
			fact[i] = (i * fact[i - 1]) % p;
		}
	}
}