Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_63.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void grat1(double *a,double *x,double *r,double *p,double *q,
00003 double *eps)
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
00010
00011
00012
00013
00014
00015
00016
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
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
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
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 }