[求助]弄的一串t分布的代码，可计算结果不对，请高手帮忙看下。

http://www.qmss.jp/e-stat/stattab/t-distribution.xls（t分布表）
代码：
//////////////////////////////////////////////////////////////
//功能：计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述：gamma[x]=Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用：
double lagam(double x);

//////////////////////////////////////////////////////////////
//功能：不完全贝塔(beta)函数
//描述：Bx[a,b]=Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}]/B[a,b]
//描述：B[a,b]=gamma[a]*gamma[b]/gamma[a+b]
//参数：a-参数,b-参数
//调用：lagam();
double lhbeta(double a,double b,double x);

//////////////////////////////////////////////////////////////
//功能：t-分布函数
//参数：n-自由度
//调用：lhbeta(),lagam();
double ljstd(double t,int n);

#include "stdio.h"
#include "math.h"
//功能：计算伽马(gamma)函数值,gamma函数积分区间为0到正无穷
//描述：Integrate[exp[-t]*t^(x-1),{t,0,∞}]
//调用：
double lagam(double x)
{
int i;
double y,t,s,u;
static double a[11]={ 0.0000677106,-0.0003442342,
0.0015397681,-0.0024467480,0.0109736958,
-0.0002109075,0.0742379071,0.0815782188,
0.4118402518,0.4227843370,1.0};
if (x<=0.0)
{
printf("err**x<=0!\n");
return(-1.0);
}
y=x;
if (y<=1.0)
{
t=1.0/(y*(y+1.0));
y=y+2.0;
}
else if (y<=2.0)
{
t=1.0/y;
y=y+1.0;
}
else if (y<=3.0)
{
t=1.0;
}
else
{
t=1.0;
while (y>3.0)
{
y=y-1.0;
t=t*y;
}
}
s=a[0];
u=y-2.0;
for (i=1; i<=10; i++)
{
s=s*u+a[i];
}
s=s*t;
return(s);
}

----------------------------------------------
//功能：不完全贝塔(beta)函数
//描述：Bx[a,b]=Integrate[t^(a-1)*(1-t)^(b-1),{t,0,x}]/B[a,b]
//描述：B[a,b]=gamma[a]*gamma[b]/gamma[a+b]
//参数：a-参数,b-参数
//调用：lagam();
double lhbeta(double a,double b,double x)
{
double y;
if (a<=0.0)
{
printf("err**a<=0!");
return(-1.0);
}
if (b<=0.0)
{
printf("err**b<=0!");
return(-1.0);
}
if ((x<0.0)||(x>1.0))
{
printf("err**x<0 or x>1 !");
return(1.0e+70);
}
if ((x==0.0)||(x==1.0))
{
y=0.0;
}
else
{
y=a*log(x)+b*log(1.0-x);
y=exp(y);
y=y*lagam(a+b)/(lagam(a)*lagam(b));
}
if (x<(a+1.0)/(a+b+2.0))
{
y=y*beta(a,b,x)/a;
}
else
{
y=1.0-y*beta(b,a,1.0-x)/b;
}
return(y);
}

static double beta(double a,double b,double x)
{ int k;
double d,p0,q0,p1,q1,s0,s1;
p0=0.0; q0=1.0; p1=1.0; q1=1.0;
for (k=1; k<=100; k++)
{
d=(a+k)*(a+b+k)*x;
d=-d/((a+k+k)*(a+k+k+1.0));
p0=p1+d*p0;
q0=q1+d*q0;
s0=p0/q0;
d=k*(b-k)*x;
d=d/((a+k+k-1.0)*(a+k+k));
p1=p0+d*p1;
q1=q0+d*q1;
s1=p1/q1;
if (fabs(s1-s0)<fabs(s1)*1.0e-07)
{
return(s1);
}
}
printf("a or b too big !");
return(s1);
}
-----------------------------------------------------------------------------
//功能：t-分布函数
//参数：n-自由度
//调用：lhbeta(),lagam();
double ljstd(double t,int n)
{
double y;
if (t<0.0)
{
t=-t;
}
y=1.0-lhbeta(n/2.0,0.5,n/(n+t*t));
return(y);
}

----------------解决方案--------------------------------------------------------
没人会吗？？？
毕业设计要用阿！！！
有一点想法也行，说说阿！！
大家一起讨论。
----------------解决方案--------------------------------------------------------

看见这些东西就头晕
计算太多了
还没有主函数

----------------解决方案--------------------------------------------------------
太深奥了！！
----------------解决方案--------------------------------------------------------

迷糊了
----------------解决方案--------------------------------------------------------