00001 #pragma ident "$Id: SpecialFunctions.cpp 1408 2008-09-24 15:18:16Z architest $"
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
00274 double erf1(const double x);
00275
00276
00277
00278 double erf2(const double x);
00279
00280
00281
00282 double erf3(const double x);
00283
00284
00285
00286 double erf4(const double x);
00287
00288
00289
00290 double erf5(const double x);
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311 double erf(const double x)
00312 {
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
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 double inf( 9e+99 );
00411
00412 if (x >= inf)
00413 {
00414 return 1.0;
00415 }
00416
00417 if (x <= -inf)
00418 {
00419 return -1.0;
00420 }
00421
00422 if ( fabs(x) < 0.84375 )
00423 {
00424 return erf1(x);
00425 }
00426 else if ( ( 0.84375 <= fabs(x) ) &&
00427 ( fabs(x) < 1.25 ) )
00428 {
00429 return erf2(x);
00430 }
00431 else if ( ( 1.25 <= fabs(x) ) &&
00432 ( fabs(x) < 2.857142) )
00433 {
00434 return erf3(x);
00435 }
00436 else if ( ( 2.857142 <= fabs(x) ) &&
00437 ( fabs(x) < 6.0 ) )
00438 {
00439 return erf4(x);
00440 }
00441
00442
00443 return erf5(x);
00444
00445 }
00446
00447
00448
00449
00450 double erf1(const double x)
00451 {
00452
00453 int i;
00454
00455
00456 frexp (x , &i);
00457
00458 if ( fabs( static_cast<double>(i) ) > 28.0 )
00459 {
00460 if( fabs( static_cast<double>(i) ) > 57.0 )
00461 {
00462 return ( 1.128379167095512586316e+0 * x );
00463 }
00464
00465 return ( x * ( 2.28379167095512586316e-01 ) );
00466 }
00467
00468 double z( x*x );
00469 double r( 1.28379167095512558561e-01
00470 + z * ( -3.25042107247001499370e-01
00471 + z * ( -2.84817495755985104766e-02
00472 + z * ( -5.77027029648944159157e-03
00473 + z * -2.37630166566501626084e-05 ) ) ) );
00474
00475
00476 double s( 1.0 + z * ( 3.97917223959155352819e-01
00477 + z * ( 6.50222499887672944485e-02
00478 + z * ( 5.08130628187576562776e-03
00479 + z * ( 1.32494738004321644526e-04
00480 + z * -3.96022827877536812320e-06 ) ) ) ) );
00481
00482 return ( x * (1.0 + r/s ) );
00483
00484 }
00485
00486
00487
00488
00489 double erf2(const double x)
00490 {
00491
00492 double s( fabs(x) - 1.0 );
00493
00494 double P( -2.36211856075265944077e-03
00495 + s * ( 4.14856118683748331666e-01
00496 + s * ( -3.72207876035701323847e-01
00497 + s * ( 3.18346619901161753674e-01
00498 + s * ( -1.10894694282396677476e-01
00499 + s * ( 3.54783043256182359371e-02
00500 + s * -2.16637559486879084300e-03 ) ) ) ) ) );
00501
00502 double Q( 1.0 + s * ( 1.06420880400844228286e-01
00503 + s * ( 5.40397917702171048937e-01
00504 + s * ( 7.18286544141962662868e-02
00505 + s * ( 1.26171219808761642112e-01
00506 + s * ( 1.36370839120290507362e-02
00507 + s * 1.19844998467991074170e-02 ) ) ) ) ) );
00508
00509 if( x >= 0.0 )
00510 {
00511 return ( 8.45062911510467529297e-01 + P/Q );
00512 }
00513 else
00514 {
00515 return ( -8.45062911510467529297e-01 - P/Q );
00516 }
00517
00518 }
00519
00520
00521
00522
00523 double erf3(const double xval)
00524 {
00525
00526 double x0(xval);
00527 double x( fabs(xval) );
00528 double s( 1.0/(x*x) );
00529
00530 double R( -9.86494403484714822705e-03
00531 + s * ( -6.93858572707181764372e-01
00532 + s * ( -1.05586262253232909814e+01
00533 + s * ( -6.23753324503260060396e+01
00534 + s * ( -1.62396669462573470355e+02
00535 + s * ( -1.84605092906711035994e+02
00536 + s * ( -8.12874355063065934246e+01
00537 + s * -9.81432934416914548592e+00 ) ) ) ) ) ) );
00538
00539 double S( 1.0 + s * ( 1.96512716674392571292e+01
00540 + s * ( 1.37657754143519042600e+02
00541 + s * ( 4.34565877475229228821e+02
00542 + s * ( 6.45387271733267880336e+02
00543 + s * ( 4.29008140027567833386e+02
00544 + s * ( 1.08635005541779435134e+02
00545 + s * ( 6.57024977031928170135e+00
00546 + s * -6.04244152148580987438e-02 ) ) ) ) ) ) ) );
00547
00548 double r( exp(-x0*x0-0.5625) * exp( (x0-x)*(x0+x)+R/S) );
00549
00550 if( x0 >= 0.0 )
00551 {
00552 return ( 1.0 - r/x );
00553 }
00554 else
00555 {
00556 return ( r/x - 1.0 );
00557 }
00558
00559 }
00560
00561
00562
00563
00564 double erf4(const double xval)
00565 {
00566
00567 double x0(xval);
00568 double x( fabs(xval) );
00569 double s( 1.0/(x*x) );
00570
00571 double R( -9.86494292470009928597e-03
00572 + s * ( -7.99283237680523006574e-01
00573 + s * ( -1.77579549177547519889e+01
00574 + s * ( -1.60636384855821916062e+02
00575 + s * ( -6.37566443368389627722e+02
00576 + s * ( -1.02509513161107724954e+03
00577 + s * -4.83519191608651397019e+02 ) ) ) ) ) );
00578
00579 double S( 1.0 + s * ( 3.03380607434824582924e+01
00580 + s * ( 3.25792512996573918826e+02
00581 + s * ( 1.53672958608443695994e+03
00582 + s * ( 3.19985821950859553908e+03
00583 + s * ( 2.55305040643316442583e+03
00584 + s * ( 4.74528541206955367215e+02
00585 + s * -2.24409524465858183362e+01 ) ) ) ) ) ) );
00586
00587 double r( exp( -x0 * x0 - 0.5625 ) * exp( (x0-x)*(x0+x) + R/S ) );
00588
00589 if( x0 >= 0.0 )
00590 {
00591 return ( 1.0 - r/x );
00592 }
00593 else
00594 {
00595 return ( r/x - 1.0 );
00596 }
00597
00598 }
00599
00600
00601
00602
00603 double erf5(const double x)
00604 {
00605
00606 double tiny( 1e-99 );
00607
00608 if ( x > 0.0 )
00609 {
00610 return ( 1.0 - tiny );
00611 }
00612 else
00613 {
00614 return ( tiny - 1.0 );
00615 }
00616
00617 }
00618
00619
00620
00621
00622 double erfc(const double x)
00623 {
00624
00625 return ( 1.0 - erf(x) );
00626
00627 }
00628
00629
00630
00631
00632
00633
00634
00635 double inverf(const double z)
00636 {
00637
00638 double inf( 9.0e+99 );
00639
00640
00641 if( z >= 1.0 ) return inf;
00642 if( z <= -1.0 ) return -inf;
00643
00644 double z2PI( z*z*PI );
00645
00646
00647
00648
00649 double x( 0.88622692545275794 * z
00650 * ( 1.0 + z2PI * ( 0.083333333333333333
00651 + z2PI * ( 0.014583333333333334
00652 + z2PI * ( 0.0031498015873015874
00653 + z2PI * ( 0.00075248704805996468
00654 + z2PI * ( 0.0001907475361251403
00655 + z2PI * ( 1.8780048076923078e-5
00656 + z2PI * ( 1.3623642420578133e-5 ) ) ) ) ) ) ) ) );
00657
00658 double delta(1.0);
00659 double threshold( 1.0e-10 );
00660 int iter( 0 );
00661
00662 while ( ( fabs(delta) > threshold ) &&
00663 ( iter < 100 ) )
00664 {
00665
00666
00667 delta = (z - gpstk::erf(x)) / (1.1283791670955126*std::exp(-x*x));
00668
00669 x = x + delta;
00670
00671 ++iter;
00672
00673 }
00674
00675 return x;
00676
00677 }
00678
00679
00680
00681
00682
00683
00684
00685 double beta(const double x, const double y)
00686 {
00687
00688 return ( gamma(x) * gamma(y) / gamma( x + y ) );
00689
00690 }
00691
00692
00693
00694
00695
00696
00697
00698 double lnbeta(double x, double y)
00699 {
00700
00701 x = fabs(x);
00702 y = fabs(y);
00703
00704 return ( lngamma(x) + lngamma(y) - lngamma( x + y ) );
00705
00706 }
00707
00708
00709
00710
00711
00712
00713
00714 double incompletebetaps(const double x, const double a, const double b)
00715 {
00716
00717
00718 double epsilon(1.0e-30);
00719 double big(1.0e99);
00720 double small(1.0e-99);
00721 double maxgam(171.0);
00722 double ai( 1.0/a );
00723 double u( (1.0-b) * x );
00724 double v( u / (a+1.0) );
00725 double t1(v);
00726 double t(u);
00727 double n(2.0);
00728 double s(0.0);
00729 double z( epsilon*ai );
00730
00731 while( std::fabs(v) > z )
00732 {
00733 u = (n-b) * x / n;
00734 t = t * u;
00735 v = t/(a+n);
00736 s = s+v;
00737 n = n+1.0;
00738 }
00739
00740 s = s + t1 + ai;
00741 u = a * std::log(x);
00742
00743
00744 if( (a+b < maxgam) &&
00745 ( std::fabs(u) < std::log(big) ) )
00746 {
00747
00748 t = gamma(a+b) / ( gamma(a) * gamma(b) );
00749 s = s * t * pow(x, a);
00750
00751 }
00752 else
00753 {
00754
00755 t = lngamma(a+b) - lngamma(a)- lngamma(b) + u + std::log(s);
00756
00757 if( t < std::log(small) )
00758 {
00759 s = 0.0;
00760 }
00761 else
00762 {
00763 s = std::exp(t);
00764 }
00765
00766 }
00767
00768 return s;
00769
00770 }
00771
00772
00773
00774
00775 double incompletebetafe(const double x, const double a, const double b)
00776 {
00777
00778
00779 double epsilon(1.0e-30);
00780 double big( 1.0e16 );
00781 double biginv( 1.0/big );
00782
00783 double k1( a );
00784 double k2( a+b );
00785 double k3( a );
00786 double k4( a+1.0 );
00787 double k5( 1.0 );
00788 double k6( b-1.0 );
00789 double k7( k4 );
00790 double k8( a+2.0 );
00791
00792 double pkm1( 1.0 );
00793 double pkm2( 0.0 );
00794 double qkm1( 1.0 );
00795 double qkm2( 1.0 );
00796
00797 double ans( 1.0 );
00798 double r( 1.0 );
00799 double thresh( 3.0 * epsilon );
00800
00801 int n( 0 );
00802
00803
00804 do
00805 {
00806
00807 double xk( -x * k1 * k2 / (k3*k4) );
00808 double pk( pkm1 + pkm2 * xk );
00809 double qk( qkm1 + qkm2 * xk );
00810
00811 pkm2 = pkm1;
00812 pkm1 = pk;
00813 qkm2 = qkm1;
00814 qkm1 = qk;
00815
00816 xk = x * k5 * k6 / (k7*k8);
00817 pk = pkm1 + pkm2 * xk;
00818 qk = qkm1 + qkm2 * xk;
00819
00820 pkm2 = pkm1;
00821 pkm1 = pk;
00822 qkm2 = qkm1;
00823 qkm1 = qk;
00824
00825 if( qk != 0.0 )
00826 {
00827 r = pk/qk;
00828 }
00829
00830 double t;
00831
00832 if( r != 0.0 )
00833 {
00834 t = std::fabs( (ans-r)/r );
00835 ans = r;
00836 }
00837 else
00838 {
00839 t = 1.0;
00840 }
00841
00842
00843 if( t < thresh )
00844 {
00845 break;
00846 }
00847
00848
00849 k1 = k1 + 1.0;
00850 k2 = k2 + 1.0;
00851 k3 = k3 + 2.0;
00852 k4 = k4 + 2.0;
00853 k5 = k5 + 1.0;
00854 k6 = k6 - 1.0;
00855 k7 = k7 + 2.0;
00856 k8 = k8 + 2.0;
00857
00858
00859 if( (std::fabs(qk) + std::fabs(pk)) > big )
00860 {
00861 pkm2 = pkm2 * biginv;
00862 pkm1 = pkm1 * biginv;
00863 qkm2 = qkm2 * biginv;
00864 qkm1 = qkm1 * biginv;
00865 }
00866
00867
00868 if( (std::fabs(qk) < biginv) ||
00869 (std::fabs(pk) < biginv ) )
00870 {
00871 pkm2 = pkm2 * big;
00872 pkm1 = pkm1 * big;
00873 qkm2 = qkm2 * big;
00874 qkm1 = qkm1 * big;
00875 }
00876
00877
00878 n = n+1;
00879
00880 }
00881 while( n < 300 );
00882
00883
00884 return ans;
00885
00886 }
00887
00888
00889
00890
00891 double incompletebetafe2(const double x, const double a, const double b)
00892 {
00893
00894
00895 double epsilon(1.0e-30);
00896 double big( 1.0e16 );
00897 double biginv( 1.0/big );
00898
00899 double k1( a );
00900 double k2( b-1.0 );
00901 double k3( a );
00902 double k4( a+1.0 );
00903 double k5( 1.0 );
00904 double k6( a+b );
00905 double k7( a+1.0 );
00906 double k8( a+2.0 );
00907
00908 double pkm1( 1.0 );
00909 double pkm2( 0.0 );
00910 double qkm1( 1.0 );
00911 double qkm2( 1.0 );
00912
00913 double z( x / (1.0-x) );
00914 double ans( 1.0 );
00915 double r( 1.0 );
00916 double thresh( 3.0 * epsilon );
00917
00918 int n( 0 );
00919
00920
00921 do
00922 {
00923
00924 double xk( -z * k1 * k2 / (k3*k4) );
00925 double pk( pkm1 + pkm2 * xk );
00926 double qk( qkm1 + qkm2 * xk );
00927
00928 pkm2 = pkm1;
00929 pkm1 = pk;
00930 qkm2 = qkm1;
00931 qkm1 = qk;
00932
00933 xk = z * k5 * k6 / (k7*k8);
00934 pk = pkm1 + pkm2 * xk;
00935 qk = qkm1 + qkm2 * xk;
00936
00937 pkm2 = pkm1;
00938 pkm1 = pk;
00939 qkm2 = qkm1;
00940 qkm1 = qk;
00941
00942 if( qk != 0.0 )
00943 {
00944 r = pk/qk;
00945 }
00946
00947 double t;
00948
00949 if( r != 0.0 )
00950 {
00951 t = std::fabs( (ans-r) / r );
00952 ans = r;
00953 }
00954 else
00955 {
00956 t = 1.0;
00957 }
00958
00959
00960 if( t < thresh )
00961 {
00962 break;
00963 }
00964
00965
00966
00967 k1 = k1 + 1.0;
00968 k2 = k2 - 1.0;
00969 k3 = k3 + 2.0;
00970 k4 = k4 + 2.0;
00971 k5 = k5 + 1.0;
00972 k6 = k6 + 1.0;
00973 k7 = k7 + 2.0;
00974 k8 = k8 + 2.0;
00975
00976
00977 if( (std::fabs(qk) + std::fabs(pk)) > big )
00978 {
00979 pkm2 = pkm2 * biginv;
00980 pkm1 = pkm1 * biginv;
00981 qkm2 = qkm2 * biginv;
00982 qkm1 = qkm1 * biginv;
00983 }
00984
00985
00986 if( (std::fabs(qk) < biginv) ||
00987 (std::fabs(pk) < biginv ) )
00988 {
00989 pkm2 = pkm2 * big;
00990 pkm1 = pkm1 * big;
00991 qkm2 = qkm2 * big;
00992 qkm1 = qkm1 * big;
00993 }
00994
00995
00996 n = n+1;
00997
00998 }
00999 while( n < 300 );
01000
01001
01002 return ans;
01003
01004 }
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014 double regIncompleteBeta(const double x, const double a, const double b)
01015 throw(InvalidParameter)
01016 {
01017
01018
01019 if( (a <= 0.0) || (b <= 0.0) )
01020 {
01021 InvalidParameter e("Function 'regIncompleteBeta()': 'a' and 'b' must \
01022 be greater than zero." );
01023 GPSTK_THROW(e);
01024 }
01025
01026 if( (x < 0.0) || (x > 1.0) )
01027 {
01028 InvalidParameter e("Function 'regIncompleteBeta()': 'x' must be \
01029 within the interval [0,1]." );
01030 GPSTK_THROW(e);
01031 }
01032
01033 if (x == 0.0) return 0.0;
01034 if (x == 1.0) return 1.0;
01035
01036
01037 int flag(0);
01038 double epsilon(1.0e-30);
01039 double xx(x);
01040 double aa(a);
01041 double bb(b);
01042
01043
01044
01045 if( (bb*xx <= 1.0) &&
01046 (xx <= 0.95) )
01047 {
01048 return incompletebetaps(xx, aa, bb);
01049 }
01050
01051 double w( 1.0-xx );
01052 double t, xc;
01053
01054 if( xx > (aa/(aa+bb)) )
01055 {
01056 flag = 1;
01057 t = aa;
01058 aa = bb;
01059 bb = t;
01060 xc = xx;
01061 xx = w;
01062 }
01063 else
01064 {
01065 xc = w;
01066 }
01067
01068 double result(0.0);
01069
01070 if( (flag==1) &&
01071 (bb*xx <= 1.0) &&
01072 (xx <= 0.95) )
01073 {
01074
01075
01076 t = incompletebetaps(xx, aa, bb);
01077
01078 if( t <= epsilon )
01079 {
01080 result = 1.0 - epsilon;
01081 }
01082 else
01083 {
01084 result = 1.0 - t;
01085 }
01086
01087 return result;
01088
01089 }
01090
01091 double y( xx * (aa+bb-2.0) - (aa-1.0) );
01092
01093 if( y < 0.0 )
01094 {
01095 w = incompletebetafe(xx, aa, bb);
01096 }
01097 else
01098 {
01099 w = incompletebetafe2(xx, aa, bb) / xc;
01100 }
01101
01102 t = std::pow(xc, bb);
01103 t = t * std::pow(xx, aa);
01104 t = t/aa;
01105 t = t*w;
01106 t = t * ( gamma(aa+bb) / ( gamma(aa) * gamma(bb) ) );
01107
01108 if( flag==1 )
01109 {
01110 if( t <= epsilon )
01111 {
01112 result = 1.0 - epsilon;
01113 }
01114 else
01115 {
01116 result = 1.0 - t;
01117 }
01118 }
01119 else
01120 {
01121 result = t;
01122 }
01123
01124 return result;
01125
01126 }
01127
01128
01129
01130 }