[gnumeric] Mathfuncs: import bessel_j and bessel_y from R.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Mathfuncs: import bessel_j and bessel_y from R.
- Date: Tue, 27 Aug 2013 20:42:27 +0000 (UTC)
commit 97840bc5d96c882c6fe297e235676c1c234cdbe9
Author: Morten Welinder <terra gnome org>
Date: Tue Aug 27 16:41:23 2013 -0400
Mathfuncs: import bessel_j and bessel_y from R.
Not connected yet.
src/mathfunc.c | 1103 +++++++++++++++++++++++++++++++++++++++++++++++++++
src/mathfunc.h | 3 +
src/outoflinedocs.c | 8 +
tools/import-R | 19 +-
4 files changed, 1131 insertions(+), 2 deletions(-)
---
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 9cf05c8..d9db981 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -77,6 +77,7 @@ static inline int imax2 (int x, int y) { return MAX (x, y); }
#define MATHLIB_STANDALONE
#define ML_ERR_return_NAN { return gnm_nan; }
static void pnorm_both (gnm_float x, gnm_float *cum, gnm_float *ccum, int i_tail, gboolean log_p);
+static gnm_float bessel_y_ex(gnm_float x, gnm_float alpha, gnm_float *by);
/* MW ---------------------------------------------------------------------- */
@@ -125,6 +126,29 @@ gnm_acoth (gnm_float x)
return gnm_atanh (1 / x);
}
+gnm_float
+gnm_gamma (gnm_float x)
+{
+ if (gnm_isnan (x))
+ return x;
+
+ if (x > 0) {
+ if (x >= G_MAXINT)
+ return gnm_pinf;
+
+ if (x == gnm_floor (x))
+ return fact ((int)x - 1);
+
+ /* We can probably do better than this. */
+ return gnm_exp (gnm_lgamma (x));
+ } else if (x == gnm_floor (x))
+ return gnm_nan;
+ else {
+ return M_PIgnum /
+ (gnm_sin (x * M_PIgnum) * gnm_gamma (1 - x));
+ }
+}
+
/* ------------------------------------------------------------------------- */
/* --- BEGIN MAGIC R SOURCE MARKER --- */
@@ -4304,6 +4328,580 @@ L230:
}
/* ------------------------------------------------------------------------ */
+/* Imported src/nmath/bessel_j.c from R. */
+/*
+ * Mathlib : A C Library of Special Functions
+ * Copyright (C) 1998-2012 Ross Ihaka and 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
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, a copy is available at
+ * http://www.r-project.org/Licenses/
+ */
+
+/* DESCRIPTION --> see below */
+
+
+/* From http://www.netlib.org/specfun/rjbesl Fortran translated by f2c,...
+ * ------------------------------=#---- Martin Maechler, ETH Zurich
+ * Additional code for nu == alpha < 0 MM
+ */
+
+#ifndef MATHLIB_STANDALONE
+#endif
+
+#define min0(x, y) (((x) <= (y)) ? (x) : (y))
+
+static void J_bessel(gnm_float *x, gnm_float *alpha, long *nb,
+ gnm_float *b, long *ncalc);
+
+gnm_float bessel_j(gnm_float x, gnm_float alpha)
+{
+ long nb, ncalc;
+ gnm_float na, *bj;
+#ifndef MATHLIB_STANDALONE
+ const void *vmax;
+#endif
+
+#ifdef IEEE_754
+ /* NaNs propagated correctly */
+ if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
+#endif
+ if (x < 0) {
+ ML_ERROR(ME_RANGE);
+ return gnm_nan;
+ }
+ na = gnm_floor(alpha);
+ if (alpha < 0) {
+ /* Using Abramowitz & Stegun 9.1.2
+ * this may not be quite optimal (CPU and accuracy wise) */
+ return(bessel_j(x, -alpha) * gnm_cos(M_PIgnum * alpha) +
+ ((alpha == na) ? 0 :
+ bessel_y(x, -alpha) * gnm_sin(M_PIgnum * alpha)));
+ }
+ nb = 1 + (long)na; /* nb-1 <= alpha < nb */
+ alpha -= (gnm_float)(nb-1);
+#ifdef MATHLIB_STANDALONE
+ bj = (gnm_float *) calloc(nb, sizeof(gnm_float));
+ if (!bj) MATHLIB_ERROR("%s", _("bessel_j allocation error"));
+#else
+ vmax = vmaxget();
+ bj = (gnm_float *) R_alloc((size_t) nb, sizeof(gnm_float));
+#endif
+ J_bessel(&x, &alpha, &nb, bj, &ncalc);
+ if(ncalc != nb) {/* error input */
+ if(ncalc < 0)
+ MATHLIB_WARNING4(_("bessel_j(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g ".
Arg. out of range?\n"),
+ x, ncalc, nb, alpha);
+ else
+ MATHLIB_WARNING2(_("bessel_j(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in result\n"),
+ x, alpha+(gnm_float)nb-1);
+ }
+ x = bj[nb-1];
+#ifdef MATHLIB_STANDALONE
+ free(bj);
+#else
+ vmaxset(vmax);
+#endif
+ return x;
+}
+
+/* modified version of bessel_j that accepts a work array instead of
+ allocating one. */
+static gnm_float bessel_j_ex(gnm_float x, gnm_float alpha, gnm_float *bj)
+{
+ long nb, ncalc;
+ gnm_float na;
+
+#ifdef IEEE_754
+ /* NaNs propagated correctly */
+ if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
+#endif
+ if (x < 0) {
+ ML_ERROR(ME_RANGE);
+ return gnm_nan;
+ }
+ na = gnm_floor(alpha);
+ if (alpha < 0) {
+ /* Using Abramowitz & Stegun 9.1.2
+ * this may not be quite optimal (CPU and accuracy wise) */
+ return(bessel_j_ex(x, -alpha, bj) * gnm_cos(M_PIgnum * alpha) +
+ ((alpha == na) ? 0 :
+ bessel_y_ex(x, -alpha, bj) * gnm_sin(M_PIgnum * alpha)));
+ }
+ nb = 1 + (long)na; /* nb-1 <= alpha < nb */
+ alpha -= (gnm_float)(nb-1);
+ J_bessel(&x, &alpha, &nb, bj, &ncalc);
+ if(ncalc != nb) {/* error input */
+ if(ncalc < 0)
+ MATHLIB_WARNING4(_("bessel_j(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g ".
Arg. out of range?\n"),
+ x, ncalc, nb, alpha);
+ else
+ MATHLIB_WARNING2(_("bessel_j(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in result\n"),
+ x, alpha+(gnm_float)nb-1);
+ }
+ x = bj[nb-1];
+ return x;
+}
+
+static void J_bessel(gnm_float *x, gnm_float *alpha, long *nb,
+ gnm_float *b, long *ncalc)
+{
+/*
+ Calculates Bessel functions J_{n+alpha} (x)
+ for non-negative argument x, and non-negative order n+alpha, n = 0,1,..,nb-1.
+
+ Explanation of variables in the calling sequence.
+
+ X - Non-negative argument for which J's are to be calculated.
+ ALPHA - Fractional part of order for which
+ J's are to be calculated. 0 <= ALPHA < 1.
+ NB - Number of functions to be calculated, NB >= 1.
+ The first function calculated is of order ALPHA, and the
+ last is of order (NB - 1 + ALPHA).
+ B - Output vector of length NB. If RJBESL
+ terminates normally (NCALC=NB), the vector B contains the
+ functions J/ALPHA/(X) through J/NB-1+ALPHA/(X).
+ NCALC - Output variable indicating possible errors.
+ Before using the vector B, the user should check that
+ NCALC=NB, i.e., all orders have been calculated to
+ the desired accuracy. See the following
+
+ ****************************************************************
+
+ Error return codes
+
+ In case of an error, NCALC != NB, and not all J's are
+ calculated to the desired accuracy.
+
+ NCALC < 0: An argument is out of range. For example,
+ NBES <= 0, ALPHA < 0 or > 1, or X is too large.
+ In this case, b[1] is set to zero, the remainder of the
+ B-vector is not calculated, and NCALC is set to
+ MIN(NB,0)-1 so that NCALC != NB.
+
+ NB > NCALC > 0: Not all requested function values could
+ be calculated accurately. This usually occurs because NB is
+ much larger than ABS(X). In this case, b[N] is calculated
+ to the desired accuracy for N <= NCALC, but precision
+ is lost for NCALC < N <= NB. If b[N] does not vanish
+ for N > NCALC (because it is too small to be represented),
+ and b[N]/b[NCALC] = 10^(-K), then only the first NSIG - K
+ significant figures of b[N] can be trusted.
+
+
+ Acknowledgement
+
+ This program is based on a program written by David J. Sookne
+ (2) that computes values of the Bessel functions J or I of float
+ argument and long order. Modifications include the restriction
+ of the computation to the J Bessel function of non-negative float
+ argument, the extension of the computation to arbitrary positive
+ order, and the elimination of most underflow.
+
+ References:
+
+ Olver, F.W.J., and Sookne, D.J. (1972)
+ "A Note on Backward Recurrence Algorithms";
+ Math. Comp. 26, 941-947.
+
+ Sookne, D.J. (1973)
+ "Bessel Functions of Real Argument and Integer Order";
+ NBS Jour. of Res. B. 77B, 125-132.
+
+ Latest modification: March 19, 1990
+
+ Author: W. J. Cody
+ Applied Mathematics Division
+ Argonne National Laboratory
+ Argonne, IL 60439
+ *******************************************************************
+ */
+
+/* ---------------------------------------------------------------------
+ Mathematical constants
+
+ PI2 = 2 / PI
+ TWOPI1 = first few significant digits of 2 * PI
+ TWOPI2 = (2*PI - TWOPI1) to working precision, i.e.,
+ TWOPI1 + TWOPI2 = 2 * PI to extra precision.
+ --------------------------------------------------------------------- */
+ const static gnm_float pi2 = GNM_const(.636619772367581343075535);
+ const static gnm_float twopi1 = GNM_const(6.28125);
+ const static gnm_float twopi2 = GNM_const(.001935307179586476925286767);
+
+/*---------------------------------------------------------------------
+ * Factorial(N)
+ *--------------------------------------------------------------------- */
+/* removed array fact */
+
+ /* Local variables */
+ long nend, intx, nbmx, i, j, k, l, m, n, nstart;
+
+ gnm_float nu, twonu, capp, capq, pold, vcos, test, vsin;
+ gnm_float p, s, t, z, alpem, halfx, aa, bb, cc, psave, plast;
+ gnm_float tover, t1, alp2em, em, en, xc, xk, xm, psavel, gnu, xin, sum;
+
+
+ /* Parameter adjustment */
+ --b;
+
+ nu = *alpha;
+ twonu = nu + nu;
+
+ /*-------------------------------------------------------------------
+ Check for out of range arguments.
+ -------------------------------------------------------------------*/
+ if (*nb > 0 && *x >= 0. && 0. <= nu && nu < 1.) {
+
+ *ncalc = *nb;
+ if(*x > xlrg_BESS_IJ) {
+ ML_ERROR(ME_RANGE);
+ /* indeed, the limit is 0,
+ * but the cutoff happens too early */
+ for(i=1; i <= *nb; i++)
+ b[i] = 0.; /*was gnm_pinf (really nonsense) */
+ return;
+ }
+ intx = (long) (*x);
+ /* Initialize result array to zero. */
+ for (i = 1; i <= *nb; ++i)
+ b[i] = 0.;
+
+ /*===================================================================
+ Branch into 3 cases :
+ 1) use 2-term ascending series for small X
+ 2) use asymptotic form for large X when NB is not too large
+ 3) use recursion otherwise
+ ===================================================================*/
+
+ if (*x < rtnsig_BESS) {
+ /* ---------------------------------------------------------------
+ Two-term ascending series for small X.
+ --------------------------------------------------------------- */
+ alpem = 1. + nu;
+
+ halfx = (*x > enmten_BESS) ? .5 * *x : 0.;
+ aa = (nu != 0.) ? gnm_pow(halfx, nu) / (nu * gnm_gamma(nu)) : 1.;
+ bb = (*x + 1. > 1.)? -halfx * halfx : 0.;
+ b[1] = aa + aa * bb / alpem;
+ if (*x != 0. && b[1] == 0.)
+ *ncalc = 0;
+
+ if (*nb != 1) {
+ if (*x <= 0.) {
+ for (n = 2; n <= *nb; ++n)
+ b[n] = 0.;
+ }
+ else {
+ /* ----------------------------------------------
+ Calculate higher order functions.
+ ---------------------------------------------- */
+ if (bb == 0.)
+ tover = (enmten_BESS + enmten_BESS) / *x;
+ else
+ tover = enmten_BESS / bb;
+ cc = halfx;
+ for (n = 2; n <= *nb; ++n) {
+ aa /= alpem;
+ alpem += 1.;
+ aa *= cc;
+ if (aa <= tover * alpem)
+ aa = 0.;
+
+ b[n] = aa + aa * bb / alpem;
+ if (b[n] == 0. && *ncalc > n)
+ *ncalc = n - 1;
+ }
+ }
+ }
+ } else if (*x > 25. && *nb <= intx + 1) {
+ /* ------------------------------------------------------------
+ Asymptotic series for X > 25 (and not too large nb)
+ ------------------------------------------------------------ */
+ xc = gnm_sqrt(pi2 / *x);
+ xin = 1 / (64 * *x * *x);
+ if (*x >= 130.) m = 4;
+ else if (*x >= 35.) m = 8;
+ else m = 11;
+ xm = 4. * (gnm_float) m;
+ /* ------------------------------------------------
+ Argument reduction for SIN and COS routines.
+ ------------------------------------------------ */
+ t = gnm_trunc(*x / (twopi1 + twopi2) + .5);
+ z = (*x - t * twopi1) - t * twopi2 - (nu + .5) / pi2;
+ vsin = gnm_sin(z);
+ vcos = gnm_cos(z);
+ gnu = twonu;
+ for (i = 1; i <= 2; ++i) {
+ s = (xm - 1. - gnu) * (xm - 1. + gnu) * xin * .5;
+ t = (gnu - (xm - 3.)) * (gnu + (xm - 3.));
+ t1= (gnu - (xm + 1.)) * (gnu + (xm + 1.));
+ k = m + m;
+ capp = s * t / fact(k);
+ capq = s * t1/ fact(k + 1);
+ xk = xm;
+ for (; k >= 4; k -= 2) {/* k + 2(j-2) == 2m */
+ xk -= 4.;
+ s = (xk - 1. - gnu) * (xk - 1. + gnu);
+ t1 = t;
+ t = (gnu - (xk - 3.)) * (gnu + (xk - 3.));
+ capp = (capp + 1. / fact(k - 2)) * s * t * xin;
+ capq = (capq + 1. / fact(k - 1)) * s * t1 * xin;
+
+ }
+ capp += 1.;
+ capq = (capq + 1.) * (gnu * gnu - 1.) * (.125 / *x);
+ b[i] = xc * (capp * vcos - capq * vsin);
+ if (*nb == 1)
+ return;
+
+ /* vsin <--> vcos */ t = vsin; vsin = -vcos; vcos = t;
+ gnu += 2.;
+ }
+ /* -----------------------------------------------
+ If NB > 2, compute J(X,ORDER+I) for I = 2, NB-1
+ ----------------------------------------------- */
+ if (*nb > 2)
+ for (gnu = twonu + 2., j = 3; j <= *nb; j++, gnu += 2.)
+ b[j] = gnu * b[j - 1] / *x - b[j - 2];
+ }
+ else {
+ /* rtnsig_BESS <= x && ( x <= 25 || intx+1 < *nb ) :
+ --------------------------------------------------------
+ Use recurrence to generate results.
+ First initialize the calculation of P*S.
+ -------------------------------------------------------- */
+ nbmx = *nb - intx;
+ n = intx + 1;
+ en = (gnm_float)(n + n) + twonu;
+ plast = 1.;
+ p = en / *x;
+ /* ---------------------------------------------------
+ Calculate general significance test.
+ --------------------------------------------------- */
+ test = ensig_BESS + ensig_BESS;
+ if (nbmx >= 3) {
+ /* ------------------------------------------------------------
+ Calculate P*S until N = NB-1. Check for possible overflow.
+ ---------------------------------------------------------- */
+ tover = enten_BESS / ensig_BESS;
+ nstart = intx + 2;
+ nend = *nb - 1;
+ en = (gnm_float) (nstart + nstart) - 2. + twonu;
+ for (k = nstart; k <= nend; ++k) {
+ n = k;
+ en += 2.;
+ pold = plast;
+ plast = p;
+ p = en * plast / *x - pold;
+ if (p > tover) {
+ /* -------------------------------------------
+ To avoid overflow, divide P*S by TOVER.
+ Calculate P*S until ABS(P) > 1.
+ -------------------------------------------*/
+ tover = enten_BESS;
+ p /= tover;
+ plast /= tover;
+ psave = p;
+ psavel = plast;
+ nstart = n + 1;
+ do {
+ ++n;
+ en += 2.;
+ pold = plast;
+ plast = p;
+ p = en * plast / *x - pold;
+ } while (p <= 1.);
+
+ bb = en / *x;
+ /* -----------------------------------------------
+ Calculate backward test and find NCALC,
+ the highest N such that the test is passed.
+ ----------------------------------------------- */
+ test = pold * plast * (.5 - .5 / (bb * bb));
+ test /= ensig_BESS;
+ p = plast * tover;
+ --n;
+ en -= 2.;
+ nend = min0(*nb,n);
+ for (l = nstart; l <= nend; ++l) {
+ pold = psavel;
+ psavel = psave;
+ psave = en * psavel / *x - pold;
+ if (psave * psavel > test) {
+ *ncalc = l - 1;
+ goto L190;
+ }
+ }
+ *ncalc = nend;
+ goto L190;
+ }
+ }
+ n = nend;
+ en = (gnm_float) (n + n) + twonu;
+ /* -----------------------------------------------------
+ Calculate special significance test for NBMX > 2.
+ -----------------------------------------------------*/
+ test = fmax2(test, gnm_sqrt(plast * ensig_BESS) * gnm_sqrt(p + p));
+ }
+ /* ------------------------------------------------
+ Calculate P*S until significance test passes. */
+ do {
+ ++n;
+ en += 2.;
+ pold = plast;
+ plast = p;
+ p = en * plast / *x - pold;
+ } while (p < test);
+
+L190:
+ /*---------------------------------------------------------------
+ Initialize the backward recursion and the normalization sum.
+ --------------------------------------------------------------- */
+ ++n;
+ en += 2.;
+ bb = 0.;
+ aa = 1. / p;
+ m = n / 2;
+ em = (gnm_float)m;
+ m = (n << 1) - (m << 2);/* = 2 n - 4 (n/2)
+ = 0 for even, 2 for odd n */
+ if (m == 0)
+ sum = 0.;
+ else {
+ alpem = em - 1. + nu;
+ alp2em = em + em + nu;
+ sum = aa * alpem * alp2em / em;
+ }
+ nend = n - *nb;
+ /* if (nend > 0) */
+ /* --------------------------------------------------------
+ Recur backward via difference equation, calculating
+ (but not storing) b[N], until N = NB.
+ -------------------------------------------------------- */
+ for (l = 1; l <= nend; ++l) {
+ --n;
+ en -= 2.;
+ cc = bb;
+ bb = aa;
+ aa = en * bb / *x - cc;
+ m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
+ if (m != 0) {
+ em -= 1.;
+ alp2em = em + em + nu;
+ if (n == 1)
+ break;
+
+ alpem = em - 1. + nu;
+ if (alpem == 0.)
+ alpem = 1.;
+ sum = (sum + aa * alp2em) * alpem / em;
+ }
+ }
+ /*--------------------------------------------------
+ Store b[NB].
+ --------------------------------------------------*/
+ b[n] = aa;
+ if (nend >= 0) {
+ if (*nb <= 1) {
+ if (nu + 1. == 1.)
+ alp2em = 1.;
+ else
+ alp2em = nu;
+ sum += b[1] * alp2em;
+ goto L250;
+ }
+ else {/*-- nb >= 2 : ---------------------------
+ Calculate and store b[NB-1].
+ ----------------------------------------*/
+ --n;
+ en -= 2.;
+ b[n] = en * aa / *x - bb;
+ if (n == 1)
+ goto L240;
+
+ m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
+ if (m != 0) {
+ em -= 1.;
+ alp2em = em + em + nu;
+ alpem = em - 1. + nu;
+ if (alpem == 0.)
+ alpem = 1.;
+ sum = (sum + b[n] * alp2em) * alpem / em;
+ }
+ }
+ }
+
+ /* if (n - 2 != 0) */
+ /* --------------------------------------------------------
+ Calculate via difference equation and store b[N],
+ until N = 2.
+ -------------------------------------------------------- */
+ for (n = n-1; n >= 2; n--) {
+ en -= 2.;
+ b[n] = en * b[n + 1] / *x - b[n + 2];
+ m = m ? 0 : 2; /* m = 2 - m failed on gcc4-20041019 */
+ if (m != 0) {
+ em -= 1.;
+ alp2em = em + em + nu;
+ alpem = em - 1. + nu;
+ if (alpem == 0.)
+ alpem = 1.;
+ sum = (sum + b[n] * alp2em) * alpem / em;
+ }
+ }
+ /* ---------------------------------------
+ Calculate b[1].
+ -----------------------------------------*/
+ b[1] = 2. * (nu + 1.) * b[2] / *x - b[3];
+
+L240:
+ em -= 1.;
+ alp2em = em + em + nu;
+ if (alp2em == 0.)
+ alp2em = 1.;
+ sum += b[1] * alp2em;
+
+L250:
+ /* ---------------------------------------------------
+ Normalize. Divide all b[N] by sum.
+ ---------------------------------------------------*/
+/* if (nu + 1. != 1.) poor test */
+ if(gnm_abs(nu) > 1e-15)
+ sum *= (gnm_gamma(nu) * gnm_pow(.5* *x, -nu));
+
+ aa = enmten_BESS;
+ if (sum > 1.)
+ aa *= sum;
+ for (n = 1; n <= *nb; ++n) {
+ if (gnm_abs(b[n]) < aa)
+ b[n] = 0.;
+ else
+ b[n] /= sum;
+ }
+ }
+
+ }
+ else {
+ /* Error return -- X, NB, or ALPHA is out of range : */
+ b[1] = 0.;
+ *ncalc = min0(*nb,0) - 1;
+ }
+}
+/* Cleaning up done by tools/import-R: */
+#undef min0
+
+/* ------------------------------------------------------------------------ */
/* Imported src/nmath/bessel_k.c from R. */
/*
* Mathlib : A C Library of Special Functions
@@ -4831,6 +5429,511 @@ L420:
}
/* ------------------------------------------------------------------------ */
+/* Imported src/nmath/bessel_y.c from R. */
+/*
+ * Mathlib : A C Library of Special Functions
+ * Copyright (C) 1998-2012 Ross Ihaka and 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
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, a copy is available at
+ * http://www.r-project.org/Licenses/
+ */
+
+/* DESCRIPTION --> see below */
+
+
+/* From http://www.netlib.org/specfun/rybesl Fortran translated by f2c,...
+ * ------------------------------=#---- Martin Maechler, ETH Zurich
+ */
+
+#ifndef MATHLIB_STANDALONE
+#endif
+
+#define min0(x, y) (((x) <= (y)) ? (x) : (y))
+
+static void Y_bessel(gnm_float *x, gnm_float *alpha, long *nb,
+ gnm_float *by, long *ncalc);
+
+gnm_float bessel_y(gnm_float x, gnm_float alpha)
+{
+ long nb, ncalc;
+ gnm_float na, *by;
+#ifndef MATHLIB_STANDALONE
+ const void *vmax;
+#endif
+
+#ifdef IEEE_754
+ /* NaNs propagated correctly */
+ if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
+#endif
+ if (x < 0) {
+ ML_ERROR(ME_RANGE);
+ return gnm_nan;
+ }
+ na = gnm_floor(alpha);
+ if (alpha < 0) {
+ /* Using Abramowitz & Stegun 9.1.2
+ * this may not be quite optimal (CPU and accuracy wise) */
+ return(bessel_y(x, -alpha) * gnm_cos(M_PIgnum * alpha) -
+ ((alpha == na) ? 0 :
+ bessel_j(x, -alpha) * gnm_sin(M_PIgnum * alpha)));
+ }
+ nb = 1+ (long)na;/* nb-1 <= alpha < nb */
+ alpha -= (gnm_float)(nb-1);
+#ifdef MATHLIB_STANDALONE
+ by = (gnm_float *) calloc(nb, sizeof(gnm_float));
+ if (!by) MATHLIB_ERROR("%s", _("bessel_y allocation error"));
+#else
+ vmax = vmaxget();
+ by = (gnm_float *) R_alloc((size_t) nb, sizeof(gnm_float));
+#endif
+ Y_bessel(&x, &alpha, &nb, by, &ncalc);
+ if(ncalc != nb) {/* error input */
+ if(ncalc == -1)
+ return gnm_pinf;
+ else if(ncalc < -1)
+ MATHLIB_WARNING4(_("bessel_y(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g
". Arg. out of range?\n"),
+ x, ncalc, nb, alpha);
+ else /* ncalc >= 0 */
+ MATHLIB_WARNING2(_("bessel_y(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in
result\n"),
+ x, alpha+(gnm_float)nb-1);
+ }
+ x = by[nb-1];
+#ifdef MATHLIB_STANDALONE
+ free(by);
+#else
+ vmaxset(vmax);
+#endif
+ return x;
+}
+
+/* modified version of bessel_y that accepts a work array instead of
+ allocating one. */
+static gnm_float bessel_y_ex(gnm_float x, gnm_float alpha, gnm_float *by)
+{
+ long nb, ncalc;
+ gnm_float na;
+
+#ifdef IEEE_754
+ /* NaNs propagated correctly */
+ if (gnm_isnan(x) || gnm_isnan(alpha)) return x + alpha;
+#endif
+ if (x < 0) {
+ ML_ERROR(ME_RANGE);
+ return gnm_nan;
+ }
+ na = gnm_floor(alpha);
+ if (alpha < 0) {
+ /* Using Abramowitz & Stegun 9.1.2
+ * this may not be quite optimal (CPU and accuracy wise) */
+ return(bessel_y_ex(x, -alpha, by) * gnm_cos(M_PIgnum * alpha) -
+ ((alpha == na) ? 0 :
+ bessel_j_ex(x, -alpha, by) * gnm_sin(M_PIgnum * alpha)));
+ }
+ nb = 1+ (long)na;/* nb-1 <= alpha < nb */
+ alpha -= (gnm_float)(nb-1);
+ Y_bessel(&x, &alpha, &nb, by, &ncalc);
+ if(ncalc != nb) {/* error input */
+ if(ncalc == -1)
+ return gnm_pinf;
+ else if(ncalc < -1)
+ MATHLIB_WARNING4(_("bessel_y(%" GNM_FORMAT_g "): ncalc (=%ld) != nb (=%ld); alpha=%" GNM_FORMAT_g
". Arg. out of range?\n"),
+ x, ncalc, nb, alpha);
+ else /* ncalc >= 0 */
+ MATHLIB_WARNING2(_("bessel_y(%" GNM_FORMAT_g ",nu=%" GNM_FORMAT_g "): precision lost in
result\n"),
+ x, alpha+(gnm_float)nb-1);
+ }
+ x = by[nb-1];
+ return x;
+}
+
+static void Y_bessel(gnm_float *x, gnm_float *alpha, long *nb,
+ gnm_float *by, long *ncalc)
+{
+/* ----------------------------------------------------------------------
+
+ This routine calculates Bessel functions Y_(N+ALPHA) (X)
+ for non-negative argument X, and non-negative order N+ALPHA.
+
+
+ Explanation of variables in the calling sequence
+
+ X - Non-negative argument for which
+ Y's are to be calculated.
+ ALPHA - Fractional part of order for which
+ Y's are to be calculated. 0 <= ALPHA < 1.0.
+ NB - Number of functions to be calculated, NB > 0.
+ The first function calculated is of order ALPHA, and the
+ last is of order (NB - 1 + ALPHA).
+ BY - Output vector of length NB. If the
+ routine terminates normally (NCALC=NB), the vector BY
+ contains the functions Y(ALPHA,X), ... , Y(NB-1+ALPHA,X),
+ If (0 < NCALC < NB), BY(I) contains correct function
+ values for I <= NCALC, and contains the ratios
+ Y(ALPHA+I-1,X)/Y(ALPHA+I-2,X) for the rest of the array.
+ NCALC - Output variable indicating possible errors.
+ Before using the vector BY, the user should check that
+ NCALC=NB, i.e., all orders have been calculated to
+ the desired accuracy. See error returns below.
+
+
+ *******************************************************************
+
+ Error returns
+
+ In case of an error, NCALC != NB, and not all Y's are
+ calculated to the desired accuracy.
+
+ NCALC < -1: An argument is out of range. For example,
+ NB <= 0, IZE is not 1 or 2, or IZE=1 and ABS(X) >=
+ XMAX. In this case, BY[0] = 0.0, the remainder of the
+ BY-vector is not calculated, and NCALC is set to
+ MIN0(NB,0)-2 so that NCALC != NB.
+ NCALC = -1: Y(ALPHA,X) >= XINF. The requested function
+ values are set to 0.0.
+ 1 < NCALC < NB: Not all requested function values could
+ be calculated accurately. BY(I) contains correct function
+ values for I <= NCALC, and and the remaining NB-NCALC
+ array elements contain 0.0.
+
+
+ Intrinsic functions required are:
+
+ DBLE, EXP, INT, MAX, MIN, REAL, SQRT
+
+
+ Acknowledgement
+
+ This program draws heavily on Temme's Algol program for Y(a,x)
+ and Y(a+1,x) and on Campbell's programs for Y_nu(x). Temme's
+ scheme is used for x < THRESH, and Campbell's scheme is used
+ in the asymptotic region. Segments of code from both sources
+ have been translated into Fortran 77, merged, and heavily modified.
+ Modifications include parameterization of machine dependencies,
+ use of a new approximation for ln(gamma(x)), and built-in
+ protection against over/underflow.
+
+ References: "Bessel functions J_nu(x) and Y_nu(x) of float
+ order and float argument," Campbell, J. B.,
+ Comp. Phy. Comm. 18, 1979, pp. 133-142.
+
+ "On the numerical evaluation of the ordinary
+ Bessel function of the second kind," Temme,
+ N. M., J. Comput. Phys. 21, 1976, pp. 343-350.
+
+ Latest modification: March 19, 1990
+
+ Modified by: W. J. Cody
+ Applied Mathematics Division
+ Argonne National Laboratory
+ Argonne, IL 60439
+ ----------------------------------------------------------------------*/
+
+
+/* ----------------------------------------------------------------------
+ Mathematical constants
+ FIVPI = 5*PI
+ PIM5 = 5*PI - 15
+ ----------------------------------------------------------------------*/
+ const static gnm_float fivpi = GNM_const(15.707963267948966192);
+ const static gnm_float pim5 = GNM_const(.70796326794896619231);
+
+ /*----------------------------------------------------------------------
+ Coefficients for Chebyshev polynomial expansion of
+ 1/gamma(1-x), abs(x) <= .5
+ ----------------------------------------------------------------------*/
+ const static gnm_float ch[21] = { GNM_const(-6.7735241822398840964e-24),
+ GNM_const(-6.1455180116049879894e-23),GNM_const(2.9017595056104745456e-21),
+ GNM_const(1.3639417919073099464e-19),GNM_const(2.3826220476859635824e-18),
+ GNM_const(-9.0642907957550702534e-18),GNM_const(-1.4943667065169001769e-15),
+ GNM_const(-3.3919078305362211264e-14),GNM_const(-1.7023776642512729175e-13),
+ GNM_const(9.1609750938768647911e-12),GNM_const(2.4230957900482704055e-10),
+ GNM_const(1.7451364971382984243e-9),GNM_const(-3.3126119768180852711e-8),
+ GNM_const(-8.6592079961391259661e-7),GNM_const(-4.9717367041957398581e-6),
+ GNM_const(7.6309597585908126618e-5),GNM_const(.0012719271366545622927),
+ GNM_const(.0017063050710955562222),GNM_const(-.07685284084478667369),
+ GNM_const(-.28387654227602353814),GNM_const(.92187029365045265648) };
+
+ /* Local variables */
+ long i, k, na;
+
+ gnm_float alfa, div, ddiv, even, gamma, term, cosmu, sinmu,
+ b, c, d, e, f, g, h, p, q, r, s, d1, d2, q0, pa,pa1, qa,qa1,
+ en, en1, nu, ex, ya,ya1, twobyx, den, odd, aye, dmu, x2, xna;
+
+ en1 = ya = ya1 = 0; /* -Wall */
+
+ ex = *x;
+ nu = *alpha;
+ if (*nb > 0 && 0. <= nu && nu < 1.) {
+ if(ex < GNM_MIN || ex > xlrg_BESS_Y) {
+ /* Warning is not really appropriate, give
+ * proper limit:
+ * ML_ERROR(ME_RANGE); */
+ *ncalc = *nb;
+ if(ex > xlrg_BESS_Y) by[0]= 0.; /*was gnm_pinf */
+ else if(ex < GNM_MIN) by[0]=gnm_ninf;
+ for(i=0; i < *nb; i++)
+ by[i] = by[0];
+ return;
+ }
+ xna = gnm_trunc(nu + .5);
+ na = (long) xna;
+ if (na == 1) {/* <==> .5 <= *alpha < 1 <==> -5. <= nu < 0 */
+ nu -= xna;
+ }
+ if (nu == -.5) {
+ p = M_SQRT_2dPI / gnm_sqrt(ex);
+ ya = p * gnm_sin(ex);
+ ya1 = -p * gnm_cos(ex);
+ } else if (ex < 3.) {
+ /* -------------------------------------------------------------
+ Use Temme's scheme for small X
+ ------------------------------------------------------------- */
+ b = ex * .5;
+ d = -gnm_log(b);
+ f = nu * d;
+ e = gnm_pow(b, -nu);
+ if (gnm_abs(nu) < M_eps_sinc)
+ c = M_1_PI;
+ else
+ c = nu / gnm_sin(nu * M_PIgnum);
+
+ /* ------------------------------------------------------------
+ Computation of gnm_sinh(f)/f
+ ------------------------------------------------------------ */
+ if (gnm_abs(f) < 1.) {
+ x2 = f * f;
+ en = 19.;
+ s = 1.;
+ for (i = 1; i <= 9; ++i) {
+ s = s * x2 / en / (en - 1.) + 1.;
+ en -= 2.;
+ }
+ } else {
+ s = (e - GNM_const(1.) / e) * .5 / f;
+ }
+ /* --------------------------------------------------------
+ Computation of 1/gamma(1-a) using Chebyshev polynomials */
+ x2 = nu * nu * 8.;
+ aye = ch[0];
+ even = 0.;
+ alfa = ch[1];
+ odd = 0.;
+ for (i = 3; i <= 19; i += 2) {
+ even = -(aye + aye + even);
+ aye = -even * x2 - aye + ch[i - 1];
+ odd = -(alfa + alfa + odd);
+ alfa = -odd * x2 - alfa + ch[i];
+ }
+ even = (even * .5 + aye) * x2 - aye + ch[20];
+ odd = (odd + alfa) * 2.;
+ gamma = odd * nu + even;
+ /* End of computation of 1/gamma(1-a)
+ ----------------------------------------------------------- */
+ g = e * gamma;
+ e = (e + GNM_const(1.) / e) * .5;
+ f = 2. * c * (odd * e + even * s * d);
+ e = nu * nu;
+ p = g * c;
+ q = M_1_PI / g;
+ c = nu * M_PI_2gnum;
+ if (gnm_abs(c) < M_eps_sinc)
+ r = 1.;
+ else
+ r = gnm_sin(c) / c;
+
+ r = M_PIgnum * c * r * r;
+ c = 1.;
+ d = -b * b;
+ h = 0.;
+ ya = f + r * q;
+ ya1 = p;
+ en = 1.;
+
+ while (gnm_abs(g / (1. + gnm_abs(ya))) +
+ gnm_abs(h / (1. + gnm_abs(ya1))) > GNM_EPSILON) {
+ f = (f * en + p + q) / (en * en - e);
+ c *= (d / en);
+ p /= en - nu;
+ q /= en + nu;
+ g = c * (f + r * q);
+ h = c * p - en * g;
+ ya += g;
+ ya1+= h;
+ en += 1.;
+ }
+ ya = -ya;
+ ya1 = -ya1 / b;
+ } else if (ex < thresh_BESS_Y) {
+ /* --------------------------------------------------------------
+ Use Temme's scheme for moderate X : 3 <= x < 16
+ -------------------------------------------------------------- */
+ c = (.5 - nu) * (.5 + nu);
+ b = ex + ex;
+ e = ex * M_1_PI * gnm_cos(nu * M_PIgnum) / GNM_EPSILON;
+ e *= e;
+ p = 1.;
+ q = -ex;
+ r = 1. + ex * ex;
+ s = r;
+ en = 2.;
+ while (r * en * en < e) {
+ en1 = en + 1.;
+ d = (en - 1. + c / en) / s;
+ p = (en + en - p * d) / en1;
+ q = (-b + q * d) / en1;
+ s = p * p + q * q;
+ r *= s;
+ en = en1;
+ }
+ f = p / s;
+ p = f;
+ g = -q / s;
+ q = g;
+L220:
+ en -= 1.;
+ if (en > 0.) {
+ r = en1 * (2. - p) - 2.;
+ s = b + en1 * q;
+ d = (en - 1. + c / en) / (r * r + s * s);
+ p = d * r;
+ q = d * s;
+ e = f + 1.;
+ f = p * e - g * q;
+ g = q * e + p * g;
+ en1 = en;
+ goto L220;
+ }
+ f = 1. + f;
+ d = f * f + g * g;
+ pa = f / d;
+ qa = -g / d;
+ d = nu + .5 - p;
+ q += ex;
+ pa1 = (pa * q - qa * d) / ex;
+ qa1 = (qa * q + pa * d) / ex;
+ b = ex - M_PI_2gnum * (nu + .5);
+ c = gnm_cos(b);
+ s = gnm_sin(b);
+ d = M_SQRT_2dPI / gnm_sqrt(ex);
+ ya = d * (pa * s + qa * c);
+ ya1 = d * (qa1 * s - pa1 * c);
+ } else { /* x > thresh_BESS_Y */
+ /* ----------------------------------------------------------
+ Use Campbell's asymptotic scheme.
+ ---------------------------------------------------------- */
+ na = 0;
+ d1 = gnm_trunc(ex / fivpi);
+ i = (long) d1;
+ dmu = ex - 15. * d1 - d1 * pim5 - (*alpha + .5) * M_PI_2gnum;
+ if (i - (i / 2 << 1) == 0) {
+ cosmu = gnm_cos(dmu);
+ sinmu = gnm_sin(dmu);
+ } else {
+ cosmu = -gnm_cos(dmu);
+ sinmu = -gnm_sin(dmu);
+ }
+ ddiv = 8. * ex;
+ dmu = *alpha;
+ den = gnm_sqrt(ex);
+ for (k = 1; k <= 2; ++k) {
+ p = cosmu;
+ cosmu = sinmu;
+ sinmu = -p;
+ d1 = (2. * dmu - 1.) * (2. * dmu + 1.);
+ d2 = 0.;
+ div = ddiv;
+ p = 0.;
+ q = 0.;
+ q0 = d1 / div;
+ term = q0;
+ for (i = 2; i <= 20; ++i) {
+ d2 += 8.;
+ d1 -= d2;
+ div += ddiv;
+ term = -term * d1 / div;
+ p += term;
+ d2 += 8.;
+ d1 -= d2;
+ div += ddiv;
+ term *= (d1 / div);
+ q += term;
+ if (gnm_abs(term) <= GNM_EPSILON) {
+ break;
+ }
+ }
+ p += 1.;
+ q += q0;
+ if (k == 1)
+ ya = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
+ else
+ ya1 = M_SQRT_2dPI * (p * cosmu - q * sinmu) / den;
+ dmu += 1.;
+ }
+ }
+ if (na == 1) {
+ h = 2. * (nu + 1.) / ex;
+ if (h > 1.) {
+ if (gnm_abs(ya1) > GNM_MAX / h) {
+ h = 0.;
+ ya = 0.;
+ }
+ }
+ h = h * ya1 - ya;
+ ya = ya1;
+ ya1 = h;
+ }
+
+ /* ---------------------------------------------------------------
+ Now have first one or two Y's
+ --------------------------------------------------------------- */
+ by[0] = ya;
+ *ncalc = 1;
+ if(*nb > 1) {
+ by[1] = ya1;
+ if (ya1 != 0.) {
+ aye = 1. + *alpha;
+ twobyx = GNM_const(2.) / ex;
+ *ncalc = 2;
+ for (i = 2; i < *nb; ++i) {
+ if (twobyx < 1.) {
+ if (gnm_abs(by[i - 1]) * twobyx >= GNM_MAX / aye)
+ goto L450;
+ } else {
+ if (gnm_abs(by[i - 1]) >= GNM_MAX / aye / twobyx)
+ goto L450;
+ }
+ by[i] = twobyx * aye * by[i - 1] - by[i - 2];
+ aye += 1.;
+ ++(*ncalc);
+ }
+ }
+ }
+L450:
+ for (i = *ncalc; i < *nb; ++i)
+ by[i] = gnm_ninf;/* was 0 */
+
+ } else {
+ by[0] = 0.;
+ *ncalc = min0(*nb,0) - 1;
+ }
+}
+
+/* Cleaning up done by tools/import-R: */
+#undef min0
+
+/* ------------------------------------------------------------------------ */
/* Imported src/nmath/dlnorm.c from R. */
/*
* Mathlib : A C Library of Special Functions
diff --git a/src/mathfunc.h b/src/mathfunc.h
index 2c15b61..358deee 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -37,6 +37,7 @@ gnm_float logspace_sub (gnm_float logx, gnm_float logy);
gnm_float stirlerr(gnm_float n);
gnm_float gnm_owent (gnm_float h, gnm_float a);
gnm_float pochhammer (gnm_float x, gnm_float n, gboolean give_log);
+gnm_float gnm_gamma (gnm_float x);
gnm_float gnm_cot (gnm_float x);
gnm_float gnm_acot (gnm_float x);
@@ -47,7 +48,9 @@ gnm_float beta (gnm_float a, gnm_float b);
gnm_float lbeta3 (gnm_float a, gnm_float b, int *sign);
gnm_float bessel_i (gnm_float x, gnm_float alpha, gnm_float expo);
+gnm_float bessel_j (gnm_float x, gnm_float alpha);
gnm_float bessel_k (gnm_float x, gnm_float alpha, gnm_float expo);
+gnm_float bessel_y (gnm_float x, gnm_float alpha);
/* "d": density. */
/* "p": distribution function. */
diff --git a/src/outoflinedocs.c b/src/outoflinedocs.c
index 847d564..841e98f 100644
--- a/src/outoflinedocs.c
+++ b/src/outoflinedocs.c
@@ -63,6 +63,14 @@
*/
/**
+ * gamma:
+ * @x: a number
+ *
+ * Returns: gamma(@x) for for positive or non-integer @x.
+ */
+
+
+/**
* beta:
* @a: a number
* @b: a number
diff --git a/tools/import-R b/tools/import-R
index 7b8b8ee..7c050fe 100644
--- a/tools/import-R
+++ b/tools/import-R
@@ -50,7 +50,9 @@ my @files =
"pcauchy.c",
"bessel.h",
"bessel_i.c",
+ "bessel_j.c",
"bessel_k.c",
+ "bessel_y.c",
);
# These are for plugin fn-R. They are so small it makes no sense to place
@@ -137,7 +139,7 @@ sub import_file {
s/\bISNAN\b/gnm_isnan/g;
s/\bR_(finite|FINITE)\b/gnm_finite/g;
- s/\b(sqrt|exp|log|pow|log1p|expm1|floor|ceil|sin|sinh|tan)(\s|$|\()/gnm_$1$2/g;
+ s/\b(sqrt|exp|log|pow|log1p|expm1|floor|ceil|sin|cos|sinh|tan)(\s|$|\()/gnm_$1$2/g;
s/\bfabs\b/gnm_abs/g;
s/\bgnm_floor\s*\(\s*([a-z]+)\s*\+\s*1e-7\s*\)/gnm_fake_floor($1)/;
@@ -152,7 +154,7 @@ sub import_file {
s~([-+]?\d+\.\d*)(\s*/\s*)([-+e0-9.]+)~GNM_const\($1\)$2$3~;
# These are made static.
-
s/^gnm_float\s+(pbeta_both|bd0|dpois_raw|chebyshev_eval|lgammacor|lbeta|pbeta_raw|dbinom_raw)/static
gnm_float $1/;
+
s/^gnm_float\s+(pbeta_both|bd0|dpois_raw|chebyshev_eval|lgammacor|lbeta|pbeta_raw|dbinom_raw|bessel_j_ex|bessel_y_ex)/static
gnm_float $1/;
s/^int\s+(chebyshev_init)/static int $1/;
# Fix a bug...
@@ -184,11 +186,24 @@ sub import_file {
s/\bgamma_cody\(empal\)/gnm_exp(lgamma1p(nu))/;
}
+ s/\bgamma_cody\b/gnm_gamma/g;
+
if (/^(static\s+)?(gnm_float|int)\s+(chebyshev_init)\s*\(/ .. /^\}/) {
next unless s/^(static\s+)?(gnm_float|int)\s+([a-zA-Z0-9_]+)\s*\(.*//;
$_ = "/* Definition of function $3 removed. */";
}
+ if ($filename =~ m|/bessel_j\.c$|) {
+ if (/^\s*const static gnm_float fact\[/) {
+ while (!/;$/) {
+ $_ .= <FIL>;
+ }
+ $_ = '/* removed array fact */';
+ } else {
+ s/\bfact\s*\[([^][]*)\]/fact($1)/g;
+ }
+ }
+
# Somewhat risky.
s/\%([-0-9.]*)([efgEFG])/\%$1\" GNM_FORMAT_$2 \"/g;
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
