comments / cosmetic

git-svn-id: https://svn.r-project.org/R/trunk@88401 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: maechler

src/nmath/bessel_i.c

@@ -231,12 +231,12 @@ static void I_bessel(double *x, double *alpha, int *nb,
     /* Local variables */
     int nend, intx, nbmx, k, l, n, nstart;
     double pold, test,	p, em, en, empal, emp2al, halfx,
-	aa, bb, cc, psave, plast, tover, psavel, sum, nu, twonu;
+	aa, bb, cc, psave, plast, tover, psavel, sum,
 
-    /*Parameter adjustments */
-    --bi;
-    nu = *alpha;
-    twonu = nu + nu;
+	nu = *alpha,
+	twonu = ldexp(nu, 1); // = 2*nu
+
+    --bi; // use 1-indexing below
 
     /*-------------------------------------------------------------------
       Check for X, NB, OR IZE out of range.
@@ -486,7 +486,7 @@ static void I_bessel(double *x, double *alpha, int *nb,
 	    empal = 1. + nu;
 #ifdef IEEE_754
 	    /* No need to check for underflow */
-	    halfx = .5 * *x;
+	    halfx = ldexp(*x, -1); // = *x / 2 = .5 * *x
 #else
 	    if (*x > enmten_BESS) */
 		halfx = .5 * *x;

src/nmath/dbinom.c

@@ -30,7 +30,8 @@
  *   This checks for argument validity, and calls dbinom_raw().
  *
  *   dbinom_raw() does the actual computation; note this is called by
- *   other functions in addition to dbinom().
+ *   many other functions in addition to dbinom():
+ *   dnbinom(), dbeta(), df(), dgeom(), dhyper()
  *     (1) dbinom_raw() has both p and q arguments, when one may be represented
  *         more accurately than the other (in particular, in df()).
  *     (2) dbinom_raw() does NOT check that inputs x and n are integers. This

src/nmath/dgamma.c

@@ -41,7 +41,6 @@
 
 double dgamma(double x, double shape, double scale, int give_log)
 {
-    double pr;
 #ifdef IEEE_754
     if (ISNAN(x) || ISNAN(shape) || ISNAN(scale))
         return x + shape + scale;
@@ -58,6 +57,7 @@ double dgamma(double x, double shape, double scale, int give_log)
 	return give_log ? -log(scale) : 1 / scale;
     }
 
+    double pr;
     if (shape < 1) {
 	pr = dpois_raw(shape, x/scale, give_log);
 	return (

src/nmath/dgeom.c

@@ -30,9 +30,7 @@
 #include "dpq.h"
 
 double dgeom(double x, double p, int give_log)
-{ 
-    double prob;
-
+{
 #ifdef IEEE_754
     if (ISNAN(x) || ISNAN(p)) return x + p;
 #endif
@@ -44,7 +42,7 @@ double dgeom(double x, double p, int give_log)
     x = R_forceint(x);
 
     /* prob = (1-p)^x, stable for small p */
-    prob = dbinom_raw(0.,x, p,1-p, give_log);
+    double prob = dbinom_raw(0.,x, p,1-p, give_log);
 
     return((give_log) ? log(p) + prob : p*prob);
 }

src/nmath/dnbeta.c

@@ -1,7 +1,7 @@
 /*
  *  Mathlib : A C Library of Special Functions
+ *  Copyright (C) 2000-2021 The R Core Team
  *  Copyright (C) 1998 Ross Ihaka
- *  Copyright (C) 2000-12 The R Core Team
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -56,10 +56,6 @@ double dnbeta(double x, double a, double b, double ncp, int give_log)
 {
     const static double eps = 1.e-15;
 
-    int kMax;
-    double k, ncp2, dx2, d, D;
-    LDOUBLE sum, term, p_k, q;
-
 #ifdef IEEE_754
     if (ISNAN(x) || ISNAN(a) || ISNAN(b) || ISNAN(ncp))
 	return x + a + b + ncp;
@@ -74,18 +70,21 @@ double dnbeta(double x, double a, double b, double ncp, int give_log)
     if(ncp == 0)
 	return dbeta(x, a, b, give_log);
 
-    /* New algorithm, starting with *largest* term : */
-    ncp2 = 0.5 * ncp;
-    dx2 = ncp2*x;
-    d = (dx2 - a - 1)/2;
-    D = d*d + dx2 * (a + b) - a;
+    /* Non-central Beta:  New algorithm, starting with *largest* term : */
+    double
+	ncp2 = ldexp(ncp, -1), // = 0.5 * ncp
+	dx2 = ncp2*x,
+	d = ldexp(dx2 - a - 1, -1), // = (...)/2
+	D = d*d + dx2 * (a + b) - a;
+    int kMax;
     if(D <= 0) {
 	kMax = 0;
     } else {
 	D = ceil(d + sqrt(D));
 	kMax = (D > 0) ? (int)D : 0;
     }
 
+    LDOUBLE sum, term, p_k, q;
     /* The starting "middle term" --- first look at it's log scale: */
     term = dbeta(x, a + kMax, b, /* log = */ TRUE);
     p_k = dpois_raw(kMax, ncp2,              TRUE);
@@ -101,7 +100,7 @@ double dnbeta(double x, double a, double b, double ncp, int give_log)
     /* Now sum from the inside out */
     sum = term = 1. /* = mid term */;
     /* middle to the left */
-    k = kMax;
+    double k = kMax;
     while(k > 0 && term > sum * eps) {
 	k--;
 	q = /* 1 / r_k = */ (k+1)*(k+a) / (k+a+b) / dx2;

src/nmath/dnchisq.c

@@ -62,6 +62,8 @@ double dnchisq(double x, double df, double ncp, int give_log)
 
     if(mid == 0) {
 	/* underflow to 0 -- maybe numerically correct; maybe can be more accurate,
+	 * TODO: above log(.):  logmid =  dpois_raw(imax, ncp2, TRUE) +  dchisq(x, dfmid, TRUE);
+	 * ----  and switch to complete log-scale summation of logterm += log(q) {etc} below
 	 * particularly when  give_log = TRUE */
 	/* Use  central-chisq approximation formula when appropriate;
 	 * ((FIXME: the optimal cutoff also depends on (x,df);  use always here? )) */

src/nmath/dnt.c

@@ -81,15 +81,15 @@ double dnt(double x, double df, double ncp, int give_log)
 
     /* If infinite df then the density is identical to a
        normal distribution with mean = ncp.  However, the formula
-       loses a lot of accuracy around df=1e9
+       loses a lot of accuracy around df=1e9 // FIXME?
     */
-    if(!R_FINITE(df) || df > 1e8)
+    if(!R_FINITE(df) || df > 1e8) 
 	return dnorm(x, ncp, 1., give_log);
 
     /* Do calculations on log scale to stabilize */
 
     /* Consider two cases: x ~= 0 or not */
-    if (fabs(x) > sqrt(df * DBL_EPSILON)) {
+    if (fabs(x) > sqrt(df * DBL_EPSILON)) { // |x| > eps * sqrt(df)
 	u = log(df) - log(fabs(x)) +
 	    log(fabs(pnt(x*sqrt((df+2)/df), df+2, ncp, 1, 0) -
 		     pnt(x, df, ncp, 1, 0)));

src/nmath/lbeta.c

@@ -71,8 +71,7 @@ double lbeta(double a, double b)
 	return lgammafn(p) + corr + p - p * log(p + q)
 		+ (q - 0.5) * log1p(-p / (p + q));
     }
-    else {
-	/* p and q are small: p <= q < 10. */
+    else { /* p and q are small: p <= q < 10. */
 	/* R change for very small args */
 	if (p < 1e-306) return lgamma(p) + (lgamma(q) - lgamma(p+q));
 	else return log(gammafn(p) * (gammafn(q) / gammafn(p + q)));

src/nmath/lgamma.c

@@ -100,7 +100,7 @@ double lgammafn_sign(double x, int *sgn)
 #endif
 	    return M_LN_SQRT_2PI + (x - 0.5) * log(x) - x + lgammacor(x);
     }
-    /* else: x < -10; y = -x */
+    /* else: x < -10; y = -x > 10 */
     sinpiy = fabs(sinpi(y));
 
     if (sinpiy == 0) { /* Negative integer argument ===

src/nmath/lgammacor.c

@@ -68,6 +68,11 @@ attribute_hidden double lgammacor(double x)
  *   xbig = 2 ^ 26.5
  *   xmax = DBL_MAX / 48 =  2^1020 / 3 */
 #define nalgm 5
+/*        NB: -- we'd need nalgm = 6 terms for full precision, but the result is
+	  ==     always used in +/- terms of considerably larger size ~ x*log(x)
+   (we could even *decrease* nalgm for larger y)
+*/
+
 #define xbig  94906265.62425156
 
     if (x < 10) // possibly consider stirlerr()