version 1.1.1.2, 2000/01/22 14:16:01 |
version 1.1.1.3, 2003/09/15 07:09:26 |
|
|
static char *RCSid = "$Id$"; |
static char *RCSid = "$Id$"; |
#endif |
#endif |
|
|
|
|
/* GNUPLOT - specfun.c */ |
/* GNUPLOT - specfun.c */ |
|
|
/*[ |
/*[ |
Line 43 static char *RCSid = "$Id$"; |
|
Line 42 static char *RCSid = "$Id$"; |
|
* |
* |
*/ |
*/ |
|
|
|
|
#include "plot.h" |
#include "plot.h" |
#include "fnproto.h" |
#include "fnproto.h" |
|
|
|
|
extern struct value stack[STACK_DEPTH]; |
extern struct value stack[STACK_DEPTH]; |
extern int s_p; |
extern int s_p; |
extern double zero; |
extern double zero; |
Line 107 static long Xa1 = 40014L; |
|
Line 106 static long Xa1 = 40014L; |
|
static long Xa2 = 40692L; |
static long Xa2 = 40692L; |
|
|
/* Local function declarations, not visible outside this file */ |
/* Local function declarations, not visible outside this file */ |
|
static int mtherr(char *name, int code); |
|
static double polevl(double x, double coef[], int N); |
|
static double p1evl(double x, double coef[], int N); |
static double confrac __PROTO((double a, double b, double x)); |
static double confrac __PROTO((double a, double b, double x)); |
static double ibeta __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 igamma __PROTO((double a, double x)); |
Line 114 static double ranf __PROTO((double init)); |
|
Line 116 static double ranf __PROTO((double init)); |
|
static double inverse_normal_func __PROTO((double p)); |
static double inverse_normal_func __PROTO((double p)); |
static double inverse_error_func __PROTO((double p)); |
static double inverse_error_func __PROTO((double p)); |
|
|
#ifndef GAMMA |
/* UNKnown arithmetic, invokes coefficients given in |
/* Provide GAMMA function for those who do not already have one */ |
* normal decimal format. Beware of range boundary |
static double lngamma __PROTO((double z)); |
* problems (MACHEP, MAXLOG, etc. in const.c) and |
static double lgamneg __PROTO((double z)); |
* roundoff problems in pow.c: |
static double lgampos __PROTO((double z)); |
* (Sun SPARCstation) |
|
*/ |
|
#define UNK 1 |
|
|
/** |
/* If you define UNK, then be sure to set BIGENDIAN properly. */ |
* from statlib, Thu Jan 23 15:02:27 EST 1992 *** |
#ifdef FLOAT_WORDS_BIGENDIAN |
* |
#define BIGENDIAN 1 |
* This file contains two algorithms for the logarithm of the gamma function. |
#else |
* Algorithm AS 245 is the faster (but longer) and gives an accuracy of about |
#define BIGENDIAN 0 |
* 10-12 significant decimal digits except for small regions around X = 1 and |
#endif |
* X = 2, where the function goes to zero. |
|
* The second algorithm is not part of the AS algorithms. It is slower but |
/* Define to support tiny denormal numbers, else undefine. */ |
* gives 14 or more significant decimal digits accuracy, except around X = 1 |
#define DENORMAL 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 |
/* Define to ask for infinity support, else undefine. */ |
* function, whereas the familiar series due to De Moivre (and usually wrongly |
#define INFINITIES 1 |
* called Stirling's approximation) is only an asymptotic approximation, as |
|
* is the true and preferable approximation due to Stirling. |
/* Define to ask for support of numbers that are Not-a-Number, |
* |
else undefine. This may automatically define INFINITIES in some files. */ |
* Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos, |
#define NANS 1 |
* C. 'A precision approximation of the gamma function', J. SIAM Numer. |
|
* Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for |
/* Define to distinguish between -0.0 and +0.0. */ |
* small regions in the vicinity of 1 and 2. |
#define MINUSZERO 1 |
* |
|
* Programmer: Alan Miller CSIRO Division of Mathematics & Statistics |
/* Variable for error reporting. See mtherr.c. */ |
* |
extern int merror; |
* Latest revision - 17 April 1988 |
|
* |
/* |
* Additions: Translated from fortran to C, code added to handle values z < 0. |
Cephes Math Library Release 2.0: April, 1987 |
* The global variable signgam contains the sign of the gamma function. |
Copyright 1984, 1987 by Stephen L. Moshier |
* |
Direct inquiries to 30 Frost Street, Cambridge, MA 02140 |
* IMPORTANT: The signgam variable contains garbage until AFTER the call to |
*/ |
* lngamma(). |
|
* |
#include <stdio.h> |
* Permission granted to distribute freely for non-commercial purposes only |
|
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl |
int merror = 0; |
|
|
|
/* Notice: the order of appearance of the following |
|
* messages is bound to the error codes defined |
|
* in mconf.h. |
*/ |
*/ |
|
static char *ermsg[7] = { |
|
"unknown", /* error code 0 */ |
|
"domain", /* error code 1 */ |
|
"singularity", /* et seq. */ |
|
"overflow", |
|
"underflow", |
|
"total loss of precision", |
|
"partial loss of precision" |
|
}; |
|
|
/* Local data, not visible outside this file |
#ifndef DOMAIN |
static double a[] = |
/* HBB 20021103: copied this from Cephes mconf.h */ |
{ |
# define DOMAIN 1 /* argument domain error */ |
0.9999999999995183E+00, |
# define SING 2 /* argument singularity */ |
0.6765203681218835E+03, |
# define OVERFLOW 3 /* overflow range error */ |
-.1259139216722289E+04, |
# define UNDERFLOW 4 /* underflow range error */ |
0.7713234287757674E+03, |
# define TLOSS 5 /* total loss of precision */ |
-.1766150291498386E+03, |
# define PLOSS 6 /* partial loss of precision */ |
0.1250734324009056E+02, |
#endif |
-.1385710331296526E+00, |
|
0.9934937113930748E-05, |
|
0.1659470187408462E-06, |
|
}; */ |
|
|
|
/* from Ray Toy */ |
|
static double GPFAR a[] = |
|
|
static int mtherr(char *name, int code) |
{ |
{ |
.99999999999980993227684700473478296744476168282198, |
|
676.52036812188509856700919044401903816411251975244084, |
|
-1259.13921672240287047156078755282840836424300664868028, |
|
771.32342877765307884865282588943070775227268469602500, |
|
-176.61502916214059906584551353999392943274507608117860, |
|
12.50734327868690481445893685327104972970563021816420, |
|
-.13857109526572011689554706984971501358032683492780, |
|
.00000998436957801957085956266828104544089848531228, |
|
.00000015056327351493115583383579667028994545044040, |
|
}; |
|
|
|
static double lgamneg(z) |
/* Display string passed by calling program, |
double z; |
* which is supposed to be the name of the |
|
* function in which the error occurred: |
|
*/ |
|
printf("\n%s ", name); |
|
|
|
/* Set global error message word */ |
|
merror = code; |
|
|
|
/* Display error message defined |
|
* by the code argument. |
|
*/ |
|
if ((code <= 0) || (code >= 7)) |
|
code = 0; |
|
printf("%s error\n", ermsg[code]); |
|
|
|
/* Return to calling |
|
* program |
|
*/ |
|
return (0); |
|
} |
|
|
|
/* polevl.c |
|
* p1evl.c |
|
* |
|
* Evaluate polynomial |
|
* |
|
* |
|
* |
|
* SYNOPSIS: |
|
* |
|
* int N; |
|
* double x, y, coef[N+1], polevl[]; |
|
* |
|
* y = polevl( x, coef, N ); |
|
* |
|
* |
|
* |
|
* DESCRIPTION: |
|
* |
|
* Evaluates polynomial of degree N: |
|
* |
|
* 2 N |
|
* y = C + C x + C x +...+ C x |
|
* 0 1 2 N |
|
* |
|
* Coefficients are stored in reverse order: |
|
* |
|
* coef[0] = C , ..., coef[N] = C . |
|
* N 0 |
|
* |
|
* The function p1evl() assumes that coef[N] = 1.0 and is |
|
* omitted from the array. Its calling arguments are |
|
* otherwise the same as polevl(). |
|
* |
|
* |
|
* SPEED: |
|
* |
|
* In the interest of speed, there are no checks for out |
|
* of bounds arithmetic. This routine is used by most of |
|
* the functions in the library. Depending on available |
|
* equipment features, the user may wish to rewrite the |
|
* program in microcode or assembly language. |
|
* |
|
*/ |
|
|
|
/* |
|
Cephes Math Library Release 2.1: December, 1988 |
|
Copyright 1984, 1987, 1988 by Stephen L. Moshier |
|
Direct inquiries to 30 Frost Street, Cambridge, MA 02140 |
|
*/ |
|
|
|
static double polevl(double x, double coef[], int N) |
{ |
{ |
double tmp; |
double ans; |
|
int i; |
|
double *p; |
|
|
/* Use reflection formula, then call lgampos() */ |
p = coef; |
tmp = sin(z * PI); |
ans = *p++; |
|
i = N; |
|
|
if (fabs(tmp) < MACHEPS) { |
do |
tmp = 0.0; |
ans = ans * x + *p++; |
} else if (tmp < 0.0) { |
while (--i); |
tmp = -tmp; |
|
signgam = -1; |
|
} |
|
return LNPI - lgampos(1.0 - z) - log(tmp); |
|
|
|
|
return (ans); |
} |
} |
|
|
static double lgampos(z) |
/* p1evl() */ |
double z; |
/* N |
|
* Evaluate polynomial when coefficient of x is 1.0. |
|
* Otherwise same as polevl. |
|
*/ |
|
|
|
static double p1evl(double x, double coef[], int N) |
{ |
{ |
double sum; |
double ans; |
double tmp; |
double *p; |
int i; |
int i; |
|
|
sum = a[0]; |
p = coef; |
for (i = 1, tmp = z; i < 9; i++) { |
ans = x + *p++; |
sum += a[i] / tmp; |
i = N - 1; |
tmp++; |
|
} |
|
|
|
return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5); |
do |
|
ans = ans * x + *p++; |
|
while (--i); |
|
|
|
return (ans); |
} |
} |
|
|
static double lngamma(z) |
#ifndef GAMMA |
double z; |
/* Provide GAMMA function for those who do not already have one */ |
|
static double lngamma __PROTO((double z)); |
|
|
|
int sgngam; |
|
/* A[]: Stirling's formula expansion of log gamma |
|
* B[], C[]: log gamma function between 2 and 3 |
|
*/ |
|
#ifdef UNK |
|
static double A[] = { |
|
8.11614167470508450300E-4, |
|
-5.95061904284301438324E-4, |
|
7.93650340457716943945E-4, |
|
-2.77777777730099687205E-3, |
|
8.33333333333331927722E-2 |
|
}; |
|
static double B[] = { |
|
-1.37825152569120859100E3, |
|
-3.88016315134637840924E4, |
|
-3.31612992738871184744E5, |
|
-1.16237097492762307383E6, |
|
-1.72173700820839662146E6, |
|
-8.53555664245765465627E5 |
|
}; |
|
static double C[] = { |
|
/* 1.00000000000000000000E0, */ |
|
-3.51815701436523470549E2, |
|
-1.70642106651881159223E4, |
|
-2.20528590553854454839E5, |
|
-1.13933444367982507207E6, |
|
-2.53252307177582951285E6, |
|
-2.01889141433532773231E6 |
|
}; |
|
/* log( sqrt( 2*pi ) ) */ |
|
static double LS2PI = 0.91893853320467274178; |
|
#define MAXLGM 2.556348e305 |
|
#endif |
|
#ifdef DEC |
|
static unsigned short A[] = { |
|
0035524, 0141201, 0034633, 0031405, |
|
0135433, 0176755, 0126007, 0045030, |
|
0035520, 0006371, 0003342, 0172730, |
|
0136066, 0005540, 0132605, 0026407, |
|
0037252, 0125252, 0125252, 0125132 |
|
}; |
|
static unsigned short B[] = { |
|
0142654, 0044014, 0077633, 0035410, |
|
0144027, 0110641, 0125335, 0144760, |
|
0144641, 0165637, 0142204, 0047447, |
|
0145215, 0162027, 0146246, 0155211, |
|
0145322, 0026110, 0010317, 0110130, |
|
0145120, 0061472, 0120300, 0025363 |
|
}; |
|
static unsigned short C[] = { |
|
/*0040200,0000000,0000000,0000000*/ |
|
0142257, 0164150, 0163630, 0112622, |
|
0143605, 0050153, 0156116, 0135272, |
|
0144527, 0056045, 0145642, 0062332, |
|
0145213, 0012063, 0106250, 0001025, |
|
0145432, 0111254, 0044577, 0115142, |
|
0145366, 0071133, 0050217, 0005122 |
|
}; |
|
/* log( sqrt( 2*pi ) ) */ |
|
static unsigned short LS2P[] = {040153, 037616, 041445, 0172645,}; |
|
#define LS2PI *(double *)LS2P |
|
#define MAXLGM 2.035093e36 |
|
#endif |
|
|
|
#ifdef IBMPC |
|
static unsigned short A[] = { |
|
0x6661, 0x2733, 0x9850, 0x3f4a, |
|
0xe943, 0xb580, 0x7fbd, 0xbf43, |
|
0x5ebb, 0x20dc, 0x019f, 0x3f4a, |
|
0xa5a1, 0x16b0, 0xc16c, 0xbf66, |
|
0x554b, 0x5555, 0x5555, 0x3fb5 |
|
}; |
|
static unsigned short B[] = { |
|
0x6761, 0x8ff3, 0x8901, 0xc095, |
|
0xb93e, 0x355b, 0xf234, 0xc0e2, |
|
0x89e5, 0xf890, 0x3d73, 0xc114, |
|
0xdb51, 0xf994, 0xbc82, 0xc131, |
|
0xf20b, 0x0219, 0x4589, 0xc13a, |
|
0x055e, 0x5418, 0x0c67, 0xc12a |
|
}; |
|
static unsigned short C[] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0x12b2, 0x1cf3, 0xfd0d, 0xc075, |
|
0xd757, 0x7b89, 0xaa0d, 0xc0d0, |
|
0x4c9b, 0xb974, 0xeb84, 0xc10a, |
|
0x0043, 0x7195, 0x6286, 0xc131, |
|
0xf34c, 0x892f, 0x5255, 0xc143, |
|
0xe14a, 0x6a11, 0xce4b, 0xc13e |
|
}; |
|
/* log( sqrt( 2*pi ) ) */ |
|
static unsigned short LS2P[] = { |
|
0xbeb5, 0xc864, 0x67f1, 0x3fed |
|
}; |
|
#define LS2PI *(double *)LS2P |
|
#define MAXLGM 2.556348e305 |
|
#endif |
|
|
|
#ifdef MIEEE |
|
static unsigned short A[] = { |
|
0x3f4a, 0x9850, 0x2733, 0x6661, |
|
0xbf43, 0x7fbd, 0xb580, 0xe943, |
|
0x3f4a, 0x019f, 0x20dc, 0x5ebb, |
|
0xbf66, 0xc16c, 0x16b0, 0xa5a1, |
|
0x3fb5, 0x5555, 0x5555, 0x554b |
|
}; |
|
static unsigned short B[] = { |
|
0xc095, 0x8901, 0x8ff3, 0x6761, |
|
0xc0e2, 0xf234, 0x355b, 0xb93e, |
|
0xc114, 0x3d73, 0xf890, 0x89e5, |
|
0xc131, 0xbc82, 0xf994, 0xdb51, |
|
0xc13a, 0x4589, 0x0219, 0xf20b, |
|
0xc12a, 0x0c67, 0x5418, 0x055e |
|
}; |
|
static unsigned short C[] = { |
|
0xc075, 0xfd0d, 0x1cf3, 0x12b2, |
|
0xc0d0, 0xaa0d, 0x7b89, 0xd757, |
|
0xc10a, 0xeb84, 0xb974, 0x4c9b, |
|
0xc131, 0x6286, 0x7195, 0x0043, |
|
0xc143, 0x5255, 0x892f, 0xf34c, |
|
0xc13e, 0xce4b, 0x6a11, 0xe14a |
|
}; |
|
/* log( sqrt( 2*pi ) ) */ |
|
static unsigned short LS2P[] = { |
|
0x3fed, 0x67f1, 0xc864, 0xbeb5 |
|
}; |
|
#define LS2PI *(double *)LS2P |
|
#define MAXLGM 2.556348e305 |
|
#endif |
|
|
|
/* static const double PI = 3.1415926535897932384626433832795028841971693993751;*/ |
|
static const double LOGPI = 1.1447298858494001741434273513530587116472948129153; |
|
|
|
int ISNAN(double x) |
{ |
{ |
signgam = 1; |
volatile double a = x; |
|
if (a != a) |
|
return 1; |
|
return 0; |
|
} |
|
int ISFINITE(double x) |
|
{ |
|
volatile double a = x; |
|
if (a < DBL_MAX) |
|
return 1; |
|
return 0; |
|
} |
|
|
if (z <= 0.0) |
double lngamma(double x) |
return lgamneg(z); |
{ |
|
double p, |
|
q, |
|
u, |
|
w, |
|
z; |
|
int i; |
|
|
|
sgngam = 1; |
|
#ifdef NANS |
|
if (ISNAN(x)) |
|
return (x); |
|
#endif |
|
|
|
#ifdef INFINITIES |
|
if (!ISFINITE((x))) |
|
return (DBL_MAX * DBL_MAX); |
|
#endif |
|
|
|
if (x < -34.0) { |
|
q = -x; |
|
w = lngamma(q); /* note this modifies sgngam! */ |
|
p = floor(q); |
|
if (p == q) { |
|
lgsing: |
|
#ifdef INFINITIES |
|
mtherr("lngamma", SING); |
|
return (DBL_MAX * DBL_MAX); |
|
#else |
|
goto loverf; |
|
#endif |
|
} |
|
i = p; |
|
if ((i & 1) == 0) |
|
sgngam = -1; |
|
else |
|
sgngam = 1; |
|
z = q - p; |
|
if (z > 0.5) { |
|
p += 1.0; |
|
z = p - q; |
|
} |
|
z = q * sin(PI * z); |
|
if (z == 0.0) |
|
goto lgsing; |
|
/* z = log(PI) - log( z ) - w;*/ |
|
z = LOGPI - log(z) - w; |
|
return (z); |
|
} |
|
if (x < 13.0) { |
|
z = 1.0; |
|
p = 0.0; |
|
u = x; |
|
while (u >= 3.0) { |
|
p -= 1.0; |
|
u = x + p; |
|
z *= u; |
|
} |
|
while (u < 2.0) { |
|
if (u == 0.0) |
|
goto lgsing; |
|
z /= u; |
|
p += 1.0; |
|
u = x + p; |
|
} |
|
if (z < 0.0) { |
|
sgngam = -1; |
|
z = -z; |
|
} else |
|
sgngam = 1; |
|
if (u == 2.0) |
|
return (log(z)); |
|
p -= 2.0; |
|
x = x + p; |
|
p = x * polevl(x, B, 5) / p1evl(x, C, 6); |
|
return (log(z) + p); |
|
} |
|
if (x > MAXLGM) { |
|
#ifdef INFINITIES |
|
return (sgngam * (DBL_MAX * DBL_MAX)); |
|
#else |
|
loverf: |
|
mtherr("lngamma", OVERFLOW); |
|
return (sgngam * MAXNUM); |
|
#endif |
|
} |
|
q = (x - 0.5) * log(x) - x + LS2PI; |
|
if (x > 1.0e8) |
|
return (q); |
|
|
|
p = 1.0 / (x * x); |
|
if (x >= 1000.0) |
|
q += ((7.9365079365079365079365e-4 * p |
|
- 2.7777777777777777777778e-3) * p |
|
+ 0.0833333333333333333333) / x; |
else |
else |
return lgampos(z); |
q += polevl(p, A, 4) / x; |
|
return (q); |
} |
} |
|
|
# define GAMMA(x) lngamma ((x)) |
#define GAMMA(x) lngamma ((x)) |
#endif /* !GAMMA */ |
#endif /* !GAMMA */ |
|
|
void f_erf() |
void f_erf() |
|
|
|
|
#endif /* BADRAND */ |
#endif /* BADRAND */ |
|
|
/** ibeta.c |
/* ** ibeta.c |
* |
* |
* DESCRIB Approximate the incomplete beta function Ix(a, b). |
* DESCRIB Approximate the incomplete beta function Ix(a, b). |
* |
* |
|
|
* REFERENCE The continued fraction expansion as given by |
* REFERENCE The continued fraction expansion as given by |
* Abramowitz and Stegun (1964) is used. |
* 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 |
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl |
|
* |
|
* Note: this function was translated from the Public Domain Fortran |
|
* version available from http://lib.stat.cmu.edu/apstat/xxx |
|
* |
*/ |
*/ |
|
|
static double ibeta(a, b, x) |
static double ibeta(a, b, x) |
|
|
return -1.0; |
return -1.0; |
} |
} |
|
|
/** igamma.c |
/* ** igamma.c |
* |
* |
* DESCRIB Approximate the incomplete gamma function P(a, x). |
* DESCRIB Approximate the incomplete gamma function P(a, x). |
* |
* |
|
|
* |
* |
* REFERENCE ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3 |
* REFERENCE ALGORITHM AS239 APPL. STATIST. (1988) VOL. 37, NO. 3 |
* |
* |
* Permission granted to distribute freely for non-commercial purposes only |
|
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl |
* Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl |
|
* |
|
* Note: this function was translated from the Public Domain Fortran |
|
* version available from http://lib.stat.cmu.edu/apstat/239 |
|
* |
*/ |
*/ |
|
|
/* Global variables, not visible outside this file */ |
/* Global variables, not visible outside this file */ |
Line 704 void f_inverse_erf() |
|
Line 1033 void f_inverse_erf() |
|
} |
} |
} |
} |
|
|
static double inverse_normal_func(p) |
/* ndtri.c |
double p; |
* |
{ |
* Inverse of Normal distribution function |
/* |
* |
Source: This routine was derived (using f2c) from the |
* |
FORTRAN subroutine MDNRIS found in |
* |
ACM Algorithm 602 obtained from netlib. |
* SYNOPSIS: |
|
* |
|
* double x, y, ndtri(); |
|
* |
|
* x = ndtri( y ); |
|
* |
|
* |
|
* |
|
* DESCRIPTION: |
|
* |
|
* Returns the argument, x, for which the area under the |
|
* Gaussian probability density function (integrated from |
|
* minus infinity to x) is equal to y. |
|
* |
|
* |
|
* For small arguments 0 < y < exp(-2), the program computes |
|
* z = sqrt( -2.0 * log(y) ); then the approximation is |
|
* x = z - log(z)/z - (1/z) P(1/z) / Q(1/z). |
|
* There are two rational functions P/Q, one for 0 < y < exp(-32) |
|
* and the other for y up to exp(-2). For larger arguments, |
|
* w = y - 0.5, and x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2)). |
|
* |
|
* |
|
* ACCURACY: |
|
* |
|
* Relative error: |
|
* arithmetic domain # trials peak rms |
|
* DEC 0.125, 1 5500 9.5e-17 2.1e-17 |
|
* DEC 6e-39, 0.135 3500 5.7e-17 1.3e-17 |
|
* IEEE 0.125, 1 20000 7.2e-16 1.3e-16 |
|
* IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17 |
|
* |
|
* |
|
* ERROR MESSAGES: |
|
* |
|
* message condition value returned |
|
* ndtri domain x <= 0 -DBL_MAX |
|
* ndtri domain x >= 1 DBL_MAX |
|
* |
|
*/ |
|
|
MDNRIS code contains the 1978 Copyright |
/* |
by IMSL, INC. . Since MDNRIS has been |
Cephes Math Library Release 2.8: June, 2000 |
submitted to netlib it may be used with |
Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier |
the restriction that it may only be |
*/ |
used for noncommercial purposes and that |
|
IMSL be acknowledged as the copyright-holder |
|
of the code. |
|
*/ |
|
|
|
/* Initialized data */ |
#ifdef UNK |
static double eps = 1e-10; |
/* sqrt(2pi) */ |
static double g0 = 1.851159e-4; |
static double s2pi = 2.50662827463100050242E0; |
static double g1 = -.002028152; |
#endif |
static double g2 = -.1498384; |
|
static double g3 = .01078639; |
|
static double h0 = .09952975; |
|
static double h1 = .5211733; |
|
static double h2 = -.06888301; |
|
static double sqrt2 = 1.414213562373095; |
|
|
|
/* Local variables */ |
#ifdef DEC |
static double a, w, x; |
static unsigned short s2p[] = {0040440, 0066230, 0177661, 0034055}; |
static double sd, wi, sn, y; |
#define s2pi *(double *)s2p |
|
#endif |
|
|
/* Note: 0.0 < p < 1.0 */ |
#ifdef IBMPC |
|
static unsigned short s2p[] = {0x2706, 0x1ff6, 0x0d93, 0x4004}; |
|
#define s2pi *(double *)s2p |
|
#endif |
|
|
/* p too small, compute y directly */ |
#ifdef MIEEE |
if (p <= eps) { |
static unsigned short s2p[] = { |
a = p + p; |
0x4004, 0x0d93, 0x1ff6, 0x2706 |
w = sqrt(-(double) log(a + (a - a * a))); |
}; |
|
#define s2pi *(double *)s2p |
|
#endif |
|
|
/* use a rational function in 1.0 / w */ |
/* approximation for 0 <= |y - 0.5| <= 3/8 */ |
wi = 1.0 / w; |
#ifdef UNK |
sn = ((g3 * wi + g2) * wi + g1) * wi; |
static double P0[5] = { |
sd = ((wi + h2) * wi + h1) * wi + h0; |
-5.99633501014107895267E1, |
y = w + w * (g0 + sn / sd); |
9.80010754185999661536E1, |
y = -y * sqrt2; |
-5.66762857469070293439E1, |
} else { |
1.39312609387279679503E1, |
x = 1.0 - (p + p); |
-1.23916583867381258016E0, |
y = inverse_error_func(x); |
}; |
y = -sqrt2 * y; |
static double Q0[8] = { |
|
/* 1.00000000000000000000E0,*/ |
|
1.95448858338141759834E0, |
|
4.67627912898881538453E0, |
|
8.63602421390890590575E1, |
|
-2.25462687854119370527E2, |
|
2.00260212380060660359E2, |
|
-8.20372256168333339912E1, |
|
1.59056225126211695515E1, |
|
-1.18331621121330003142E0, |
|
}; |
|
#endif |
|
#ifdef DEC |
|
static unsigned short P0[20] = { |
|
0141557, 0155170, 0071360, 0120550, |
|
0041704, 0000214, 0172417, 0067307, |
|
0141542, 0132204, 0040066, 0156723, |
|
0041136, 0163161, 0157276, 0007747, |
|
0140236, 0116374, 0073666, 0051764, |
|
}; |
|
static unsigned short Q0[32] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0040372, 0026256, 0110403, 0123707, |
|
0040625, 0122024, 0020277, 0026661, |
|
0041654, 0134161, 0124134, 0007244, |
|
0142141, 0073162, 0133021, 0131371, |
|
0042110, 0041235, 0043516, 0057767, |
|
0141644, 0011417, 0036155, 0137305, |
|
0041176, 0076556, 0004043, 0125430, |
|
0140227, 0073347, 0152776, 0067251, |
|
}; |
|
#endif |
|
#ifdef IBMPC |
|
static unsigned short P0[20] = { |
|
0x142d, 0x0e5e, 0xfb4f, 0xc04d, |
|
0xedd9, 0x9ea1, 0x8011, 0x4058, |
|
0xdbba, 0x8806, 0x5690, 0xc04c, |
|
0xc1fd, 0x3bd7, 0xdcce, 0x402b, |
|
0xca7e, 0x8ef6, 0xd39f, 0xbff3, |
|
}; |
|
static unsigned short Q0[36] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0x74f9, 0xd220, 0x4595, 0x3fff, |
|
0xe5b6, 0x8417, 0xb482, 0x4012, |
|
0x81d4, 0x350b, 0x970e, 0x4055, |
|
0x365f, 0x56c2, 0x2ece, 0xc06c, |
|
0xcbff, 0xa8e9, 0x0853, 0x4069, |
|
0xb7d9, 0xe78d, 0x8261, 0xc054, |
|
0x7563, 0xc104, 0xcfad, 0x402f, |
|
0xcdd5, 0xfabf, 0xeedc, 0xbff2, |
|
}; |
|
#endif |
|
#ifdef MIEEE |
|
static unsigned short P0[20] = { |
|
0xc04d, 0xfb4f, 0x0e5e, 0x142d, |
|
0x4058, 0x8011, 0x9ea1, 0xedd9, |
|
0xc04c, 0x5690, 0x8806, 0xdbba, |
|
0x402b, 0xdcce, 0x3bd7, 0xc1fd, |
|
0xbff3, 0xd39f, 0x8ef6, 0xca7e, |
|
}; |
|
static unsigned short Q0[32] = { |
|
/*0x3ff0,0x0000,0x0000,0x0000,*/ |
|
0x3fff, 0x4595, 0xd220, 0x74f9, |
|
0x4012, 0xb482, 0x8417, 0xe5b6, |
|
0x4055, 0x970e, 0x350b, 0x81d4, |
|
0xc06c, 0x2ece, 0x56c2, 0x365f, |
|
0x4069, 0x0853, 0xa8e9, 0xcbff, |
|
0xc054, 0x8261, 0xe78d, 0xb7d9, |
|
0x402f, 0xcfad, 0xc104, 0x7563, |
|
0xbff2, 0xeedc, 0xfabf, 0xcdd5, |
|
}; |
|
#endif |
|
|
|
/* Approximation for interval z = sqrt(-2 log y ) between 2 and 8 |
|
* i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14. |
|
*/ |
|
#ifdef UNK |
|
static double P1[9] = { |
|
4.05544892305962419923E0, |
|
3.15251094599893866154E1, |
|
5.71628192246421288162E1, |
|
4.40805073893200834700E1, |
|
1.46849561928858024014E1, |
|
2.18663306850790267539E0, |
|
-1.40256079171354495875E-1, |
|
-3.50424626827848203418E-2, |
|
-8.57456785154685413611E-4, |
|
}; |
|
static double Q1[8] = { |
|
/* 1.00000000000000000000E0,*/ |
|
1.57799883256466749731E1, |
|
4.53907635128879210584E1, |
|
4.13172038254672030440E1, |
|
1.50425385692907503408E1, |
|
2.50464946208309415979E0, |
|
-1.42182922854787788574E-1, |
|
-3.80806407691578277194E-2, |
|
-9.33259480895457427372E-4, |
|
}; |
|
#endif |
|
#ifdef DEC |
|
static unsigned short P1[36] = { |
|
0040601, 0143074, 0150744, 0073326, |
|
0041374, 0031554, 0113253, 0146016, |
|
0041544, 0123272, 0012463, 0176771, |
|
0041460, 0051160, 0103560, 0156511, |
|
0041152, 0172624, 0117772, 0030755, |
|
0040413, 0170713, 0151545, 0176413, |
|
0137417, 0117512, 0022154, 0131671, |
|
0137017, 0104257, 0071432, 0007072, |
|
0135540, 0143363, 0063137, 0036166, |
|
}; |
|
static unsigned short Q1[32] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0041174, 0075325, 0004736, 0120326, |
|
0041465, 0110044, 0047561, 0045567, |
|
0041445, 0042321, 0012142, 0030340, |
|
0041160, 0127074, 0166076, 0141051, |
|
0040440, 0046055, 0040745, 0150400, |
|
0137421, 0114146, 0067330, 0010621, |
|
0137033, 0175162, 0025555, 0114351, |
|
0135564, 0122773, 0145750, 0030357, |
|
}; |
|
#endif |
|
#ifdef IBMPC |
|
static unsigned short P1[36] = { |
|
0x8edb, 0x9a3c, 0x38c7, 0x4010, |
|
0x7982, 0x92d5, 0x866d, 0x403f, |
|
0x7fbf, 0x42a6, 0x94d7, 0x404c, |
|
0x1ba9, 0x10ee, 0x0a4e, 0x4046, |
|
0x463e, 0x93ff, 0x5eb2, 0x402d, |
|
0xbfa1, 0x7a6c, 0x7e39, 0x4001, |
|
0x9677, 0x448d, 0xf3e9, 0xbfc1, |
|
0x41c7, 0xee63, 0xf115, 0xbfa1, |
|
0xe78f, 0x6ccb, 0x18de, 0xbf4c, |
|
}; |
|
static unsigned short Q1[32] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0xd41b, 0xa13b, 0x8f5a, 0x402f, |
|
0x296f, 0x89ee, 0xb204, 0x4046, |
|
0x461c, 0x228c, 0xa89a, 0x4044, |
|
0xd845, 0x9d87, 0x15c7, 0x402e, |
|
0xba20, 0xa83c, 0x0985, 0x4004, |
|
0x0232, 0xcddb, 0x330c, 0xbfc2, |
|
0xb31d, 0x456d, 0x7f4e, 0xbfa3, |
|
0x061e, 0x797d, 0x94bf, 0xbf4e, |
|
}; |
|
#endif |
|
#ifdef MIEEE |
|
static unsigned short P1[36] = { |
|
0x4010, 0x38c7, 0x9a3c, 0x8edb, |
|
0x403f, 0x866d, 0x92d5, 0x7982, |
|
0x404c, 0x94d7, 0x42a6, 0x7fbf, |
|
0x4046, 0x0a4e, 0x10ee, 0x1ba9, |
|
0x402d, 0x5eb2, 0x93ff, 0x463e, |
|
0x4001, 0x7e39, 0x7a6c, 0xbfa1, |
|
0xbfc1, 0xf3e9, 0x448d, 0x9677, |
|
0xbfa1, 0xf115, 0xee63, 0x41c7, |
|
0xbf4c, 0x18de, 0x6ccb, 0xe78f, |
|
}; |
|
static unsigned short Q1[32] = { |
|
/*0x3ff0,0x0000,0x0000,0x0000,*/ |
|
0x402f, 0x8f5a, 0xa13b, 0xd41b, |
|
0x4046, 0xb204, 0x89ee, 0x296f, |
|
0x4044, 0xa89a, 0x228c, 0x461c, |
|
0x402e, 0x15c7, 0x9d87, 0xd845, |
|
0x4004, 0x0985, 0xa83c, 0xba20, |
|
0xbfc2, 0x330c, 0xcddb, 0x0232, |
|
0xbfa3, 0x7f4e, 0x456d, 0xb31d, |
|
0xbf4e, 0x94bf, 0x797d, 0x061e, |
|
}; |
|
#endif |
|
|
|
/* Approximation for interval z = sqrt(-2 log y ) between 8 and 64 |
|
* i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890. |
|
*/ |
|
|
|
#ifdef UNK |
|
static double P2[9] = { |
|
3.23774891776946035970E0, |
|
6.91522889068984211695E0, |
|
3.93881025292474443415E0, |
|
1.33303460815807542389E0, |
|
2.01485389549179081538E-1, |
|
1.23716634817820021358E-2, |
|
3.01581553508235416007E-4, |
|
2.65806974686737550832E-6, |
|
6.23974539184983293730E-9, |
|
}; |
|
static double Q2[8] = { |
|
/* 1.00000000000000000000E0,*/ |
|
6.02427039364742014255E0, |
|
3.67983563856160859403E0, |
|
1.37702099489081330271E0, |
|
2.16236993594496635890E-1, |
|
1.34204006088543189037E-2, |
|
3.28014464682127739104E-4, |
|
2.89247864745380683936E-6, |
|
6.79019408009981274425E-9, |
|
}; |
|
#endif |
|
#ifdef DEC |
|
static unsigned short P2[36] = { |
|
0040517, 0033507, 0036236, 0125641, |
|
0040735, 0044616, 0014473, 0140133, |
|
0040574, 0012567, 0114535, 0102541, |
|
0040252, 0120340, 0143474, 0150135, |
|
0037516, 0051057, 0115361, 0031211, |
|
0036512, 0131204, 0101511, 0125144, |
|
0035236, 0016627, 0043160, 0140216, |
|
0033462, 0060512, 0060141, 0010641, |
|
0031326, 0062541, 0101304, 0077706, |
|
}; |
|
static unsigned short Q2[32] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0040700, 0143322, 0132137, 0040501, |
|
0040553, 0101155, 0053221, 0140257, |
|
0040260, 0041071, 0052573, 0010004, |
|
0037535, 0066472, 0177261, 0162330, |
|
0036533, 0160475, 0066666, 0036132, |
|
0035253, 0174533, 0027771, 0044027, |
|
0033502, 0016147, 0117666, 0063671, |
|
0031351, 0047455, 0141663, 0054751, |
|
}; |
|
#endif |
|
#ifdef IBMPC |
|
static unsigned short P2[36] = { |
|
0xd574, 0xe793, 0xe6e8, 0x4009, |
|
0x780b, 0xc327, 0xa931, 0x401b, |
|
0xb0ac, 0xf32b, 0x82ae, 0x400f, |
|
0x9a0c, 0x18e7, 0x541c, 0x3ff5, |
|
0x2651, 0xf35e, 0xca45, 0x3fc9, |
|
0x354d, 0x9069, 0x5650, 0x3f89, |
|
0x1812, 0xe8ce, 0xc3b2, 0x3f33, |
|
0x2234, 0x4c0c, 0x4c29, 0x3ec6, |
|
0x8ff9, 0x3058, 0xccac, 0x3e3a, |
|
}; |
|
static unsigned short Q2[32] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0xe828, 0x568b, 0x18da, 0x4018, |
|
0x3816, 0xaad2, 0x704d, 0x400d, |
|
0x6200, 0x2aaf, 0x0847, 0x3ff6, |
|
0x3c9b, 0x5fd6, 0xada7, 0x3fcb, |
|
0xc78b, 0xadb6, 0x7c27, 0x3f8b, |
|
0x2903, 0x65ff, 0x7f2b, 0x3f35, |
|
0xccf7, 0xf3f6, 0x438c, 0x3ec8, |
|
0x6b3d, 0xb876, 0x29e5, 0x3e3d, |
|
}; |
|
#endif |
|
#ifdef MIEEE |
|
static unsigned short P2[36] = { |
|
0x4009, 0xe6e8, 0xe793, 0xd574, |
|
0x401b, 0xa931, 0xc327, 0x780b, |
|
0x400f, 0x82ae, 0xf32b, 0xb0ac, |
|
0x3ff5, 0x541c, 0x18e7, 0x9a0c, |
|
0x3fc9, 0xca45, 0xf35e, 0x2651, |
|
0x3f89, 0x5650, 0x9069, 0x354d, |
|
0x3f33, 0xc3b2, 0xe8ce, 0x1812, |
|
0x3ec6, 0x4c29, 0x4c0c, 0x2234, |
|
0x3e3a, 0xccac, 0x3058, 0x8ff9, |
|
}; |
|
static unsigned short Q2[32] = { |
|
/*0x3ff0,0x0000,0x0000,0x0000,*/ |
|
0x4018, 0x18da, 0x568b, 0xe828, |
|
0x400d, 0x704d, 0xaad2, 0x3816, |
|
0x3ff6, 0x0847, 0x2aaf, 0x6200, |
|
0x3fcb, 0xada7, 0x5fd6, 0x3c9b, |
|
0x3f8b, 0x7c27, 0xadb6, 0xc78b, |
|
0x3f35, 0x7f2b, 0x65ff, 0x2903, |
|
0x3ec8, 0x438c, 0xf3f6, 0xccf7, |
|
0x3e3d, 0x29e5, 0xb876, 0x6b3d, |
|
}; |
|
#endif |
|
|
|
static double inverse_normal_func(double y0) |
|
{ |
|
double x, |
|
y, |
|
z, |
|
y2, |
|
x0, |
|
x1; |
|
int code; |
|
|
|
if (y0 <= 0.0) { |
|
mtherr("inverse_normal_func", DOMAIN); |
|
return (-DBL_MAX); |
} |
} |
return (y); |
if (y0 >= 1.0) { |
|
mtherr("inverse_normal_func", DOMAIN); |
|
return (DBL_MAX); |
|
} |
|
code = 1; |
|
y = y0; |
|
if (y > (1.0 - 0.13533528323661269189)) { /* 0.135... = exp(-2) */ |
|
y = 1.0 - y; |
|
code = 0; |
|
} |
|
if (y > 0.13533528323661269189) { |
|
y = y - 0.5; |
|
y2 = y * y; |
|
x = y + y * (y2 * polevl(y2, P0, 4) / p1evl(y2, Q0, 8)); |
|
x = x * s2pi; |
|
return (x); |
|
} |
|
x = sqrt(-2.0 * log(y)); |
|
x0 = x - log(x) / x; |
|
|
|
z = 1.0 / x; |
|
if (x < 8.0) /* y > exp(-32) = 1.2664165549e-14 */ |
|
x1 = z * polevl(z, P1, 8) / p1evl(z, Q1, 8); |
|
else |
|
x1 = z * polevl(z, P2, 8) / p1evl(z, Q2, 8); |
|
x = x0 - x1; |
|
if (code != 0) |
|
x = -x; |
|
return (x); |
} |
} |
|
/* ndtr.c |
|
* |
|
* Normal distribution function |
|
* |
|
* |
|
* |
|
* SYNOPSIS: |
|
* |
|
* double x, y, ndtr(); |
|
* |
|
* y = ndtr( x ); |
|
* |
|
* |
|
* |
|
* DESCRIPTION: |
|
* |
|
* Returns the area under the Gaussian probability density |
|
* function, integrated from minus infinity to x: |
|
* |
|
* x |
|
* - |
|
* 1 | | 2 |
|
* ndtr(x) = --------- | exp( - t /2 ) dt |
|
* sqrt(2pi) | | |
|
* - |
|
* -inf. |
|
* |
|
* = ( 1 + erf(z) ) / 2 |
|
* = erfc(z) / 2 |
|
* |
|
* where z = x/sqrt(2). Computation is via the functions |
|
* erf and erfc. |
|
* |
|
* |
|
* ACCURACY: |
|
* |
|
* Relative error: |
|
* arithmetic domain # trials peak rms |
|
* DEC -13,0 8000 2.1e-15 4.8e-16 |
|
* IEEE -13,0 30000 3.4e-14 6.7e-15 |
|
* |
|
* |
|
* ERROR MESSAGES: |
|
* |
|
* message condition value returned |
|
* erfc underflow x > 37.519379347 0.0 |
|
* |
|
*/ |
|
|
|
/* erf.c |
|
* |
|
* Error function |
|
* |
|
* |
|
* |
|
* SYNOPSIS: |
|
* |
|
* double x, y, erf(); |
|
* |
|
* y = erf( x ); |
|
* |
|
* |
|
* |
|
* DESCRIPTION: |
|
* |
|
* The integral is |
|
* |
|
* x |
|
* - |
|
* 2 | | 2 |
|
* erf(x) = -------- | exp( - t ) dt. |
|
* sqrt(pi) | | |
|
* - |
|
* 0 |
|
* |
|
* The magnitude of x is limited to 9.231948545 for DEC |
|
* arithmetic; 1 or -1 is returned outside this range. |
|
* |
|
* For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise |
|
* erf(x) = 1 - erfc(x). |
|
* |
|
* |
|
* |
|
* ACCURACY: |
|
* |
|
* Relative error: |
|
* arithmetic domain # trials peak rms |
|
* DEC 0,1 14000 4.7e-17 1.5e-17 |
|
* IEEE 0,1 30000 3.7e-16 1.0e-16 |
|
* |
|
*/ |
|
|
static double inverse_error_func(p) |
/* erfc.c |
double p; |
* |
|
* Complementary error function |
|
* |
|
* |
|
* |
|
* SYNOPSIS: |
|
* |
|
* double x, y, erfc(); |
|
* |
|
* y = erfc( x ); |
|
* |
|
* |
|
* |
|
* DESCRIPTION: |
|
* |
|
* |
|
* 1 - erf(x) = |
|
* |
|
* inf. |
|
* - |
|
* 2 | | 2 |
|
* erfc(x) = -------- | exp( - t ) dt |
|
* sqrt(pi) | | |
|
* - |
|
* x |
|
* |
|
* |
|
* For small x, erfc(x) = 1 - erf(x); otherwise rational |
|
* approximations are computed. |
|
* |
|
* |
|
* |
|
* ACCURACY: |
|
* |
|
* Relative error: |
|
* arithmetic domain # trials peak rms |
|
* DEC 0, 9.2319 12000 5.1e-16 1.2e-16 |
|
* IEEE 0,26.6417 30000 5.7e-14 1.5e-14 |
|
* |
|
* |
|
* ERROR MESSAGES: |
|
* |
|
* message condition value returned |
|
* erfc underflow x > 9.231948545 (DEC) 0.0 |
|
* |
|
* |
|
*/ |
|
|
|
/* |
|
Cephes Math Library Release 2.8: June, 2000 |
|
Copyright 1984, 1987, 1988, 1992, 2000 by Stephen L. Moshier |
|
*/ |
|
|
|
static const double SQRTH = 0.70710678118654752440084436210484903928483593768845; |
|
|
|
#ifdef UNK |
|
static double P[] = { |
|
2.46196981473530512524E-10, |
|
5.64189564831068821977E-1, |
|
7.46321056442269912687E0, |
|
4.86371970985681366614E1, |
|
1.96520832956077098242E2, |
|
5.26445194995477358631E2, |
|
9.34528527171957607540E2, |
|
1.02755188689515710272E3, |
|
5.57535335369399327526E2 |
|
}; |
|
static double Q[] = { |
|
/* 1.00000000000000000000E0,*/ |
|
1.32281951154744992508E1, |
|
8.67072140885989742329E1, |
|
3.54937778887819891062E2, |
|
9.75708501743205489753E2, |
|
1.82390916687909736289E3, |
|
2.24633760818710981792E3, |
|
1.65666309194161350182E3, |
|
5.57535340817727675546E2 |
|
}; |
|
static double R[] = { |
|
5.64189583547755073984E-1, |
|
1.27536670759978104416E0, |
|
5.01905042251180477414E0, |
|
6.16021097993053585195E0, |
|
7.40974269950448939160E0, |
|
2.97886665372100240670E0 |
|
}; |
|
static double S[] = { |
|
/* 1.00000000000000000000E0,*/ |
|
2.26052863220117276590E0, |
|
9.39603524938001434673E0, |
|
1.20489539808096656605E1, |
|
1.70814450747565897222E1, |
|
9.60896809063285878198E0, |
|
3.36907645100081516050E0 |
|
}; |
|
static double T[] = { |
|
9.60497373987051638749E0, |
|
9.00260197203842689217E1, |
|
2.23200534594684319226E3, |
|
7.00332514112805075473E3, |
|
5.55923013010394962768E4 |
|
}; |
|
static double U[] = { |
|
/* 1.00000000000000000000E0,*/ |
|
3.35617141647503099647E1, |
|
5.21357949780152679795E2, |
|
4.59432382970980127987E3, |
|
2.26290000613890934246E4, |
|
4.92673942608635921086E4 |
|
}; |
|
|
|
#define UTHRESH 37.519379347 |
|
#endif |
|
|
|
#ifdef DEC |
|
static unsigned short P[] = { |
|
0030207, 0054445, 0011173, 0021706, |
|
0040020, 0067272, 0030661, 0122075, |
|
0040756, 0151236, 0173053, 0067042, |
|
0041502, 0106175, 0062555, 0151457, |
|
0042104, 0102525, 0047401, 0003667, |
|
0042403, 0116176, 0011446, 0075303, |
|
0042551, 0120723, 0061641, 0123275, |
|
0042600, 0070651, 0007264, 0134516, |
|
0042413, 0061102, 0167507, 0176625 |
|
}; |
|
static unsigned short Q[] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0041123, 0123257, 0165741, 0017142, |
|
0041655, 0065027, 0173413, 0115450, |
|
0042261, 0074011, 0021573, 0004150, |
|
0042563, 0166530, 0013662, 0007200, |
|
0042743, 0176427, 0162443, 0105214, |
|
0043014, 0062546, 0153727, 0123772, |
|
0042717, 0012470, 0006227, 0067424, |
|
0042413, 0061103, 0003042, 0013254 |
|
}; |
|
static unsigned short R[] = { |
|
0040020, 0067272, 0101024, 0155421, |
|
0040243, 0037467, 0056706, 0026462, |
|
0040640, 0116017, 0120665, 0034315, |
|
0040705, 0020162, 0143350, 0060137, |
|
0040755, 0016234, 0134304, 0130157, |
|
0040476, 0122700, 0051070, 0015473 |
|
}; |
|
static unsigned short S[] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0040420, 0126200, 0044276, 0070413, |
|
0041026, 0053051, 0007302, 0063746, |
|
0041100, 0144203, 0174051, 0061151, |
|
0041210, 0123314, 0126343, 0177646, |
|
0041031, 0137125, 0051431, 0033011, |
|
0040527, 0117362, 0152661, 0066201 |
|
}; |
|
static unsigned short T[] = { |
|
0041031, 0126770, 0170672, 0166101, |
|
0041664, 0006522, 0072360, 0031770, |
|
0043013, 0100025, 0162641, 0126671, |
|
0043332, 0155231, 0161627, 0076200, |
|
0044131, 0024115, 0021020, 0117343 |
|
}; |
|
static unsigned short U[] = { |
|
/*0040200,0000000,0000000,0000000,*/ |
|
0041406, 0037461, 0177575, 0032714, |
|
0042402, 0053350, 0123061, 0153557, |
|
0043217, 0111227, 0032007, 0164217, |
|
0043660, 0145000, 0004013, 0160114, |
|
0044100, 0071544, 0167107, 0125471 |
|
}; |
|
#define UTHRESH 14.0 |
|
#endif |
|
|
|
#ifdef IBMPC |
|
static unsigned short P[] = { |
|
0x6479, 0xa24f, 0xeb24, 0x3df0, |
|
0x3488, 0x4636, 0x0dd7, 0x3fe2, |
|
0x6dc4, 0xdec5, 0xda53, 0x401d, |
|
0xba66, 0xacad, 0x518f, 0x4048, |
|
0x20f7, 0xa9e0, 0x90aa, 0x4068, |
|
0xcf58, 0xc264, 0x738f, 0x4080, |
|
0x34d8, 0x6c74, 0x343a, 0x408d, |
|
0x972a, 0x21d6, 0x0e35, 0x4090, |
|
0xffb3, 0x5de8, 0x6c48, 0x4081 |
|
}; |
|
static unsigned short Q[] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0x23cc, 0xfd7c, 0x74d5, 0x402a, |
|
0x7365, 0xfee1, 0xad42, 0x4055, |
|
0x610d, 0x246f, 0x2f01, 0x4076, |
|
0x41d0, 0x02f6, 0x7dab, 0x408e, |
|
0x7151, 0xfca4, 0x7fa2, 0x409c, |
|
0xf4ff, 0xdafa, 0x8cac, 0x40a1, |
|
0xede2, 0x0192, 0xe2a7, 0x4099, |
|
0x42d6, 0x60c4, 0x6c48, 0x4081 |
|
}; |
|
static unsigned short R[] = { |
|
0x9b62, 0x5042, 0x0dd7, 0x3fe2, |
|
0xc5a6, 0xebb8, 0x67e6, 0x3ff4, |
|
0xa71a, 0xf436, 0x1381, 0x4014, |
|
0x0c0c, 0x58dd, 0xa40e, 0x4018, |
|
0x960e, 0x9718, 0xa393, 0x401d, |
|
0x0367, 0x0a47, 0xd4b8, 0x4007 |
|
}; |
|
static unsigned short S[] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0xce21, 0x0917, 0x1590, 0x4002, |
|
0x4cfd, 0x21d8, 0xcac5, 0x4022, |
|
0x2c4d, 0x7f05, 0x1910, 0x4028, |
|
0x7ff5, 0x959c, 0x14d9, 0x4031, |
|
0x26c1, 0xaa63, 0x37ca, 0x4023, |
|
0x2d90, 0x5ab6, 0xf3de, 0x400a |
|
}; |
|
static unsigned short T[] = { |
|
0x5d88, 0x1e37, 0x35bf, 0x4023, |
|
0x067f, 0x4e9e, 0x81aa, 0x4056, |
|
0x35b7, 0xbcb4, 0x7002, 0x40a1, |
|
0xef90, 0x3c72, 0x5b53, 0x40bb, |
|
0x13dc, 0xa442, 0x2509, 0x40eb |
|
}; |
|
static unsigned short U[] = { |
|
/*0x0000,0x0000,0x0000,0x3ff0,*/ |
|
0xa6ba, 0x3fef, 0xc7e6, 0x4040, |
|
0x3aee, 0x14c6, 0x4add, 0x4080, |
|
0xfd12, 0xe680, 0xf252, 0x40b1, |
|
0x7c0a, 0x0101, 0x1940, 0x40d6, |
|
0xf567, 0x9dc8, 0x0e6c, 0x40e8 |
|
}; |
|
#define UTHRESH 37.519379347 |
|
#endif |
|
|
|
#ifdef MIEEE |
|
static unsigned short P[] = { |
|
0x3df0, 0xeb24, 0xa24f, 0x6479, |
|
0x3fe2, 0x0dd7, 0x4636, 0x3488, |
|
0x401d, 0xda53, 0xdec5, 0x6dc4, |
|
0x4048, 0x518f, 0xacad, 0xba66, |
|
0x4068, 0x90aa, 0xa9e0, 0x20f7, |
|
0x4080, 0x738f, 0xc264, 0xcf58, |
|
0x408d, 0x343a, 0x6c74, 0x34d8, |
|
0x4090, 0x0e35, 0x21d6, 0x972a, |
|
0x4081, 0x6c48, 0x5de8, 0xffb3 |
|
}; |
|
static unsigned short Q[] = { |
|
0x402a, 0x74d5, 0xfd7c, 0x23cc, |
|
0x4055, 0xad42, 0xfee1, 0x7365, |
|
0x4076, 0x2f01, 0x246f, 0x610d, |
|
0x408e, 0x7dab, 0x02f6, 0x41d0, |
|
0x409c, 0x7fa2, 0xfca4, 0x7151, |
|
0x40a1, 0x8cac, 0xdafa, 0xf4ff, |
|
0x4099, 0xe2a7, 0x0192, 0xede2, |
|
0x4081, 0x6c48, 0x60c4, 0x42d6 |
|
}; |
|
static unsigned short R[] = { |
|
0x3fe2, 0x0dd7, 0x5042, 0x9b62, |
|
0x3ff4, 0x67e6, 0xebb8, 0xc5a6, |
|
0x4014, 0x1381, 0xf436, 0xa71a, |
|
0x4018, 0xa40e, 0x58dd, 0x0c0c, |
|
0x401d, 0xa393, 0x9718, 0x960e, |
|
0x4007, 0xd4b8, 0x0a47, 0x0367 |
|
}; |
|
static unsigned short S[] = { |
|
0x4002, 0x1590, 0x0917, 0xce21, |
|
0x4022, 0xcac5, 0x21d8, 0x4cfd, |
|
0x4028, 0x1910, 0x7f05, 0x2c4d, |
|
0x4031, 0x14d9, 0x959c, 0x7ff5, |
|
0x4023, 0x37ca, 0xaa63, 0x26c1, |
|
0x400a, 0xf3de, 0x5ab6, 0x2d90 |
|
}; |
|
static unsigned short T[] = { |
|
0x4023, 0x35bf, 0x1e37, 0x5d88, |
|
0x4056, 0x81aa, 0x4e9e, 0x067f, |
|
0x40a1, 0x7002, 0xbcb4, 0x35b7, |
|
0x40bb, 0x5b53, 0x3c72, 0xef90, |
|
0x40eb, 0x2509, 0xa442, 0x13dc |
|
}; |
|
static unsigned short U[] = { |
|
0x4040, 0xc7e6, 0x3fef, 0xa6ba, |
|
0x4080, 0x4add, 0x14c6, 0x3aee, |
|
0x40b1, 0xf252, 0xe680, 0xfd12, |
|
0x40d6, 0x1940, 0x0101, 0x7c0a, |
|
0x40e8, 0x0e6c, 0x9dc8, 0xf567 |
|
}; |
|
#define UTHRESH 37.519379347 |
|
#endif |
|
/* |
|
double ndtr(double a) |
{ |
{ |
/* |
double x, |
Source: This routine was derived (using f2c) from the |
y, |
FORTRAN subroutine MERFI found in |
z; |
ACM Algorithm 602 obtained from netlib. |
|
|
|
MDNRIS code contains the 1978 Copyright |
x = a * SQRTH; |
by IMSL, INC. . Since MERFI has been |
z = fabs(x); |
submitted to netlib, it may be used with |
|
the restriction that it may only be |
|
used for noncommercial purposes and that |
|
IMSL be acknowledged as the copyright-holder |
|
of the code. |
|
*/ |
|
|
|
|
if (z < SQRTH) |
|
y = 0.5 + 0.5 * erf(x); |
|
|
|
else { |
|
y = 0.5 * erfc(z); |
|
|
/* Initialized data */ |
if (x > 0) |
static double a1 = -.5751703; |
y = 1.0 - y; |
static double a2 = -1.896513; |
} |
static double a3 = -.05496261; |
|
static double b0 = -.113773; |
|
static double b1 = -3.293474; |
|
static double b2 = -2.374996; |
|
static double b3 = -1.187515; |
|
static double c0 = -.1146666; |
|
static double c1 = -.1314774; |
|
static double c2 = -.2368201; |
|
static double c3 = .05073975; |
|
static double d0 = -44.27977; |
|
static double d1 = 21.98546; |
|
static double d2 = -7.586103; |
|
static double e0 = -.05668422; |
|
static double e1 = .3937021; |
|
static double e2 = -.3166501; |
|
static double e3 = .06208963; |
|
static double f0 = -6.266786; |
|
static double f1 = 4.666263; |
|
static double f2 = -2.962883; |
|
static double g0 = 1.851159e-4; |
|
static double g1 = -.002028152; |
|
static double g2 = -.1498384; |
|
static double g3 = .01078639; |
|
static double h0 = .09952975; |
|
static double h1 = .5211733; |
|
static double h2 = -.06888301; |
|
|
|
/* Local variables */ |
return (y); |
static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn; |
} |
|
*/ |
|
double erf(double); |
|
|
x = p; |
double erfc(double a) |
|
{ |
|
double p, |
|
q, |
|
x, |
|
y, |
|
z; |
|
|
/* determine sign of x */ |
if (a < 0.0) |
if (x > 0) |
x = -a; |
sigma = 1.0; |
|
else |
else |
sigma = -1.0; |
x = a; |
|
|
/* Note: -1.0 < x < 1.0 */ |
if (x < 1.0) |
|
return (1.0 - erf(a)); |
|
|
z = fabs(x); |
z = -a * a; |
|
|
/* z between 0.0 and 0.85, approx. f by a |
if (z < DBL_MIN_10_EXP) { |
rational function in z */ |
under: |
|
mtherr("erfc", UNDERFLOW); |
|
if (a < 0) |
|
return (2.0); |
|
else |
|
return (0.0); |
|
} |
|
z = exp(z); |
|
|
if (z <= 0.85) { |
if (x < 8.0) { |
z2 = z * z; |
p = polevl(x, P, 8); |
f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2 |
q = p1evl(x, Q, 8); |
/ (b2 + z2 + a3 / (b3 + z2)))); |
|
|
|
/* z greater than 0.85 */ |
|
} else { |
} else { |
a = 1.0 - z; |
p = polevl(x, R, 5); |
b = z; |
q = p1evl(x, S, 6); |
|
} |
|
y = (z * p) / q; |
|
|
/* reduced argument is in (0.85,1.0), |
if (a < 0) |
obtain the transformed variable */ |
y = 2.0 - y; |
|
|
w = sqrt(-(double) log(a + a * b)); |
if (y == 0.0) |
|
goto under; |
|
|
/* w greater than 4.0, approx. f by a |
return (y); |
rational function in 1.0 / w */ |
} |
|
|
if (w >= 4.0) { |
double erf(double x) |
wi = 1.0 / w; |
{ |
sn = ((g3 * wi + g2) * wi + g1) * wi; |
double y, |
sd = ((wi + h2) * wi + h1) * wi + h0; |
z; |
f = w + w * (g0 + sn / sd); |
|
|
|
/* w between 2.5 and 4.0, approx. |
if (fabs(x) > 1.0) |
f by a rational function in w */ |
return (1.0 - erfc(x)); |
|
z = x * x; |
|
y = x * polevl(z, T, 4) / p1evl(z, U, 5); |
|
return (y); |
|
|
} else if (w < 4.0 && w > 2.5) { |
} |
sn = ((e3 * w + e2) * w + e1) * w; |
|
sd = ((w + f2) * w + f1) * w + f0; |
|
f = w + w * (e0 + sn / sd); |
|
|
|
/* w between 1.13222 and 2.5, approx. f by |
static double inverse_error_func(double y) |
a rational function in w */ |
{ |
} else if (w <= 2.5 && w > 1.13222) { |
double x = 0.0; /* The output */ |
sn = ((c3 * w + c2) * w + c1) * w; |
double z = 0.0; /* Intermadiate variable */ |
sd = ((w + d2) * w + d1) * w + d0; |
double y0 = 0.7; /* Central range variable */ |
f = w + w * (c0 + sn / sd); |
|
} |
/* Coefficients in rational approximations. */ |
|
double a[4] = {0.886226899, -1.645349621, 0.914624893, -0.140543331}; |
|
double b[4] = {-2.118377725, 1.442710462, -0.329097515, 0.012229801}; |
|
double c[4] = {-1.970840454, -1.624906493, 3.429567803, 1.641345311}; |
|
double d[2] = {3.543889200, 1.637067800}; |
|
|
|
if ((y < -1.0) || (1.0 < y)) { |
|
printf("inverse_error_func: The value out of the range of the function"); |
|
x = log(-1.0); |
|
} else if ((y == -1.0) || (1.0 == y)) { |
|
x = -y * log(0.0); |
|
} else if ((-1.0 < y) && (y < -y0)) { |
|
z = sqrt(-log((1.0 + y) / 2.0)); |
|
x = -(((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0); |
|
} else { |
|
if ((-y0 < y) && (y < y0)) { |
|
z = y * y; |
|
x = y * (((a[3] * z + a[2]) * z + a[1]) * z + a[0]) / |
|
((((b[3] * z + b[3]) * z + b[1]) * z + b[0]) * z + 1.0); |
|
} else if ((y0 < y) && (y < 1.0)) { |
|
z = sqrt(-log((1.0 - y) / 2.0)); |
|
x = (((c[3] * z + c[2]) * z + c[1]) * z + c[0]) / ((d[1] * z + d[0]) * z + 1.0); |
|
} |
|
/* Two steps of Newton-Raphson correction to full accuracy. */ |
|
x = x - (erf(x) - y) / (2.0 / sqrt(PI) * exp(-x * x)); |
|
x = x - (erf(x) - y) / (2.0 / sqrt(PI) * exp(-x * x)); |
} |
} |
y = sigma * f; |
return (x); |
return (y); |
|
} |
} |