Doxygen Source Code Documentation
cdf_63.c File Reference
#include "cdflib.h"Go to the source code of this file.
Functions | |
| void | grat1 (double *a, double *x, double *r, double *p, double *q, double *eps) |
Function Documentation
|
||||||||||||||||||||||||||||
|
Definition at line 2 of file cdf_63.c. References a, c, erf1(), erfc1(), gam1(), l, p, q, r, and rexp(). Referenced by bgrat().
00004 {
00005 static int K2 = 0;
00006 static double a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,t,tol,w,z,T1,T3;
00007 /*
00008 ..
00009 .. Executable Statements ..
00010 */
00011 /*
00012 -----------------------------------------------------------------------
00013 EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
00014 P(A,X) AND Q(A,X)
00015 IT IS ASSUMED THAT A .LE. 1. EPS IS THE TOLERANCE TO BE USED.
00016 THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
00017 -----------------------------------------------------------------------
00018 */
00019 if(*a**x == 0.0e0) goto S120;
00020 if(*a == 0.5e0) goto S100;
00021 if(*x < 1.1e0) goto S10;
00022 goto S60;
00023 S10:
00024 /*
00025 TAYLOR SERIES FOR P(A,X)/X**A
00026 */
00027 an = 3.0e0;
00028 c = *x;
00029 sum = *x/(*a+3.0e0);
00030 tol = 0.1e0**eps/(*a+1.0e0);
00031 S20:
00032 an += 1.0e0;
00033 c = -(c*(*x/an));
00034 t = c/(*a+an);
00035 sum += t;
00036 if(fabs(t) > tol) goto S20;
00037 j = *a**x*((sum/6.0e0-0.5e0/(*a+2.0e0))**x+1.0e0/(*a+1.0e0));
00038 z = *a*log(*x);
00039 h = gam1(a);
00040 g = 1.0e0+h;
00041 if(*x < 0.25e0) goto S30;
00042 if(*a < *x/2.59e0) goto S50;
00043 goto S40;
00044 S30:
00045 if(z > -.13394e0) goto S50;
00046 S40:
00047 w = exp(z);
00048 *p = w*g*(0.5e0+(0.5e0-j));
00049 *q = 0.5e0+(0.5e0-*p);
00050 return;
00051 S50:
00052 l = rexp(&z);
00053 w = 0.5e0+(0.5e0+l);
00054 *q = (w*j-l)*g-h;
00055 if(*q < 0.0e0) goto S90;
00056 *p = 0.5e0+(0.5e0-*q);
00057 return;
00058 S60:
00059 /*
00060 CONTINUED FRACTION EXPANSION
00061 */
00062 a2nm1 = a2n = 1.0e0;
00063 b2nm1 = *x;
00064 b2n = *x+(1.0e0-*a);
00065 c = 1.0e0;
00066 S70:
00067 a2nm1 = *x*a2n+c*a2nm1;
00068 b2nm1 = *x*b2n+c*b2nm1;
00069 am0 = a2nm1/b2nm1;
00070 c += 1.0e0;
00071 cma = c-*a;
00072 a2n = a2nm1+cma*a2n;
00073 b2n = b2nm1+cma*b2n;
00074 an0 = a2n/b2n;
00075 if(fabs(an0-am0) >= *eps*an0) goto S70;
00076 *q = *r*an0;
00077 *p = 0.5e0+(0.5e0-*q);
00078 return;
00079 S80:
00080 /*
00081 SPECIAL CASES
00082 */
00083 *p = 0.0e0;
00084 *q = 1.0e0;
00085 return;
00086 S90:
00087 *p = 1.0e0;
00088 *q = 0.0e0;
00089 return;
00090 S100:
00091 if(*x >= 0.25e0) goto S110;
00092 T1 = sqrt(*x);
00093 *p = erf1(&T1);
00094 *q = 0.5e0+(0.5e0-*p);
00095 return;
00096 S110:
00097 T3 = sqrt(*x);
00098 *q = erfc1(&K2,&T3);
00099 *p = 0.5e0+(0.5e0-*q);
00100 return;
00101 S120:
00102 if(*x <= *a) goto S80;
00103 goto S90;
00104 } /* END */
|
