00001 #pragma ident "$Id: SpecialFunctions.cpp 3140 2012-06-18 15:03:02Z susancummins $"
00002
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "SpecialFunctions.hpp"
00033
00034
00035 using namespace std;
00036
00037 namespace gpstk
00038 {
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 double gamma(double val)
00051 {
00052
00053 double inf( 9.0e+99 );
00054
00055 if( val == 0.0 )
00056 {
00057 return inf;
00058 }
00059
00060 if ( ( val < 0.0 ) &&
00061 ( floor(val) == val ) )
00062 {
00063 return inf;
00064 }
00065
00066
00067 int g(7);
00068 const double lanczos_coef[] = { 0.99999999999980993,
00069 676.5203681218851,
00070 -1259.1392167224028,
00071 771.32342877765313,
00072 -176.61502916214059,
00073 12.507343278686905,
00074 -0.13857109526572012,
00075 9.9843695780195716e-6,
00076 1.5056327351493116e-7 };
00077
00078 if(val < 0.5)
00079 {
00080 return ( PI / (sin(PI*val)*gamma(1.0-val)) );
00081 }
00082 else
00083 {
00084 val -= 1.0;
00085
00086 double x( lanczos_coef[0] );
00087
00088 for(int i = 1; i<g+2; i++)
00089 {
00090 x += lanczos_coef[i]/(val+(double)i);
00091 }
00092
00093 double t (val + static_cast<double>(g) + 0.5);
00094
00095 return ( 2.5066282746310002 * pow( t, (val+0.5) ) * exp(-t) * x );
00096
00097 }
00098
00099 }
00100
00101
00102
00103
00104
00105
00106
00107
00108 double lngamma(double val)
00109 {
00110
00111 double inf( 9.0e+99 );
00112
00113 if( val <= 0.0 )
00114 {
00115 return inf;
00116 }
00117
00118
00119 int g(7);
00120 const double lanczos_coef[] = { 0.99999999999980993,
00121 676.5203681218851,
00122 -1259.1392167224028,
00123 771.32342877765313,
00124 -176.61502916214059,
00125 12.507343278686905,
00126 -0.13857109526572012,
00127 9.9843695780195716e-6,
00128 1.5056327351493116e-7 };
00129
00130 if(val < 0.5)
00131 {
00132 return ( 1.1447298858494002 - (log(sin(PI*val)) + lngamma(1.0-val)) );
00133 }
00134 else
00135 {
00136 val -= 1.0;
00137
00138 double x( lanczos_coef[0] );
00139
00140 for(int i = 1; i<g+2; i++)
00141 {
00142 x += lanczos_coef[i]/(val+(double)i);
00143 }
00144
00145 double t (val + static_cast<double>(g) + 0.5);
00146
00147 return ( 0.918938533204672741781 + (val+0.5)*log(t) + (-t) + log(x) );
00148 }
00149
00150
00151 }
00152
00153
00154
00155
00156 double kummerFunc(const double a, const double z);
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175 double kummerFunc(const double aval, const double zval)
00176 {
00177
00178 double eps(1.0e-15);
00179
00180 double z( abs(zval) );
00181 double a( abs(aval) );
00182
00183 double den(a);
00184
00185 double s(1.0);
00186
00187 double coef(1.0);
00188
00189 while( abs(coef) > eps)
00190 {
00191 coef = coef*z;
00192 den += 1.0;
00193 coef = coef/den;
00194 s += coef;
00195 }
00196
00197 return s;
00198
00199 }
00200
00201
00202
00203
00204 double lower_gamma(const double a, const double z)
00205 {
00206
00207 double zp( abs(z) );
00208 double ap( abs(a) );
00209
00210 double s( kummerFunc(ap, zp) );
00211
00212 return exp(log(zp)*ap) * exp(-zp) * s / ap;
00213
00214 }
00215
00216
00217
00218
00219 double upper_gamma(const double a, const double z)
00220 {
00221
00222 return ( gamma(a) - lower_gamma(a, z) );
00223
00224 }
00225
00226
00227
00228
00229 double gammaP(const double a, const double z)
00230 {
00231 return ( lower_gamma(a,z) / gamma(a) );
00232 }
00233
00234
00235
00236
00237 double gammaQ(const double a, const double z)
00238 {
00239 return ( 1.0 - gammaP(a,z) );
00240 }
00241
00242
00243
00244
00245
00246
00247
00248
00249 double factorial(const int n)
00250 {
00251
00252
00253 if( n < 0 ) return 0.0;
00254
00255 double facttable[] = { 1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0,
00256 40320.0, 362880.0, 3628800.0, 39916800.0,
00257 479001600.0, 6227020800.0, 87178291200.0,
00258 1307674368000.0 };
00259
00260
00261 if( (n >= 0) && (n <= 15) )
00262 {
00263 return facttable[n];
00264 }
00265
00266
00267 return gamma(n+1.0);
00268
00269 }
00270
00271
00272
00273 double factorial(const double d)
00274 {
00275
00276 double n(d);
00277
00278 double returnValue(0.0);
00279
00280 if (n < 0.0)
00281 {
00282 return returnValue;
00283
00284
00285
00286 }
00287 else if ((n == 0.0)||(n == 1.0))
00288 {
00289 returnValue = 1.0;
00290
00291 } else
00292 {
00293 returnValue = n * factorial(n - 1.0);
00294 }
00295
00296 return returnValue;
00297 }
00298
00299
00300
00301
00302 double erf1(const double x);
00303
00304
00305
00306 double erf2(const double x);
00307
00308
00309
00310 double erf3(const double x);
00311
00312
00313
00314 double erf4(const double x);
00315
00316
00317
00318 double erf5(const double x);
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 double erf(const double x)
00340 {
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438 double inf( 9e+99 );
00439
00440 if (x >= inf)
00441 {
00442 return 1.0;
00443 }
00444
00445 if (x <= -inf)
00446 {
00447 return -1.0;
00448 }
00449
00450 if ( fabs(x) < 0.84375 )
00451 {
00452 return erf1(x);
00453 }
00454 else if ( ( 0.84375 <= fabs(x) ) &&
00455 ( fabs(x) < 1.25 ) )
00456 {
00457 return erf2(x);
00458 }
00459 else if ( ( 1.25 <= fabs(x) ) &&
00460 ( fabs(x) < 2.857142) )
00461 {
00462 return erf3(x);
00463 }
00464 else if ( ( 2.857142 <= fabs(x) ) &&
00465 ( fabs(x) < 6.0 ) )
00466 {
00467 return erf4(x);
00468 }
00469
00470
00471 return erf5(x);
00472
00473 }
00474
00475
00476
00477
00478 double erf1(const double x)
00479 {
00480
00481 int i;
00482
00483
00484 frexp (x , &i);
00485
00486 if ( fabs( static_cast<double>(i) ) > 28.0 )
00487 {
00488 if( fabs( static_cast<double>(i) ) > 57.0 )
00489 {
00490 return ( 1.128379167095512586316e+0 * x );
00491 }
00492
00493 return ( x * ( 2.28379167095512586316e-01 ) );
00494 }
00495
00496 double z( x*x );
00497 double r( 1.28379167095512558561e-01
00498 + z * ( -3.25042107247001499370e-01
00499 + z * ( -2.84817495755985104766e-02
00500 + z * ( -5.77027029648944159157e-03
00501 + z * -2.37630166566501626084e-05 ) ) ) );
00502
00503
00504 double s( 1.0 + z * ( 3.97917223959155352819e-01
00505 + z * ( 6.50222499887672944485e-02
00506 + z * ( 5.08130628187576562776e-03
00507 + z * ( 1.32494738004321644526e-04
00508 + z * -3.96022827877536812320e-06 ) ) ) ) );
00509
00510 return ( x * (1.0 + r/s ) );
00511
00512 }
00513
00514
00515
00516
00517 double erf2(const double x)
00518 {
00519
00520 double s( fabs(x) - 1.0 );
00521
00522 double P( -2.36211856075265944077e-03
00523 + s * ( 4.14856118683748331666e-01
00524 + s * ( -3.72207876035701323847e-01
00525 + s * ( 3.18346619901161753674e-01
00526 + s * ( -1.10894694282396677476e-01
00527 + s * ( 3.54783043256182359371e-02
00528 + s * -2.16637559486879084300e-03 ) ) ) ) ) );
00529
00530 double Q( 1.0 + s * ( 1.06420880400844228286e-01
00531 + s * ( 5.40397917702171048937e-01
00532 + s * ( 7.18286544141962662868e-02
00533 + s * ( 1.26171219808761642112e-01
00534 + s * ( 1.36370839120290507362e-02
00535 + s * 1.19844998467991074170e-02 ) ) ) ) ) );
00536
00537 if( x >= 0.0 )
00538 {
00539 return ( 8.45062911510467529297e-01 + P/Q );
00540 }
00541 else
00542 {
00543 return ( -8.45062911510467529297e-01 - P/Q );
00544 }
00545
00546 }
00547
00548
00549
00550
00551 double erf3(const double xval)
00552 {
00553
00554 double x0(xval);
00555 double x( fabs(xval) );
00556 double s( 1.0/(x*x) );
00557
00558 double R( -9.86494403484714822705e-03
00559 + s * ( -6.93858572707181764372e-01
00560 + s * ( -1.05586262253232909814e+01
00561 + s * ( -6.23753324503260060396e+01
00562 + s * ( -1.62396669462573470355e+02
00563 + s * ( -1.84605092906711035994e+02
00564 + s * ( -8.12874355063065934246e+01
00565 + s * -9.81432934416914548592e+00 ) ) ) ) ) ) );
00566
00567 double S( 1.0 + s * ( 1.96512716674392571292e+01
00568 + s * ( 1.37657754143519042600e+02
00569 + s * ( 4.34565877475229228821e+02
00570 + s * ( 6.45387271733267880336e+02
00571 + s * ( 4.29008140027567833386e+02
00572 + s * ( 1.08635005541779435134e+02
00573 + s * ( 6.57024977031928170135e+00
00574 + s * -6.04244152148580987438e-02 ) ) ) ) ) ) ) );
00575
00576 double r( exp(-x0*x0-0.5625) * exp( (x0-x)*(x0+x)+R/S) );
00577
00578 if( x0 >= 0.0 )
00579 {
00580 return ( 1.0 - r/x );
00581 }
00582 else
00583 {
00584 return ( r/x - 1.0 );
00585 }
00586
00587 }
00588
00589
00590
00591
00592 double erf4(const double xval)
00593 {
00594
00595 double x0(xval);
00596 double x( fabs(xval) );
00597 double s( 1.0/(x*x) );
00598
00599 double R( -9.86494292470009928597e-03
00600 + s * ( -7.99283237680523006574e-01
00601 + s * ( -1.77579549177547519889e+01
00602 + s * ( -1.60636384855821916062e+02
00603 + s * ( -6.37566443368389627722e+02
00604 + s * ( -1.02509513161107724954e+03
00605 + s * -4.83519191608651397019e+02 ) ) ) ) ) );
00606
00607 double S( 1.0 + s * ( 3.03380607434824582924e+01
00608 + s * ( 3.25792512996573918826e+02
00609 + s * ( 1.53672958608443695994e+03
00610 + s * ( 3.19985821950859553908e+03
00611 + s * ( 2.55305040643316442583e+03
00612 + s * ( 4.74528541206955367215e+02
00613 + s * -2.24409524465858183362e+01 ) ) ) ) ) ) );
00614
00615 double r( exp( -x0 * x0 - 0.5625 ) * exp( (x0-x)*(x0+x) + R/S ) );
00616
00617 if( x0 >= 0.0 )
00618 {
00619 return ( 1.0 - r/x );
00620 }
00621 else
00622 {
00623 return ( r/x - 1.0 );
00624 }
00625
00626 }
00627
00628
00629
00630
00631 double erf5(const double x)
00632 {
00633
00634 double tiny( 1e-99 );
00635
00636 if ( x > 0.0 )
00637 {
00638 return ( 1.0 - tiny );
00639 }
00640 else
00641 {
00642 return ( tiny - 1.0 );
00643 }
00644
00645 }
00646
00647
00648
00649
00650 double erfc(const double x)
00651 {
00652
00653 return ( 1.0 - erf(x) );
00654
00655 }
00656
00657
00658
00659
00660
00661
00662
00663 double inverf(const double z)
00664 {
00665
00666 double inf( 9.0e+99 );
00667
00668
00669 if( z >= 1.0 ) return inf;
00670 if( z <= -1.0 ) return -inf;
00671
00672 double z2PI( z*z*PI );
00673
00674
00675
00676
00677 double x( 0.88622692545275794 * z
00678 * ( 1.0 + z2PI * ( 0.083333333333333333
00679 + z2PI * ( 0.014583333333333334
00680 + z2PI * ( 0.0031498015873015874
00681 + z2PI * ( 0.00075248704805996468
00682 + z2PI * ( 0.0001907475361251403
00683 + z2PI * ( 1.8780048076923078e-5
00684 + z2PI * ( 1.3623642420578133e-5 ) ) ) ) ) ) ) ) );
00685
00686 double delta(1.0);
00687 double threshold( 1.0e-10 );
00688 int iter( 0 );
00689
00690 while ( ( fabs(delta) > threshold ) &&
00691 ( iter < 100 ) )
00692 {
00693
00694
00695 delta = (z - gpstk::erf(x)) / (1.1283791670955126*std::exp(-x*x));
00696
00697 x = x + delta;
00698
00699 ++iter;
00700
00701 }
00702
00703 return x;
00704
00705 }
00706
00707
00708
00709
00710
00711
00712
00713 double beta(const double x, const double y)
00714 {
00715
00716 return ( gamma(x) * gamma(y) / gamma( x + y ) );
00717
00718 }
00719
00720
00721
00722
00723
00724
00725
00726 double lnbeta(double x, double y)
00727 {
00728
00729 x = fabs(x);
00730 y = fabs(y);
00731
00732 return ( lngamma(x) + lngamma(y) - lngamma( x + y ) );
00733
00734 }
00735
00736
00737
00738
00739
00740
00741
00742 double incompletebetaps(const double x, const double a, const double b)
00743 {
00744
00745
00746 double epsilon(1.0e-30);
00747 double big(1.0e99);
00748 double small(1.0e-99);
00749 double maxgam(171.0);
00750 double ai( 1.0/a );
00751 double u( (1.0-b) * x );
00752 double v( u / (a+1.0) );
00753 double t1(v);
00754 double t(u);
00755 double n(2.0);
00756 double s(0.0);
00757 double z( epsilon*ai );
00758
00759 while( std::fabs(v) > z )
00760 {
00761 u = (n-b) * x / n;
00762 t = t * u;
00763 v = t/(a+n);
00764 s = s+v;
00765 n = n+1.0;
00766 }
00767
00768 s = s + t1 + ai;
00769 u = a * std::log(x);
00770
00771
00772 if( (a+b < maxgam) &&
00773 ( std::fabs(u) < std::log(big) ) )
00774 {
00775
00776 t = gamma(a+b) / ( gamma(a) * gamma(b) );
00777 s = s * t * pow(x, a);
00778
00779 }
00780 else
00781 {
00782
00783 t = lngamma(a+b) - lngamma(a)- lngamma(b) + u + std::log(s);
00784
00785 if( t < std::log(small) )
00786 {
00787 s = 0.0;
00788 }
00789 else
00790 {
00791 s = std::exp(t);
00792 }
00793
00794 }
00795
00796 return s;
00797
00798 }
00799
00800
00801
00802
00803 double incompletebetafe(const double x, const double a, const double b)
00804 {
00805
00806
00807 double epsilon(1.0e-30);
00808 double big( 1.0e16 );
00809 double biginv( 1.0/big );
00810
00811 double k1( a );
00812 double k2( a+b );
00813 double k3( a );
00814 double k4( a+1.0 );
00815 double k5( 1.0 );
00816 double k6( b-1.0 );
00817 double k7( k4 );
00818 double k8( a+2.0 );
00819
00820 double pkm1( 1.0 );
00821 double pkm2( 0.0 );
00822 double qkm1( 1.0 );
00823 double qkm2( 1.0 );
00824
00825 double ans( 1.0 );
00826 double r( 1.0 );
00827 double thresh( 3.0 * epsilon );
00828
00829 int n( 0 );
00830
00831
00832 do
00833 {
00834
00835 double xk( -x * k1 * k2 / (k3*k4) );
00836 double pk( pkm1 + pkm2 * xk );
00837 double qk( qkm1 + qkm2 * xk );
00838
00839 pkm2 = pkm1;
00840 pkm1 = pk;
00841 qkm2 = qkm1;
00842 qkm1 = qk;
00843
00844 xk = x * k5 * k6 / (k7*k8);
00845 pk = pkm1 + pkm2 * xk;
00846 qk = qkm1 + qkm2 * xk;
00847
00848 pkm2 = pkm1;
00849 pkm1 = pk;
00850 qkm2 = qkm1;
00851 qkm1 = qk;
00852
00853 if( qk != 0.0 )
00854 {
00855 r = pk/qk;
00856 }
00857
00858 double t;
00859
00860 if( r != 0.0 )
00861 {
00862 t = std::fabs( (ans-r)/r );
00863 ans = r;
00864 }
00865 else
00866 {
00867 t = 1.0;
00868 }
00869
00870
00871 if( t < thresh )
00872 {
00873 break;
00874 }
00875
00876
00877 k1 = k1 + 1.0;
00878 k2 = k2 + 1.0;
00879 k3 = k3 + 2.0;
00880 k4 = k4 + 2.0;
00881 k5 = k5 + 1.0;
00882 k6 = k6 - 1.0;
00883 k7 = k7 + 2.0;
00884 k8 = k8 + 2.0;
00885
00886
00887 if( (std::fabs(qk) + std::fabs(pk)) > big )
00888 {
00889 pkm2 = pkm2 * biginv;
00890 pkm1 = pkm1 * biginv;
00891 qkm2 = qkm2 * biginv;
00892 qkm1 = qkm1 * biginv;
00893 }
00894
00895
00896 if( (std::fabs(qk) < biginv) ||
00897 (std::fabs(pk) < biginv ) )
00898 {
00899 pkm2 = pkm2 * big;
00900 pkm1 = pkm1 * big;
00901 qkm2 = qkm2 * big;
00902 qkm1 = qkm1 * big;
00903 }
00904
00905
00906 n = n+1;
00907
00908 }
00909 while( n < 300 );
00910
00911
00912 return ans;
00913
00914 }
00915
00916
00917
00918
00919 double incompletebetafe2(const double x, const double a, const double b)
00920 {
00921
00922
00923 double epsilon(1.0e-30);
00924 double big( 1.0e16 );
00925 double biginv( 1.0/big );
00926
00927 double k1( a );
00928 double k2( b-1.0 );
00929 double k3( a );
00930 double k4( a+1.0 );
00931 double k5( 1.0 );
00932 double k6( a+b );
00933 double k7( a+1.0 );
00934 double k8( a+2.0 );
00935
00936 double pkm1( 1.0 );
00937 double pkm2( 0.0 );
00938 double qkm1( 1.0 );
00939 double qkm2( 1.0 );
00940
00941 double z( x / (1.0-x) );
00942 double ans( 1.0 );
00943 double r( 1.0 );
00944 double thresh( 3.0 * epsilon );
00945
00946 int n( 0 );
00947
00948
00949 do
00950 {
00951
00952 double xk( -z * k1 * k2 / (k3*k4) );
00953 double pk( pkm1 + pkm2 * xk );
00954 double qk( qkm1 + qkm2 * xk );
00955
00956 pkm2 = pkm1;
00957 pkm1 = pk;
00958 qkm2 = qkm1;
00959 qkm1 = qk;
00960
00961 xk = z * k5 * k6 / (k7*k8);
00962 pk = pkm1 + pkm2 * xk;
00963 qk = qkm1 + qkm2 * xk;
00964
00965 pkm2 = pkm1;
00966 pkm1 = pk;
00967 qkm2 = qkm1;
00968 qkm1 = qk;
00969
00970 if( qk != 0.0 )
00971 {
00972 r = pk/qk;
00973 }
00974
00975 double t;
00976
00977 if( r != 0.0 )
00978 {
00979 t = std::fabs( (ans-r) / r );
00980 ans = r;
00981 }
00982 else
00983 {
00984 t = 1.0;
00985 }
00986
00987
00988 if( t < thresh )
00989 {
00990 break;
00991 }
00992
00993
00994
00995 k1 = k1 + 1.0;
00996 k2 = k2 - 1.0;
00997 k3 = k3 + 2.0;
00998 k4 = k4 + 2.0;
00999 k5 = k5 + 1.0;
01000 k6 = k6 + 1.0;
01001 k7 = k7 + 2.0;
01002 k8 = k8 + 2.0;
01003
01004
01005 if( (std::fabs(qk) + std::fabs(pk)) > big )
01006 {
01007 pkm2 = pkm2 * biginv;
01008 pkm1 = pkm1 * biginv;
01009 qkm2 = qkm2 * biginv;
01010 qkm1 = qkm1 * biginv;
01011 }
01012
01013
01014 if( (std::fabs(qk) < biginv) ||
01015 (std::fabs(pk) < biginv ) )
01016 {
01017 pkm2 = pkm2 * big;
01018 pkm1 = pkm1 * big;
01019 qkm2 = qkm2 * big;
01020 qkm1 = qkm1 * big;
01021 }
01022
01023
01024 n = n+1;
01025
01026 }
01027 while( n < 300 );
01028
01029
01030 return ans;
01031
01032 }
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 double regIncompleteBeta(const double x, const double a, const double b)
01043 throw(InvalidParameter)
01044 {
01045
01046
01047 if( (a <= 0.0) || (b <= 0.0) )
01048 {
01049 InvalidParameter e("Function 'regIncompleteBeta()': 'a' and 'b' must \
01050 be greater than zero." );
01051 GPSTK_THROW(e);
01052 }
01053
01054 if( (x < 0.0) || (x > 1.0) )
01055 {
01056 InvalidParameter e("Function 'regIncompleteBeta()': 'x' must be \
01057 within the interval [0,1]." );
01058 GPSTK_THROW(e);
01059 }
01060
01061 if (x == 0.0) return 0.0;
01062 if (x == 1.0) return 1.0;
01063
01064
01065 int flag(0);
01066 double epsilon(1.0e-30);
01067 double xx(x);
01068 double aa(a);
01069 double bb(b);
01070
01071
01072
01073 if( (bb*xx <= 1.0) &&
01074 (xx <= 0.95) )
01075 {
01076 return incompletebetaps(xx, aa, bb);
01077 }
01078
01079 double w( 1.0-xx );
01080 double t, xc;
01081
01082 if( xx > (aa/(aa+bb)) )
01083 {
01084 flag = 1;
01085 t = aa;
01086 aa = bb;
01087 bb = t;
01088 xc = xx;
01089 xx = w;
01090 }
01091 else
01092 {
01093 xc = w;
01094 }
01095
01096 double result(0.0);
01097
01098 if( (flag==1) &&
01099 (bb*xx <= 1.0) &&
01100 (xx <= 0.95) )
01101 {
01102
01103
01104 t = incompletebetaps(xx, aa, bb);
01105
01106 if( t <= epsilon )
01107 {
01108 result = 1.0 - epsilon;
01109 }
01110 else
01111 {
01112 result = 1.0 - t;
01113 }
01114
01115 return result;
01116
01117 }
01118
01119 double y( xx * (aa+bb-2.0) - (aa-1.0) );
01120
01121 if( y < 0.0 )
01122 {
01123 w = incompletebetafe(xx, aa, bb);
01124 }
01125 else
01126 {
01127 w = incompletebetafe2(xx, aa, bb) / xc;
01128 }
01129
01130 t = std::pow(xc, bb);
01131 t = t * std::pow(xx, aa);
01132 t = t/aa;
01133 t = t*w;
01134 t = t * ( gamma(aa+bb) / ( gamma(aa) * gamma(bb) ) );
01135
01136 if( flag==1 )
01137 {
01138 if( t <= epsilon )
01139 {
01140 result = 1.0 - epsilon;
01141 }
01142 else
01143 {
01144 result = 1.0 - t;
01145 }
01146 }
01147 else
01148 {
01149 result = t;
01150 }
01151
01152 return result;
01153
01154 }
01155
01156
01157
01158 }