Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_10.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void bratio(double *a,double *b,double *x,double *y,double *w,
00003 double *w1,int *ierr)
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 {
00039 static int K1 = 1;
00040 static double a0,b0,eps,lambda,t,x0,y0,z;
00041 static int ierr1,ind,n;
00042 static double T2,T3,T4,T5;
00043
00044
00045
00046
00047
00048
00049
00050
00051 eps = spmpar(&K1);
00052 *w = *w1 = 0.0e0;
00053 if(*a < 0.0e0 || *b < 0.0e0) goto S270;
00054 if(*a == 0.0e0 && *b == 0.0e0) goto S280;
00055 if(*x < 0.0e0 || *x > 1.0e0) goto S290;
00056 if(*y < 0.0e0 || *y > 1.0e0) goto S300;
00057 z = *x+*y-0.5e0-0.5e0;
00058 if(fabs(z) > 3.0e0*eps) goto S310;
00059 *ierr = 0;
00060 if(*x == 0.0e0) goto S210;
00061 if(*y == 0.0e0) goto S230;
00062 if(*a == 0.0e0) goto S240;
00063 if(*b == 0.0e0) goto S220;
00064 eps = fifdmax1(eps,1.e-15);
00065 if(fifdmax1(*a,*b) < 1.e-3*eps) goto S260;
00066 ind = 0;
00067 a0 = *a;
00068 b0 = *b;
00069 x0 = *x;
00070 y0 = *y;
00071 if(fifdmin1(a0,b0) > 1.0e0) goto S40;
00072
00073
00074
00075 if(*x <= 0.5e0) goto S10;
00076 ind = 1;
00077 a0 = *b;
00078 b0 = *a;
00079 x0 = *y;
00080 y0 = *x;
00081 S10:
00082 if(b0 < fifdmin1(eps,eps*a0)) goto S90;
00083 if(a0 < fifdmin1(eps,eps*b0) && b0*x0 <= 1.0e0) goto S100;
00084 if(fifdmax1(a0,b0) > 1.0e0) goto S20;
00085 if(a0 >= fifdmin1(0.2e0,b0)) goto S110;
00086 if(pow(x0,a0) <= 0.9e0) goto S110;
00087 if(x0 >= 0.3e0) goto S120;
00088 n = 20;
00089 goto S140;
00090 S20:
00091 if(b0 <= 1.0e0) goto S110;
00092 if(x0 >= 0.3e0) goto S120;
00093 if(x0 >= 0.1e0) goto S30;
00094 if(pow(x0*b0,a0) <= 0.7e0) goto S110;
00095 S30:
00096 if(b0 > 15.0e0) goto S150;
00097 n = 20;
00098 goto S140;
00099 S40:
00100
00101
00102
00103 if(*a > *b) goto S50;
00104 lambda = *a-(*a+*b)**x;
00105 goto S60;
00106 S50:
00107 lambda = (*a+*b)**y-*b;
00108 S60:
00109 if(lambda >= 0.0e0) goto S70;
00110 ind = 1;
00111 a0 = *b;
00112 b0 = *a;
00113 x0 = *y;
00114 y0 = *x;
00115 lambda = fabs(lambda);
00116 S70:
00117 if(b0 < 40.0e0 && b0*x0 <= 0.7e0) goto S110;
00118 if(b0 < 40.0e0) goto S160;
00119 if(a0 > b0) goto S80;
00120 if(a0 <= 100.0e0) goto S130;
00121 if(lambda > 0.03e0*a0) goto S130;
00122 goto S200;
00123 S80:
00124 if(b0 <= 100.0e0) goto S130;
00125 if(lambda > 0.03e0*b0) goto S130;
00126 goto S200;
00127 S90:
00128
00129
00130
00131 *w = fpser(&a0,&b0,&x0,&eps);
00132 *w1 = 0.5e0+(0.5e0-*w);
00133 goto S250;
00134 S100:
00135 *w1 = apser(&a0,&b0,&x0,&eps);
00136 *w = 0.5e0+(0.5e0-*w1);
00137 goto S250;
00138 S110:
00139 *w = bpser(&a0,&b0,&x0,&eps);
00140 *w1 = 0.5e0+(0.5e0-*w);
00141 goto S250;
00142 S120:
00143 *w1 = bpser(&b0,&a0,&y0,&eps);
00144 *w = 0.5e0+(0.5e0-*w1);
00145 goto S250;
00146 S130:
00147 T2 = 15.0e0*eps;
00148 *w = bfrac(&a0,&b0,&x0,&y0,&lambda,&T2);
00149 *w1 = 0.5e0+(0.5e0-*w);
00150 goto S250;
00151 S140:
00152 *w1 = bup(&b0,&a0,&y0,&x0,&n,&eps);
00153 b0 += (double)n;
00154 S150:
00155 T3 = 15.0e0*eps;
00156 bgrat(&b0,&a0,&y0,&x0,w1,&T3,&ierr1);
00157 *w = 0.5e0+(0.5e0-*w1);
00158 goto S250;
00159 S160:
00160 n = b0;
00161 b0 -= (double)n;
00162 if(b0 != 0.0e0) goto S170;
00163 n -= 1;
00164 b0 = 1.0e0;
00165 S170:
00166 *w = bup(&b0,&a0,&y0,&x0,&n,&eps);
00167 if(x0 > 0.7e0) goto S180;
00168 *w += bpser(&a0,&b0,&x0,&eps);
00169 *w1 = 0.5e0+(0.5e0-*w);
00170 goto S250;
00171 S180:
00172 if(a0 > 15.0e0) goto S190;
00173 n = 20;
00174 *w += bup(&a0,&b0,&x0,&y0,&n,&eps);
00175 a0 += (double)n;
00176 S190:
00177 T4 = 15.0e0*eps;
00178 bgrat(&a0,&b0,&x0,&y0,w,&T4,&ierr1);
00179 *w1 = 0.5e0+(0.5e0-*w);
00180 goto S250;
00181 S200:
00182 T5 = 100.0e0*eps;
00183 *w = basym(&a0,&b0,&lambda,&T5);
00184 *w1 = 0.5e0+(0.5e0-*w);
00185 goto S250;
00186 S210:
00187
00188
00189
00190 if(*a == 0.0e0) goto S320;
00191 S220:
00192 *w = 0.0e0;
00193 *w1 = 1.0e0;
00194 return;
00195 S230:
00196 if(*b == 0.0e0) goto S330;
00197 S240:
00198 *w = 1.0e0;
00199 *w1 = 0.0e0;
00200 return;
00201 S250:
00202 if(ind == 0) return;
00203 t = *w;
00204 *w = *w1;
00205 *w1 = t;
00206 return;
00207 S260:
00208
00209
00210
00211 *w = *b/(*a+*b);
00212 *w1 = *a/(*a+*b);
00213 return;
00214 S270:
00215
00216
00217
00218 *ierr = 1;
00219 return;
00220 S280:
00221 *ierr = 2;
00222 return;
00223 S290:
00224 *ierr = 3;
00225 return;
00226 S300:
00227 *ierr = 4;
00228 return;
00229 S310:
00230 *ierr = 5;
00231 return;
00232 S320:
00233 *ierr = 6;
00234 return;
00235 S330:
00236 *ierr = 7;
00237 return;
00238 }