Math-CDF
view release on metacpan or search on metacpan
cdflib/dcdflib.c view on Meta::CPAN
i99999 = 1;
goto S270;
S10:
fb = *fx;
*xlo = *xhi;
a = *x = *xlo;
/*
GET-FUNCTION-VALUE
*/
i99999 = 2;
goto S270;
S20:
/*
Check that F(ZXLO) < 0 < F(ZXHI) or
F(ZXLO) > 0 > F(ZXHI)
*/
if(!(fb < 0.0e0)) goto S40;
if(!(*fx < 0.0e0)) goto S30;
*status = -1;
*qleft = *fx < fb;
*qhi = 0;
return;
S40:
S30:
if(!(fb > 0.0e0)) goto S60;
if(!(*fx > 0.0e0)) goto S50;
*status = -1;
*qleft = *fx > fb;
*qhi = 1;
return;
S60:
S50:
fa = *fx;
first = 1;
S70:
c = a;
fc = fa;
ext = 0;
S80:
if(!(fabs(fc) < fabs(fb))) goto S100;
if(!(c != a)) goto S90;
d = a;
fd = fa;
S90:
a = b;
fa = fb;
*xlo = c;
b = *xlo;
fb = fc;
c = a;
fc = fa;
S100:
tol = ftol(*xlo);
m = (c+b)*.5e0;
mb = m-b;
if(!(fabs(mb) > tol)) goto S240;
if(!(ext > 3)) goto S110;
w = mb;
goto S190;
S110:
tol = fifdsign(tol,mb);
p = (b-a)*fb;
if(!first) goto S120;
q = fa-fb;
first = 0;
goto S130;
S120:
fdb = (fd-fb)/(d-b);
fda = (fd-fa)/(d-a);
p = fda*p;
q = fdb*fa-fda*fb;
S130:
if(!(p < 0.0e0)) goto S140;
p = -p;
q = -q;
S140:
if(ext == 3) p *= 2.0e0;
if(!(p*1.0e0 == 0.0e0 || p <= q*tol)) goto S150;
w = tol;
goto S180;
S150:
if(!(p < mb*q)) goto S160;
w = p/q;
goto S170;
S160:
w = mb;
S190:
S180:
S170:
d = a;
fd = fa;
a = b;
fa = fb;
b += w;
*xlo = b;
*x = *xlo;
/*
GET-FUNCTION-VALUE
*/
i99999 = 3;
goto S270;
S200:
fb = *fx;
if(!(fc*fb >= 0.0e0)) goto S210;
goto S70;
S210:
if(!(w == mb)) goto S220;
ext = 0;
goto S230;
S220:
ext += 1;
S230:
goto S80;
S240:
*xhi = c;
qrzero = fc >= 0.0e0 && fb <= 0.0e0 || fc < 0.0e0 && fb >= 0.0e0;
if(!qrzero) goto S250;
*status = 0;
goto S260;
S250:
*status = -1;
cdflib/dcdflib.c view on Meta::CPAN
*/
{
static double c = .564189583547756e0;
static double a[5] = {
.771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
.479137145607681e-01,.128379167095513e+00
};
static double b[3] = {
.301048631703895e-02,.538971687740286e-01,.375795757275549e+00
};
static double p[8] = {
-1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
4.51918953711873e+02,3.00459261020162e+02
};
static double q[8] = {
1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7.90950925327898e+02,3.00459260956983e+02
};
static double r[5] = {
2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
4.65807828718470e+00,2.82094791773523e-01
};
static double s[4] = {
9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
1.80124575948747e+01
};
static double erf1,ax,bot,t,top,x2;
/*
..
.. Executable Statements ..
*/
ax = fabs(*x);
if(ax > 0.5e0) goto S10;
t = *x**x;
top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
erf1 = *x*(top/bot);
return erf1;
S10:
if(ax > 4.0e0) goto S20;
top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
7];
bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
7];
erf1 = 0.5e0+(0.5e0-exp(-(*x**x))*top/bot);
if(*x < 0.0e0) erf1 = -erf1;
return erf1;
S20:
if(ax >= 5.8e0) goto S30;
x2 = *x**x;
t = 1.0e0/x2;
top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
erf1 = (c-top/(x2*bot))/ax;
erf1 = 0.5e0+(0.5e0-exp(-x2)*erf1);
if(*x < 0.0e0) erf1 = -erf1;
return erf1;
S30:
erf1 = fifdsign(1.0e0,*x);
return erf1;
}
double erfc1(int *ind,double *x)
/*
-----------------------------------------------------------------------
EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
ERFC1(IND,X) = ERFC(X) IF IND = 0
ERFC1(IND,X) = EXP(X*X)*ERFC(X) OTHERWISE
-----------------------------------------------------------------------
*/
{
static double c = .564189583547756e0;
static double a[5] = {
.771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
.479137145607681e-01,.128379167095513e+00
};
static double b[3] = {
.301048631703895e-02,.538971687740286e-01,.375795757275549e+00
};
static double p[8] = {
-1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
4.51918953711873e+02,3.00459261020162e+02
};
static double q[8] = {
1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
7.90950925327898e+02,3.00459260956983e+02
};
static double r[5] = {
2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
4.65807828718470e+00,2.82094791773523e-01
};
static double s[4] = {
9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
1.80124575948747e+01
};
static int K1 = 1;
static double erfc1,ax,bot,e,t,top,w;
/*
..
.. Executable Statements ..
*/
/*
ABS(X) .LE. 0.5
*/
ax = fabs(*x);
if(ax > 0.5e0) goto S10;
t = *x**x;
top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
if(*ind != 0) erfc1 = exp(t)*erfc1;
return erfc1;
S10:
/*
0.5 .LT. ABS(X) .LE. 4
*/
if(ax > 4.0e0) goto S20;
cdflib/dcdflib.c view on Meta::CPAN
-0.322232431088e0,-1.000000000000e0,-0.342242088547e0,-0.204231210245e-1,
-0.453642210148e-4
};
static int K1 = 5;
static double stvaln,sign,y,z;
/*
..
.. Executable Statements ..
*/
if(!(*p <= 0.5e0)) goto S10;
sign = -1.0e0;
z = *p;
goto S20;
S10:
sign = 1.0e0;
z = 1.0e0-*p;
S20:
y = sqrt(-(2.0e0*log(z)));
stvaln = y+devlpl(xnum,&K1,&y)/devlpl(xden,&K1,&y);
stvaln = sign*stvaln;
return stvaln;
}
/************************************************************************
FIFDINT:
Truncates a double precision number to an integer and returns the
value in a double.
************************************************************************/
double fifdint(double a)
/* a - number to be truncated */
{
long temp;
temp = (long)(a);
return (double)(temp);
}
/************************************************************************
FIFDMAX1:
returns the maximum of two numbers a and b
************************************************************************/
double fifdmax1(double a,double b)
/* a - first number */
/* b - second number */
{
if (a < b) return b;
else return a;
}
/************************************************************************
FIFDMIN1:
returns the minimum of two numbers a and b
************************************************************************/
double fifdmin1(double a,double b)
/* a - first number */
/* b - second number */
{
if (a < b) return a;
else return b;
}
/************************************************************************
FIFDSIGN:
transfers the sign of the variable "sign" to the variable "mag"
************************************************************************/
double fifdsign(double mag,double sign)
/* mag - magnitude */
/* sign - sign to be transfered */
{
if (mag < 0) mag = -mag;
if (sign < 0) mag = -mag;
return mag;
}
/************************************************************************
FIFIDINT:
Truncates a double precision number to a long integer
************************************************************************/
long fifidint(double a)
/* a - number to be truncated */
{
return (long)(a);
}
/************************************************************************
FIFMOD:
returns the modulo of a and b
************************************************************************/
long fifmod(long a,long b)
/* a - numerator */
/* b - denominator */
{
return a % b;
}
/************************************************************************
FTNSTOP:
Prints msg to standard error and then exits
************************************************************************/
void ftnstop(char* msg)
/* msg - error message */
{
if (msg != NULL) fprintf(stderr,"%s\n",msg);
exit(EXIT_FAILURE); /* EXIT_FAILURE from stdlib.h, or use an int */
}
( run in 0.536 second using v1.01-cache-2.11-cpan-71847e10f99 )