X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fkfunc.c;h=a637b6ca70afd768e736117b0dbbc377d36d2318;hb=f69357a6a9e701f12d1f01c52cc0d2893bb2ac7a;hp=b534a31e597d9255b59264304d7b95deec4253b2;hpb=4c6157aec247da8e65994231dccb312b021ba069;p=samtools.git

diff --git a/bcftools/kfunc.c b/bcftools/kfunc.c
index b534a31..a637b6c 100644
--- a/bcftools/kfunc.c
+++ b/bcftools/kfunc.c
@@ -112,16 +112,18 @@ double kf_gammaq(double s, double z)
 }
 
 /* Regularized incomplete beta function. The method is taken from
- * Numerical Recipe in C, 2nd edition, section 6.4. The following page
- * calculates the incomplete beta function, which equals
+ * Numerical Recipe in C, 2nd edition, section 6.4. The following web
+ * page calculates the incomplete beta function, which equals
  * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
  *
  *   http://www.danielsoper.com/statcalc/calc36.aspx
  */
 static double kf_betai_aux(double a, double b, double x)
 {
-	double C, D, f, z;
+	double C, D, f;
 	int j;
+	if (x == 0.) return 0.;
+	if (x == 1.) return 1.;
 	f = 1.; C = f; D = 0.;
 	// Modified Lentz's algorithm for computing continued fraction
 	for (j = 1; j < 200; ++j) {
@@ -138,8 +140,7 @@ static double kf_betai_aux(double a, double b, double x)
 		f *= d;
 		if (fabs(d - 1.) < KF_GAMMA_EPS) break;
 	}
-	z = (x == 0. || x == 1.)? 0. : exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x));
-	return z / a / f;
+	return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
 }
 double kf_betai(double a, double b, double x)
 {