Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_20.c
Go to the documentation of this file.00001 #include "cdflib.h"
00002 void cdfgam(int *which,double *p,double *q,double *x,double *shape,
00003 double *scale,int *status,double *bound)
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
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 {
00095 #define tol (1.0e-8)
00096 #define atol (1.0e-50)
00097 #define zero (1.0e-300)
00098 #define inf 1.0e300
00099 static int K1 = 1;
00100 static double K5 = 0.5e0;
00101 static double K6 = 5.0e0;
00102 static double xx,fx,xscale,cum,ccum,pq,porq;
00103 static int ierr;
00104 static unsigned long qhi,qleft,qporq;
00105 static double T2,T3,T4,T7,T8,T9;
00106
00107
00108
00109
00110
00111
00112
00113 if(!(*which < 1 || *which > 4)) goto S30;
00114 if(!(*which < 1)) goto S10;
00115 *bound = 1.0e0;
00116 goto S20;
00117 S10:
00118 *bound = 4.0e0;
00119 S20:
00120 *status = -1;
00121 return;
00122 S30:
00123 if(*which == 1) goto S70;
00124
00125
00126
00127 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00128 if(!(*p < 0.0e0)) goto S40;
00129 *bound = 0.0e0;
00130 goto S50;
00131 S40:
00132 *bound = 1.0e0;
00133 S50:
00134 *status = -2;
00135 return;
00136 S70:
00137 S60:
00138 if(*which == 1) goto S110;
00139
00140
00141
00142 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00143 if(!(*q <= 0.0e0)) goto S80;
00144 *bound = 0.0e0;
00145 goto S90;
00146 S80:
00147 *bound = 1.0e0;
00148 S90:
00149 *status = -3;
00150 return;
00151 S110:
00152 S100:
00153 if(*which == 2) goto S130;
00154
00155
00156
00157 if(!(*x < 0.0e0)) goto S120;
00158 *bound = 0.0e0;
00159 *status = -4;
00160 return;
00161 S130:
00162 S120:
00163 if(*which == 3) goto S150;
00164
00165
00166
00167 if(!(*shape <= 0.0e0)) goto S140;
00168 *bound = 0.0e0;
00169 *status = -5;
00170 return;
00171 S150:
00172 S140:
00173 if(*which == 4) goto S170;
00174
00175
00176
00177 if(!(*scale <= 0.0e0)) goto S160;
00178 *bound = 0.0e0;
00179 *status = -6;
00180 return;
00181 S170:
00182 S160:
00183 if(*which == 1) goto S210;
00184
00185
00186
00187 pq = *p+*q;
00188 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
00189 if(!(pq < 0.0e0)) goto S180;
00190 *bound = 0.0e0;
00191 goto S190;
00192 S180:
00193 *bound = 1.0e0;
00194 S190:
00195 *status = 3;
00196 return;
00197 S210:
00198 S200:
00199 if(*which == 1) goto S240;
00200
00201
00202
00203 qporq = *p <= *q;
00204 if(!qporq) goto S220;
00205 porq = *p;
00206 goto S230;
00207 S220:
00208 porq = *q;
00209 S240:
00210 S230:
00211
00212
00213
00214 if(1 == *which) {
00215
00216
00217
00218 *status = 0;
00219 xscale = *x**scale;
00220 cumgam(&xscale,shape,p,q);
00221 if(porq > 1.5e0) *status = 10;
00222 }
00223 else if(2 == *which) {
00224
00225
00226
00227 T2 = -1.0e0;
00228 gaminv(shape,&xx,&T2,p,q,&ierr);
00229 if(ierr < 0.0e0) {
00230 *status = 10;
00231 return;
00232 }
00233 else {
00234 *x = xx/ *scale;
00235 *status = 0;
00236 }
00237 }
00238 else if(3 == *which) {
00239
00240
00241
00242 *shape = 5.0e0;
00243 xscale = *x**scale;
00244 T3 = zero;
00245 T4 = inf;
00246 T7 = atol;
00247 T8 = tol;
00248 dstinv(&T3,&T4,&K5,&K5,&K6,&T7,&T8);
00249 *status = 0;
00250 dinvr(status,shape,&fx,&qleft,&qhi);
00251 S250:
00252 if(!(*status == 1)) goto S290;
00253 cumgam(&xscale,shape,&cum,&ccum);
00254 if(!qporq) goto S260;
00255 fx = cum-*p;
00256 goto S270;
00257 S260:
00258 fx = ccum-*q;
00259 S270:
00260 if(!(qporq && cum > 1.5e0 || !qporq && ccum > 1.5e0)) goto S280;
00261 *status = 10;
00262 return;
00263 S280:
00264 dinvr(status,shape,&fx,&qleft,&qhi);
00265 goto S250;
00266 S290:
00267 if(!(*status == -1)) goto S320;
00268 if(!qleft) goto S300;
00269 *status = 1;
00270 *bound = zero;
00271 goto S310;
00272 S300:
00273 *status = 2;
00274 *bound = inf;
00275 S320:
00276 S310:
00277 ;
00278 }
00279 else if(4 == *which) {
00280
00281
00282
00283 T9 = -1.0e0;
00284 gaminv(shape,&xx,&T9,p,q,&ierr);
00285 if(ierr < 0.0e0) {
00286 *status = 10;
00287 return;
00288 }
00289 else {
00290 *scale = xx/ *x;
00291 *status = 0;
00292 }
00293 }
00294 return;
00295 #undef tol
00296 #undef atol
00297 #undef zero
00298 #undef inf
00299 }