#ifndef lint
static char *RCSid() { return RCSid("$Id: specfun.c,v 1.7.4.1 2000/06/22 12:57:39 broeker Exp $"); }
#endif
/* GNUPLOT - specfun.c */
/*[
* Copyright 1986 - 1993, 1998 Thomas Williams, Colin Kelley
*
* Permission to use, copy, and distribute this software and its
* documentation for any purpose with or without fee is hereby granted,
* provided that the above copyright notice appear in all copies and
* that both that copyright notice and this permission notice appear
* in supporting documentation.
*
* Permission to modify the software is granted, but not the right to
* distribute the complete modified source code. Modifications are to
* be distributed as patches to the released version. Permission to
* distribute binaries produced by compiling modified sources is granted,
* provided you
* 1. distribute the corresponding source modifications from the
* released version in the form of a patch file along with the binaries,
* 2. add special version identification to distinguish your version
* in addition to the base release version number,
* 3. provide your name and address as the primary contact for the
* support of your modified version, and
* 4. retain our contact information in regard to use of the base
* software.
* Permission to distribute the released version of the source code along
* with corresponding source modifications in the form of a patch file is
* granted with same provisions 2 through 4 for binary distributions.
*
* This software is provided "as is" without express or implied warranty
* to the extent permitted by applicable law.
]*/
/*
* AUTHORS
*
* Original Software:
* Jos van der Woude, jvdwoude@hut.nl
*
*/
#include "specfun.h"
#include "internal.h"
/* #include "setshow.h" */
#include "util.h"
#define ITMAX 200
#ifdef FLT_EPSILON
# define MACHEPS FLT_EPSILON /* 1.0E-08 */
#else
# define MACHEPS 1.0E-08
#endif
/* AS239 value, e^-88 = 2^-127 */
#define MINEXP -88.0
#ifdef FLT_MAX
# define OFLOW FLT_MAX /* 1.0E+37 */
#else
# define OFLOW 1.0E+37
#endif
/* AS239 value for igamma(a,x>=XBIG) = 1.0 */
#define XBIG 1.0E+08
/*
* Mathematical constants
*/
#define LNPI 1.14472988584940016
#define LNSQRT2PI 0.9189385332046727
#define PNT68 0.6796875
#define SQRT_TWO 1.41421356237309504880168872420969809 /* JG */
/* Prefer lgamma */
#ifndef GAMMA
# ifdef HAVE_LGAMMA
# define GAMMA(x) lgamma (x)
# elif defined(HAVE_GAMMA)
# define GAMMA(x) gamma (x)
# else
# undef GAMMA
# endif
#endif
#ifndef GAMMA
int signgam = 0;
#else
extern int signgam; /* this is not always declared in math.h */
#endif
/* Global variables, not visible outside this file */
static long Xm1 = 2147483563L;
static long Xm2 = 2147483399L;
static long Xa1 = 40014L;
static long Xa2 = 40692L;
/* Local function declarations, not visible outside this file */
static double confrac __PROTO((double a, double b, double x));
static double ibeta __PROTO((double a, double b, double x));
static double igamma __PROTO((double a, double x));
static double ranf __PROTO((double init));
static double inverse_normal_func __PROTO((double p));
static double inverse_error_func __PROTO((double p));
#ifndef GAMMA
/* Provide GAMMA function for those who do not already have one */
static double lngamma __PROTO((double z));
static double lgamneg __PROTO((double z));
static double lgampos __PROTO((double z));
/**
* from statlib, Thu Jan 23 15:02:27 EST 1992 ***
*
* This file contains two algorithms for the logarithm of the gamma function.
* Algorithm AS 245 is the faster (but longer) and gives an accuracy of about
* 10-12 significant decimal digits except for small regions around X = 1 and
* X = 2, where the function goes to zero.
* The second algorithm is not part of the AS algorithms. It is slower but
* gives 14 or more significant decimal digits accuracy, except around X = 1
* and X = 2. The Lanczos series from which this algorithm is derived is
* interesting in that it is a convergent series approximation for the gamma
* function, whereas the familiar series due to De Moivre (and usually wrongly
* called Stirling's approximation) is only an asymptotic approximation, as
* is the true and preferable approximation due to Stirling.
*
* Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos,
* C. 'A precision approximation of the gamma function', J. SIAM Numer.
* Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for
* small regions in the vicinity of 1 and 2.
*
* Programmer: Alan Miller CSIRO Division of Mathematics & Statistics
*
* Latest revision - 17 April 1988
*
* Additions: Translated from fortran to C, code added to handle values z < 0.
* The global variable signgam contains the sign of the gamma function.
*
* IMPORTANT: The signgam variable contains garbage until AFTER the call to
* lngamma().
*
* Permission granted to distribute freely for non-commercial purposes only
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
*/
/* Local data, not visible outside this file
static double a[] =
{
0.9999999999995183E+00,
0.6765203681218835E+03,
-.1259139216722289E+04,
0.7713234287757674E+03,
-.1766150291498386E+03,
0.1250734324009056E+02,
-.1385710331296526E+00,
0.9934937113930748E-05,
0.1659470187408462E-06,
}; */
/* from Ray Toy */
static double GPFAR a[] =
{
.99999999999980993227684700473478296744476168282198,
676.52036812188509856700919044401903816411251975244084,
-1259.13921672240287047156078755282840836424300664868028,
771.32342877765307884865282588943070775227268469602500,
-176.61502916214059906584551353999392943274507608117860,
12.50734327868690481445893685327104972970563021816420,
-.13857109526572011689554706984971501358032683492780,
.00000998436957801957085956266828104544089848531228,
.00000015056327351493115583383579667028994545044040,
};
static double
lgamneg(z)
double z;
{
double tmp;
/* Use reflection formula, then call lgampos() */
tmp = sin(z * M_PI);
if (fabs(tmp) < MACHEPS) {
tmp = 0.0;
} else if (tmp < 0.0) {
tmp = -tmp;
signgam = -1;
}
return LNPI - lgampos(1.0 - z) - log(tmp);
}
static double
lgampos(z)
double z;
{
double sum;
double tmp;
int i;
sum = a[0];
for (i = 1, tmp = z; i < 9; i++) {
sum += a[i] / tmp;
tmp++;
}
return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5);
}
static double
lngamma(z)
double z;
{
signgam = 1;
if (z = 0 ? 1 : 0;
x = igamma(0.5, (x) * (x));
if (x == -1.0) {
undefined = TRUE;
x = 0.0;
} else {
if (fsign == 0)
x = -x;
}
}
#endif
push(Gcomplex(&a, x, 0.0));
}
void
f_erfc()
{
struct value a;
double x;
x = real(pop(&a));
#ifdef HAVE_ERFC
x = erfc(x);
#else
{
int fsign;
fsign = x >= 0 ? 1 : 0;
x = igamma(0.5, (x) * (x));
if (x == 1.0) {
undefined = TRUE;
x = 0.0;
} else {
x = fsign > 0 ? 1.0 - x : 1.0 + x;
}
}
#endif
push(Gcomplex(&a, x, 0.0));
}
void
f_ibeta()
{
struct value a;
double x;
double arg1;
double arg2;
x = real(pop(&a));
arg2 = real(pop(&a));
arg1 = real(pop(&a));
x = ibeta(arg1, arg2, x);
if (x == -1.0) {
undefined = TRUE;
push(Ginteger(&a, 0));
} else
push(Gcomplex(&a, x, 0.0));
}
void
f_igamma()
{
struct value a;
double x;
double arg1;
x = real(pop(&a));
arg1 = real(pop(&a));
x = igamma(arg1, x);
if (x == -1.0) {
undefined = TRUE;
push(Ginteger(&a, 0));
} else
push(Gcomplex(&a, x, 0.0));
}
void
f_gamma()
{
register double y;
struct value a;
y = GAMMA(real(pop(&a)));
if (y > 88.0) {
undefined = TRUE;
push(Ginteger(&a, 0));
} else
push(Gcomplex(&a, signgam * gp_exp(y), 0.0));
}
void
f_lgamma()
{
struct value a;
push(Gcomplex(&a, GAMMA(real(pop(&a))), 0.0));
}
#ifndef BADRAND
void
f_rand()
{
struct value a;
push(Gcomplex(&a, ranf(real(pop(&a))), 0.0));
}
#else /* BADRAND */
/* Use only to observe the effect of a "bad" random number generator. */
void
f_rand()
{
struct value a;
static unsigned int y = 0;
unsigned int maxran = 1000;
(void) real(pop(&a));
y = (781 * y + 387) % maxran;
push(Gcomplex(&a, (double) y / maxran, 0.0));
}
#endif /* BADRAND */
/* ** ibeta.c
*
* DESCRIB Approximate the incomplete beta function Ix(a, b).
*
* _
* |(a + b) /x (a-1) (b-1)
* Ix(a, b) = -_-------_--- * | t * (1 - t) dt (a,b > 0)
* |(a) * |(b) /0
*
*
*
* CALL p = ibeta(a, b, x)
*
* double a > 0
* double b > 0
* double x [0, 1]
*
* WARNING none
*
* RETURN double p [0, 1]
* -1.0 on error condition
*
* XREF lngamma()
*
* BUGS none
*
* REFERENCE The continued fraction expansion as given by
* Abramowitz and Stegun (1964) is used.
*
* Permission granted to distribute freely for non-commercial purposes only
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
*/
static double
ibeta(a, b, x)
double a, b, x;
{
/* Test for admissibility of arguments */
if (a 0
* double x >= 0
*
* WARNING none
*
* RETURN double p [0, 1]
* -1.0 on error condition
*
* XREF lngamma()
*
* BUGS Values 0 = a + 2.0)) {
/* Use a continued fraction expansion */
double rn;
double rnold;
aa = 1.0 - a;
b = aa + x + 1.0;
pn1 = 1.0;
pn2 = x;
pn3 = x + 1.0;
pn4 = x * b;
rnold = pn3 / pn4;
for (i = 1; i