ADMB Documentation  Fournier-pthread.1088
 All Classes Namespaces Files Functions Variables Typedefs Friends Defines
df13incbet.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df13incbet.cpp 494 2012-06-13 20:41:16Z johnoel $
00003  *
00004  * Based on code from Cephes C and C++ language special functions math library
00005  * and adapeted for automatic differentiation.
00006  * http://www.moshier.net/#Cephes or
00007  * http://www.netlib.org/cephes/
00008  *
00009  * Author: David Fournier
00010  * Copyright (c) 2009-2012 ADMB Foundation
00011  */
00016 #include <df13fun.h>
00017 
00018 static int mtherr(char *s, int n);
00019 static df1_three_variable igam(const df1_three_variable & a,
00020              const df1_three_variable & x);
00021 double pseries(const double & _a, const double & _b, const double & _x);
00022 double lgam(double x);
00023 
00024 static double MAXLOG = 200;
00025 static double MINLOG = -200;
00026 static double big = 4.503599627370496e15;
00027 static double biginv = 2.22044604925031308085e-16;
00028 static double MACHEP = 2.22045e-16;
00029 /*              gamma.c
00030  *
00031  *  Gamma function
00032  *
00033  *
00034  *
00035  * SYNOPSIS:
00036  *
00037  * double x, y, gamma();
00038  * extern int sgngam;
00039  *
00040  * y = gamma( x );
00041  *
00042  *
00043  *
00044  * DESCRIPTION:
00045  *
00046  * Returns gamma function of the argument.  The result is
00047  * correctly signed, and the sign (+1 or -1) is also
00048  * returned in a global (extern) variable named sgngam.
00049  * This variable is also filled in by the logarithmic gamma
00050  * function lgam().
00051  *
00052  * Arguments |x| <= 34 are reduced by recurrence and the function
00053  * approximated by a rational function of degree 6/7 in the
00054  * interval (2,3).  Large arguments are handled by Stirling's
00055  * formula. Large negative arguments are made positive using
00056  * a reflection formula.  
00057  *
00058  *
00059  * ACCURACY:
00060  *
00061  *                      Relative error:
00062  * arithmetic   domain     # trials      peak         rms
00063  *    DEC      -34, 34      10000       1.3e-16     2.5e-17
00064  *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
00065  *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
00066  *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
00067  *
00068  * Error for arguments outside the test range will be larger
00069  * owing to error amplification by the exponential function.
00070  *
00071  */
00072 /*              lgam()
00073  *
00074  *  Natural logarithm of gamma function
00075  *
00076  *
00077  *
00078  * SYNOPSIS:
00079  *
00080  * double x, y, lgam();
00081  * extern int sgngam;
00082  *
00083  * y = lgam( x );
00084  *
00085  *
00086  *
00087  * DESCRIPTION:
00088  *
00089  * Returns the base e (2.718...) logarithm of the absolute
00090  * value of the gamma function of the argument.
00091  * The sign (+1 or -1) of the gamma function is returned in a
00092  * global (extern) variable named sgngam.
00093  *
00094  * For arguments greater than 13, the logarithm of the gamma
00095  * function is approximated by the logarithmic version of
00096  * Stirling's formula using a polynomial approximation of
00097  * degree 4. Arguments between -33 and +33 are reduced by
00098  * recurrence to the interval [2,3] of a rational approximation.
00099  * The cosecant reflection formula is employed for arguments
00100  * less than -33.
00101  *
00102  * Arguments greater than MAXLGM return MAXNUM and an error
00103  * message.  MAXLGM = 2.035093e36 for DEC
00104  * arithmetic or 2.556348e305 for IEEE arithmetic.
00105  *
00106  *
00107  *
00108  * ACCURACY:
00109  *
00110  *
00111  * arithmetic      domain        # trials     peak         rms
00112  *    DEC     0, 3                  7000     5.2e-17     1.3e-17
00113  *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
00114  *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
00115  *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
00116  * The error criterion was relative when the function magnitude
00117  * was greater than one but absolute when it was less than one.
00118  *
00119  * The following test used the relative error criterion, though
00120  * at certain points the relative error could be much higher than
00121  * indicated.
00122  *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
00123  *
00124  */
00125 
00126 /*              gamma.c */
00127 /*  gamma function  */
00128 
00129 static double P[] = {
00130    1.60119522476751861407E-4,
00131    1.19135147006586384913E-3,
00132    1.04213797561761569935E-2,
00133    4.76367800457137231464E-2,
00134    2.07448227648435975150E-1,
00135    4.94214826801497100753E-1,
00136    9.99999999999999996796E-1
00137 };
00138 
00139 static double Q[] = {
00140    -2.31581873324120129819E-5,
00141    5.39605580493303397842E-4,
00142    -4.45641913851797240494E-3,
00143    1.18139785222060435552E-2,
00144    3.58236398605498653373E-2,
00145    -2.34591795718243348568E-1,
00146    7.14304917030273074085E-2,
00147    1.00000000000000000320E0
00148 };
00149 
00150 static double MAXGAM = 171.624376956302725;
00151 static double LOGPI = 1.14472988584940017414;
00152 
00153 /* Stirling's formula for the gamma function */
00154 static double STIR[5] = {
00155    7.87311395793093628397E-4,
00156    -2.29549961613378126380E-4,
00157    -2.68132617805781232825E-3,
00158    3.47222221605458667310E-3,
00159    8.33333333333482257126E-2,
00160 };
00161 static double MAXSTIR = 143.01608;
00162 static double SQTPI = 2.50662827463100050242E0;
00163 
00164 static int sgngam = 0;
00165 static double MAXNUM = 1.7976931348623158E+308;
00166 
00167 static double polevl(double, void *, int);
00168 static  df1_three_variable polevl(const df1_three_variable &, void *, int);
00169 static  double p1evl(double, void *, int);
00170 static  df1_three_variable p1evl(const df1_three_variable &, void *, int);
00171 //static double lgam(double);
00172 
00173 const double MYINF = 1.7976931348623158E+308;
00174 
00180 static df1_three_variable stirf(const df1_three_variable & _x)
00181 {
00182    ADUNCONST(df1_three_variable, x) df1_three_variable y, w, v;
00183 
00184    w = 1.0 / x;
00185    w = 1.0 + w * polevl(w, STIR, 4);
00186    y = exp(x);
00187    if (value(x) > MAXSTIR)
00188    {        /* Avoid overflow in pow() */
00189       v = pow(x, 0.5 * x - 0.25);
00190       y = v * (v / y);
00191    } else
00192    {
00193       y = pow(x, x - 0.5) / y;
00194    }
00195    y = SQTPI * y * w;
00196    return (y);
00197 }
00198 
00199 static double stirf(double _x)
00200 {
00201    double x = _x;
00202    double y, w, v;
00203 
00204    w = 1.0 / x;
00205    w = 1.0 + w * polevl(w, STIR, 4);
00206    y = exp(x);
00207    if (x > MAXSTIR)
00208    {        /* Avoid overflow in pow() */
00209       v = pow(x, 0.5 * x - 0.25);
00210       y = v * (v / y);
00211    } else
00212    {
00213       y = pow(x, x - 0.5) / y;
00214    }
00215    y = SQTPI * y * w;
00216    return (y);
00217 }
00218 
00219 
00224 static df1_three_variable gamma(const df1_three_variable & xx1)
00225 {
00226    df1_three_variable x;
00227    x = xx1;
00228    df1_three_variable MYBIG;
00229    MYBIG = 1.e+300;
00230    df1_three_variable p, q, z;
00231    df1_three_variable zero;
00232    zero = 0.0;
00233    int i;
00234 
00235    sgngam = 1;
00236 
00237 #ifdef NANS
00238    if (isnan(x))
00239       return (x);
00240 #endif
00241 #ifdef INFINITIES
00242 #ifdef NANS
00243    if (value(x) == MYINF)
00244       return (x);
00245    if (value(x) == -MYINF)
00246       return (NAN);
00247 #else
00248    if (!isfinite(value(x)))
00249       return (x);
00250 #endif
00251 #endif
00252    q = fabs(x);
00253 
00254    if (value(q) > 33.0)
00255    {
00256       if (value(x) < 0.0)
00257       {
00258    p = floor(value(q));
00259    if (value(p) == value(q))
00260    {
00261 #ifdef NANS
00262     gamnan:
00263       cerr << "gamma DOMAIN" << endl;
00264       return (zero);
00265 #else
00266       goto goverf;
00267 #endif
00268    }
00269    i = value(p);
00270    if ((i & 1) == 0)
00271       sgngam = -1;
00272    z = q - p;
00273    if (value(z) > 0.5)
00274    {
00275       p += 1.0;
00276       z = q - p;
00277    }
00278    z = q * sin(PI * z);
00279    if (value(z) == 0.0)
00280    {
00281 #ifdef INFINITIES
00282       return (sgngam * MYINF);
00283 #else
00284     goverf:
00285       mtherr("gamma", OVERFLOW);
00286       return (sgngam * MYBIG);
00287       //return( sgngam * MAXNUM);
00288 #endif
00289    }
00290    z = fabs(z);
00291    z = PI / (z * stirf(q));
00292       } else
00293       {
00294    z = stirf(x);
00295       }
00296       return (sgngam * z);
00297    }
00298 
00299    z = 1.0;
00300    while (value(x) >= 3.0)
00301    {
00302       x -= 1.0;
00303       z *= x;
00304    }
00305 
00306    while (value(x) < 0.0)
00307    {
00308       if (value(x) > -1.E-9)
00309    goto SMALL;
00310       z /= x;
00311       x += 1.0;
00312    }
00313 
00314    while (value(x) < 2.0)
00315    {
00316       if (value(x) < 1.e-9)
00317    goto SMALL;
00318       z /= x;
00319       x += 1.0;
00320    }
00321 
00322    if (value(x) == 2.0)
00323       return (z);
00324 
00325    x -= 2.0;
00326    p = polevl(x, P, 6);
00327    q = polevl(x, Q, 7);
00328    return (z * p / q);
00329 
00330  SMALL:
00331    if (value(x) == 0.0)
00332    {
00333 #ifdef INFINITIES
00334 #ifdef NANS
00335       goto gamnan;
00336 #else
00337       return (MYINF);
00338 #endif
00339 #else
00340       cerr << "gamma SING " << endl;
00341       return (MYBIG);
00342 #endif
00343    } else
00344       return (z / ((1.0 + 0.5772156649015329 * x) * x));
00345 }
00346 
00347 
00348 double gamma(double xx1)
00349 {
00350    double x = xx1;
00351    double MYBIG = 1.e+300;
00352    double p, q, z;
00353    double zero;
00354    zero = 0.0;
00355    int i;
00356 
00357    sgngam = 1;
00358 
00359 #ifdef NANS
00360    if (isnan(x))
00361       return (x);
00362 #endif
00363 #ifdef INFINITIES
00364 #ifdef NANS
00365    if (x == MYINF)
00366       return (x);
00367    if (x == -MYINF)
00368       return (NAN);
00369 #else
00370    if (!isfinite(x))
00371       return (x);
00372 #endif
00373 #endif
00374    q = fabs(x);
00375 
00376    if (q > 33.0)
00377    {
00378       if (x < 0.0)
00379       {
00380    p = floor(q);
00381    if (p == q)
00382    {
00383 #ifdef NANS
00384     gamnan:
00385       cerr << "gamma DOMAIN" << endl;
00386       return (zero);
00387 #else
00388       goto goverf;
00389 #endif
00390    }
00391    i = p;
00392    if ((i & 1) == 0)
00393       sgngam = -1;
00394    z = q - p;
00395    if (z > 0.5)
00396    {
00397       p += 1.0;
00398       z = q - p;
00399    }
00400    z = q * sin(PI * z);
00401    if (z == 0.0)
00402    {
00403 #ifdef INFINITIES
00404       return (sgngam * MYINF);
00405 #else
00406     goverf:
00407       mtherr("gamma", OVERFLOW);
00408       return (sgngam * MYBIG);
00409       //return( sgngam * MAXNUM);
00410 #endif
00411    }
00412    z = fabs(z);
00413    z = PI / (z * stirf(q));
00414       } else
00415       {
00416    z = stirf(x);
00417       }
00418       return (sgngam * z);
00419    }
00420 
00421    z = 1.0;
00422    while (x >= 3.0)
00423    {
00424       x -= 1.0;
00425       z *= x;
00426    }
00427 
00428    while (x < 0.0)
00429    {
00430       if (x > -1.E-9)
00431    goto SMALL;
00432       z /= x;
00433       x += 1.0;
00434    }
00435 
00436    while (x < 2.0)
00437    {
00438       if (x < 1.e-9)
00439    goto SMALL;
00440       z /= x;
00441       x += 1.0;
00442    }
00443 
00444    if (x == 2.0)
00445       return (z);
00446 
00447    x -= 2.0;
00448    p = polevl(x, P, 6);
00449    q = polevl(x, Q, 7);
00450    return (z * p / q);
00451 
00452  SMALL:
00453    if (x == 0.0)
00454    {
00455 #ifdef INFINITIES
00456 #ifdef NANS
00457       goto gamnan;
00458 #else
00459       return (MYINF);
00460 #endif
00461 #else
00462       cerr << "gamma SING " << endl;
00463       return (MYBIG);
00464 #endif
00465    } else
00466       return (z / ((1.0 + 0.5772156649015329 * x) * x));
00467 }
00468 
00469 
00470 
00471 
00472 /* A[]: Stirling's formula expansion of log gamma
00473  * B[], C[]: log gamma function between 2 and 3
00474  */
00475 static double A[] = {
00476    8.11614167470508450300E-4,
00477    -5.95061904284301438324E-4,
00478    7.93650340457716943945E-4,
00479    -2.77777777730099687205E-3,
00480    8.33333333333331927722E-2
00481 };
00482 
00483 static double B[] = {
00484    -1.37825152569120859100E3,
00485    -3.88016315134637840924E4,
00486    -3.31612992738871184744E5,
00487    -1.16237097492762307383E6,
00488    -1.72173700820839662146E6,
00489    -8.53555664245765465627E5
00490 };
00491 
00492 static double C[] = {
00493 /* 1.00000000000000000000E0, */
00494    -3.51815701436523470549E2,
00495    -1.70642106651881159223E4,
00496    -2.20528590553854454839E5,
00497    -1.13933444367982507207E6,
00498    -2.53252307177582951285E6,
00499    -2.01889141433532773231E6
00500 };
00501 
00502 /* log( sqrt( 2*pi ) ) */
00503 static double LS2PI = 0.91893853320467274178;
00504 static double MAXLGM = 2.556348e305;
00505 
00506 /* Logarithm of gamma function */
00507 
00508 int operator <(const df1_three_variable & x, double n)
00509 {
00510    return value(x) < n;
00511 }
00512 
00513 int operator >(const df1_three_variable & x, double n)
00514 {
00515    return value(x) > n;
00516 }
00517 
00518 int operator >=(const df1_three_variable & x, double n)
00519 {
00520    return value(x) >= n;
00521 }
00522 
00523 int operator ==(const df1_three_variable & x, const df1_three_variable & n)
00524 {
00525    return value(x) == value(n);
00526 }
00527 
00528 int operator ==(const df1_three_variable & x, double n)
00529 {
00530    return value(x) == n;
00531 }
00532 
00533 int operator ==(double x, const df1_three_variable & n)
00534 {
00535    return x == value(n);
00536 }
00537 
00538 int operator <(const df1_three_variable & x, const df1_three_variable & n)
00539 {
00540    return value(x) < value(n);
00541 }
00542 
00543 int operator >(const df1_three_variable & x, const df1_three_variable & n)
00544 {
00545    return value(x) > value(n);
00546 }
00547 
00552 df1_three_variable lgam(const df1_three_variable & _x)
00553 {
00554    df1_three_variable x, p, q, u, w, z, p1;
00555    x = _x;
00556    int i;
00557 
00558    sgngam = 1;
00559 #ifdef NANS
00560    if (isnan(x))
00561       return (x);
00562 #endif
00563 
00564 #ifdef INFINITIES
00565    if (!isfinite(x))
00566       return (MYINF);
00567 #endif
00568 
00569    if (x < -34.0)
00570    {
00571       q = -x;
00572       w = lgam(q);    /* note this modifies sgngam! */
00573       p = floor(value(q));
00574       if (p == q)
00575       {
00576        lgsing:
00577 #ifdef INFINITIES
00578    mtherr("lgam", SING);
00579    return (MYINF);
00580 #else
00581    goto loverf;
00582 #endif
00583       }
00584       i = value(p);
00585       if ((i & 1) == 0)
00586    sgngam = -1;
00587       else
00588    sgngam = 1;
00589       z = q - p;
00590       if (z > 0.5)
00591       {
00592    p += 1.0;
00593    z = p - q;
00594       }
00595       z = q * sin(PI * z);
00596       if (z == 0.0)
00597    goto lgsing;
00598 /*  z = log(PI) - log( z ) - w;*/
00599       z = LOGPI - log(z) - w;
00600       return (z);
00601    }
00602 
00603    if (x < 13.0)
00604    {
00605       z = 1.0;
00606       p = 0.0;
00607       u = x;
00608       while (u >= 3.0)
00609       {
00610    p -= 1.0;
00611    u = x + p;
00612    z *= u;
00613       }
00614       while (u < 2.0)
00615       {
00616    if (u == 0.0)
00617       goto lgsing;
00618    z /= u;
00619    p += 1.0;
00620    u = x + p;
00621       }
00622       if (z < 0.0)
00623       {
00624    sgngam = -1;
00625    z = -z;
00626       } else
00627    sgngam = 1;
00628       if (u == 2.0)
00629    return (log(z));
00630       p -= 2.0;
00631       x = x + p;
00632       p = x * polevl(x, B, 5) / p1evl(x, C, 6);
00633       return (log(z) + p);
00634    }
00635 
00636    if (x > MAXLGM)
00637    {
00638 #ifdef INFINITIES
00639       return (sgngam * MYINF);
00640 #else
00641     loverf:
00642       mtherr("lgam", OVERFLOW);
00643       df1_three_variable tmp;
00644       tmp = sgngam * MAXNUM;
00645       return tmp;
00646 #endif
00647    }
00648 
00649    q = (x - 0.5) * log(x) - x + LS2PI;
00650    if (x > 1.0e8)
00651       return (q);
00652 
00653    p = 1.0 / (x * x);
00654    if (x >= 1000.0)
00655       q += ((7.9365079365079365079365e-4 * p
00656        - 2.7777777777777777777778e-3) * p
00657       + 0.0833333333333333333333) / x;
00658    else
00659       q += polevl(p, A, 4) / x;
00660    return (q);
00661 }
00662 
00663 /*              polevl.c
00664  *              p1evl.c
00665  *
00666  *  Evaluate polynomial
00667  *
00668  *
00669  *
00670  * SYNOPSIS:
00671  *
00672  * int N;
00673  * double x, y, coef[N+1], polevl[];
00674  *
00675  * y = polevl( x, coef, N );
00676  *
00677  *
00678  *
00679  * DESCRIPTION:
00680  *
00681  * Evaluates polynomial of degree N:
00682  *
00683  *                     2          N
00684  * y  =  C  + C x + C x  +...+ C x
00685  *        0    1     2          N
00686  *
00687  * Coefficients are stored in reverse order:
00688  *
00689  * coef[0] = C  , ..., coef[N] = C  .
00690  *            N                   0
00691  *
00692  *  The function p1evl() assumes that coef[N] = 1.0 and is
00693  * omitted from the array.  Its calling arguments are
00694  * otherwise the same as polevl().
00695  *
00696  *
00697  * SPEED:
00698  *
00699  * In the interest of speed, there are no checks for out
00700  * of bounds arithmetic.  This routine is used by most of
00701  * the functions in the library.  Depending on available
00702  * equipment features, the user may wish to rewrite the
00703  * program in microcode or assembly language.
00704  *
00705  */
00706 
00711 static double polevl(double x, void *_coef, int N)
00712 {
00713    double *coef = (double *) (_coef);
00714    double ans;
00715    int i;
00716    double *p;
00717 
00718    p = coef;
00719    ans = *p++;
00720    i = N;
00721 
00722    do
00723       ans = ans * x + *p++;
00724    while (--i);
00725 
00726    return (ans);
00727 }
00728 
00733 static df1_three_variable polevl(const df1_three_variable & x, void *_coef,
00734          int N)
00735 {
00736    double *coef = (double *) (_coef);
00737    df1_three_variable ans;
00738    int i;
00739    double *p;
00740 
00741    p = coef;
00742    ans = *p++;
00743    i = N;
00744 
00745    do
00746       ans = ans * x + *p++;
00747    while (--i);
00748 
00749    return (ans);
00750 }
00751 
00756 static double p1evl(double x, void *_coef, int N)
00757 {
00758    double *coef = (double *) (_coef);
00759    double ans;
00760    double *p;
00761    int i;
00762 
00763    p = coef;
00764    ans = x + *p++;
00765    i = N - 1;
00766 
00767    do
00768       ans = ans * x + *p++;
00769    while (--i);
00770 
00771    return (ans);
00772 }
00773 
00778 static df1_three_variable p1evl(const df1_three_variable & x, void *_coef,
00779         int N)
00780 {
00781    double *coef = (double *) (_coef);
00782    df1_three_variable ans;
00783    double *p;
00784    int i;
00785 
00786    p = coef;
00787    ans = x + *p++;
00788    i = N - 1;
00789 
00790    do
00791       ans = ans * x + *p++;
00792    while (--i);
00793 
00794    return (ans);
00795 }
00796 
00801 static int mtherr(char *s, int n)
00802 {
00803    cerr << "Error code " << n << "in " << *s << endl;
00804    ad_exit(1);
00805    return 0;
00806 }
00807 
00808 df1_three_variable incbet(const df1_three_variable & _aa,
00809          const df1_three_variable & _bb,
00810          const df1_three_variable & _xx);
00811 
00816 dvariable incbet(const dvariable & _a, const dvariable & _b,
00817      const dvariable & _x)
00818 {
00819    ADUNCONST(dvariable, a) ADUNCONST(dvariable, b) ADUNCONST(dvariable, x)
00820    // these three lines will automatically put in the
00821    // derivatvie glue
00822    init_df1_three_variable vx(x);
00823    init_df1_three_variable va(a);
00824    init_df1_three_variable vb(b);
00825    df1_three_variable vy;
00826    vy = incbet(va, vb, vx);
00827    dvariable z;
00828    z = vy;
00829    return z;
00830 }
00831 
00832   dvariable betai(const dvariable& _a,const dvariable& _b,const dvariable& _x)
00833   {
00834     ADUNCONST(dvariable, a) ADUNCONST(dvariable, b) ADUNCONST(dvariable, x)
00835     return incbet(a,b,x);
00836   }
00837 
00843 double incbcf(double _a, double _b, double _x)
00844 {
00845    double a;
00846    double b;
00847    double x;
00848    a = _a;
00849    b = _b;
00850    x = _x;
00851    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
00852    double k1, k2, k3, k4, k5, k6, k7, k8;
00853    double r, t, ans, thresh;
00854    int n;
00855 
00856    k1 = a;
00857    k2 = a + b;
00858    k3 = a;
00859    k4 = a + 1.0;
00860    k5 = 1.0;
00861    k6 = b - 1.0;
00862    k7 = k4;
00863    k8 = a + 2.0;
00864 
00865    pkm2 = 0.0;
00866    qkm2 = 1.0;
00867    pkm1 = 1.0;
00868    qkm1 = 1.0;
00869    ans = 1.0;
00870    r = 1.0;
00871    n = 0;
00872    thresh = 3.0 * MACHEP;
00873    do
00874    {
00875 
00876       xk = -(x * k1 * k2) / (k3 * k4);
00877       pk = pkm1 + pkm2 * xk;
00878       qk = qkm1 + qkm2 * xk;
00879       pkm2 = pkm1;
00880       pkm1 = pk;
00881       qkm2 = qkm1;
00882       qkm1 = qk;
00883 
00884       xk = (x * k5 * k6) / (k7 * k8);
00885       pk = pkm1 + pkm2 * xk;
00886       qk = qkm1 + qkm2 * xk;
00887       pkm2 = pkm1;
00888       pkm1 = pk;
00889       qkm2 = qkm1;
00890       qkm1 = qk;
00891 
00892       if (qk != 0)
00893    r = pk / qk;
00894       if (r != 0)
00895       {
00896    t = fabs((ans - r) / r);
00897    ans = r;
00898       } else
00899    t = 1.0;
00900 
00901       if (t < thresh)
00902    goto cdone;
00903 
00904       k1 += 1.0;
00905       k2 += 1.0;
00906       k3 += 2.0;
00907       k4 += 2.0;
00908       k5 += 1.0;
00909       k6 -= 1.0;
00910       k7 += 2.0;
00911       k8 += 2.0;
00912 
00913       if ((fabs(qk) + fabs(pk)) > big)
00914       {
00915    pkm2 *= biginv;
00916    pkm1 *= biginv;
00917    qkm2 *= biginv;
00918    qkm1 *= biginv;
00919       }
00920       if ((fabs(qk) < biginv) || (fabs(pk) < biginv))
00921       {
00922    pkm2 *= big;
00923    pkm1 *= big;
00924    qkm2 *= big;
00925    qkm1 *= big;
00926       }
00927    }
00928    while (++n < 300);
00929 
00930  cdone:
00931    return (ans);
00932 }
00933 
00939 double incbd(double _a, double _b, double _x)
00940 {
00941    double a=_a;
00942    double b=_b;
00943    double x=_x;
00944 
00945    double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
00946    double k1, k2, k3, k4, k5, k6, k7, k8;
00947    double r, t, ans, z, thresh;
00948    int n;
00949 
00950    k1 = a;
00951    k2 = b - 1.0;
00952    k3 = a;
00953    k4 = a + 1.0;
00954    k5 = 1.0;
00955    k6 = a + b;
00956    k7 = a + 1.0;;
00957    k8 = a + 2.0;
00958 
00959    pkm2 = 0.0;
00960    qkm2 = 1.0;
00961    pkm1 = 1.0;
00962    qkm1 = 1.0;
00963    z = x / (1.0 - x);
00964    ans = 1.0;
00965    r = 1.0;
00966    n = 0;
00967    thresh = 3.0 * MACHEP;
00968    do
00969    {
00970 
00971       xk = -(z * k1 * k2) / (k3 * k4);
00972       pk = pkm1 + pkm2 * xk;
00973       qk = qkm1 + qkm2 * xk;
00974       pkm2 = pkm1;
00975       pkm1 = pk;
00976       qkm2 = qkm1;
00977       qkm1 = qk;
00978 
00979       xk = (z * k5 * k6) / (k7 * k8);
00980       pk = pkm1 + pkm2 * xk;
00981       qk = qkm1 + qkm2 * xk;
00982       pkm2 = pkm1;
00983       pkm1 = pk;
00984       qkm2 = qkm1;
00985       qkm1 = qk;
00986 
00987       if (qk != 0)
00988    r = pk / qk;
00989       if (r != 0)
00990       {
00991    t = fabs((ans - r) / r);
00992    ans = r;
00993       } else
00994    t = 1.0;
00995 
00996       if (t < thresh)
00997    goto cdone;
00998 
00999       k1 += 1.0;
01000       k2 -= 1.0;
01001       k3 += 2.0;
01002       k4 += 2.0;
01003       k5 += 1.0;
01004       k6 += 1.0;
01005       k7 += 2.0;
01006       k8 += 2.0;
01007 
01008       if ((fabs(qk) + fabs(pk)) > big)
01009       {
01010    pkm2 *= biginv;
01011    pkm1 *= biginv;
01012    qkm2 *= biginv;
01013    qkm1 *= biginv;
01014       }
01015       if ((fabs(qk) < biginv) || (fabs(pk) < biginv))
01016       {
01017    pkm2 *= big;
01018    pkm1 *= big;
01019    qkm2 *= big;
01020    qkm1 *= big;
01021       }
01022    }
01023    while (++n < 300);
01024  cdone:
01025    return (ans);
01026 }
01027 
01028 
01029 
01034 double incbet(double _aa,double _bb,double _xx)
01035 {
01036    double aa;
01037    double bb;
01038    double xx;
01039    aa = _aa;
01040    bb = _bb;
01041    xx = _xx;
01042 
01043    double a, b, t, x, xc, w, y;
01044 
01045    int flag;
01046    double one;
01047    double zero;
01048    one = 1.0;
01049    zero = 0.0;
01050 
01051    if (aa <= 0.0 || bb <= 0.0)
01052       goto domerr;
01053 
01054    if ((xx <= 0.0) || (xx >= 1.0))
01055    {
01056       if (xx == 0.0)
01057    return (zero);
01058       if (xx == 1.0)
01059    return (one);
01060     domerr:
01061       cerr << "incbet DOMAIN " << endl;
01062 
01063       return (zero);
01064    }
01065 
01066    flag = 0;
01067    if ( (bb * xx) <= 1.0 && xx <= 0.95)
01068    {
01069       t = pseries(aa, bb, xx);
01070       goto done;
01071    }
01072 
01073    w = 1.0 - xx;
01074 
01075    // Reverse a and b if x is greater than the mean.
01076    if (xx > aa / (aa + bb))
01077    {
01078       flag = 1;
01079       a = bb;
01080       b = aa;
01081       xc = xx;
01082       x = w;
01083    } else
01084    {
01085       a = aa;
01086       b = bb;
01087       xc = w;
01088       x = xx;
01089    }
01090 
01091    if (flag == 1 && (b * x) <= 1.0 && x <= 0.95)
01092    {
01093       t = pseries(a, b, x);
01094       goto done;
01095    }
01096 
01097    // Choose expansion for better convergence.
01098    y = x * (a + b - 2.0) - (a - 1.0);
01099    if (y < 0.0)
01100       w = incbcf(a, b, x);
01101    else
01102       w = incbd(a, b, x) / xc;
01103 
01104    // Multiply w by the factor
01105    // a      b   _             _     _
01106    // x  (1-x)   | (a+b) / ( a | (a) | (b) ) .
01107 
01108    y = a * log(x);
01109    t = b * log(xc);
01110    if ((a + b) < MAXGAM && fabs(y) < MAXLOG
01111        && fabs(t) < MAXLOG)
01112    {
01113       t = pow(xc, b);
01114       t *= pow(x, a);
01115       t /= a;
01116       t *= w;
01117       t *= gamma(a + b) / (gamma(a) * gamma(b));
01118       goto done;
01119    }
01120    // Resort to logarithms.
01121    y += t + lgam(a + b) - lgam(a) - lgam(b);
01122    y += log(w / a);
01123    if (y < MINLOG)
01124       t = 0.0;
01125    else
01126       t = exp(y);
01127  
01128  done:
01129 
01130    if (flag == 1)
01131    {
01132       if (t <= MACHEP)
01133    t = 1.0 - MACHEP;
01134       else
01135    t = 1.0 - t;
01136    }
01137    return (t);
01138 }
01139 
01140 double betai(double _aa, double _bb, double _xx)
01141 {
01142   double ret = incbet(_aa,_bb,_xx);
01143   return ret;
01144 }
01145 
01146 
01147 /*              igam.c
01148  *
01149  *  Incomplete gamma integral
01150  *
01151  *
01152  *
01153  * SYNOPSIS:
01154  *
01155  * double a, x, y, igam();
01156  *
01157  * y = igam( a, x );
01158  *
01159  * DESCRIPTION:
01160  *
01161  * The function is defined by
01162  *
01163  *                           x
01164  *                            -
01165  *                   1       | |  -t  a-1
01166  *  igam(a,x)  =   -----     |   e   t   dt.
01167  *                  -      | |
01168  *                 | (a)    -
01169  *                           0
01170  *
01171  *
01172  * In this implementation both arguments must be positive.
01173  * The integral is evaluated by either a power series or
01174  * continued fraction expansion, depending on the relative
01175  * values of a and x.
01176  *
01177  * ACCURACY:
01178  *
01179  *                      Relative error:
01180  * arithmetic   domain     # trials      peak         rms
01181  *    IEEE      0,30       200000       3.6e-14     2.9e-15
01182  *    IEEE      0,100      300000       9.9e-14     1.5e-14
01183  */
01184 /*             igamc()
01185  *
01186  *  Complemented incomplete gamma integral
01187  *
01188  *
01189  *
01190  * SYNOPSIS:
01191  *
01192  * double a, x, y, igamc();
01193  *
01194  * y = igamc( a, x );
01195  *
01196  * DESCRIPTION:
01197  *
01198  * The function is defined by
01199  *
01200  *
01201  *  igamc(a,x)   =   1 - igam(a,x)
01202  *
01203  *                            inf.
01204  *                              -
01205  *                     1       | |  -t  a-1
01206  *               =   -----     |   e   t   dt.
01207  *                    -      | |
01208  *                   | (a)    -
01209  *                             x
01210  *
01211  *
01212  * In this implementation both arguments must be positive.
01213  * The integral is evaluated by either a power series or
01214  * continued fraction expansion, depending on the relative
01215  * values of a and x.
01216  *
01217  * ACCURACY:
01218  *
01219  * Tested at random a, x.
01220  *                a         x                      Relative error:
01221  * arithmetic   domain   domain     # trials      peak         rms
01222  *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
01223  *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
01224  */
01225 
01226 static df1_three_variable igam(const df1_three_variable & a,
01227              const df1_three_variable & x);
01228 
01233 df1_three_variable igamc(const df1_three_variable & _a,
01234        const df1_three_variable & _x)
01235 {
01236    ADUNCONST(df1_three_variable, a)
01237       ADUNCONST(df1_three_variable, x)
01238       df1_three_variable ans, ax, c, yc, r, t, y, z;
01239    df1_three_variable pk, pkm1, pkm2, qk, qkm1, qkm2;
01240 
01241    if ((value(x) <= 0) || (value(a) <= 0))
01242    {
01243       df1_three_variable tmp;
01244       tmp = 0.0;
01245       return (tmp);
01246    }
01247 
01248    if ((value(x) < 1.0) || (value(x) < value(a)))
01249       return (1.0 - igam(a, x));
01250 
01251    ax = a * log(x) - x - lgam(a);
01252    if (value(ax) < -MAXLOG)
01253    {
01254       cerr << "igamc UNDERFLOW " << endl;
01255       ad_exit(1);
01256    }
01257    ax = exp(ax);
01258 
01259    /* continued fraction */
01260    y = 1.0 - a;
01261    z = x + y + 1.0;
01262    c = 0.0;
01263    pkm2 = 1.0;
01264    qkm2 = x;
01265    pkm1 = x + 1.0;
01266    qkm1 = z * x;
01267    ans = pkm1 / qkm1;
01268 
01269    do
01270    {
01271       c += 1.0;
01272       y += 1.0;
01273       z += 2.0;
01274       yc = y * c;
01275       pk = pkm1 * z - pkm2 * yc;
01276       qk = qkm1 * z - qkm2 * yc;
01277       if (value(qk) != 0)
01278       {
01279    r = pk / qk;
01280    t = fabs((ans - r) / r);
01281    ans = r;
01282       } else
01283    t = 1.0;
01284       pkm2 = pkm1;
01285       pkm1 = pk;
01286       qkm2 = qkm1;
01287       qkm1 = qk;
01288       if (fabs(value(pk)) > big)
01289       {
01290    pkm2 *= biginv;
01291    pkm1 *= biginv;
01292    qkm2 *= biginv;
01293    qkm1 *= biginv;
01294       }
01295    }
01296    while (t > MACHEP);
01297 
01298    return (ans * ax);
01299 }
01300 
01301   // 
01302   // 
01303   // 
01304   // /* left tail of incomplete gamma function:
01305   //  *
01306   //  *          inf.      k
01307   //  *   a  -x   -       x
01308   //  *  x  e     >   ----------
01309   //  *           -     -
01310   //  *          k=0   | (a+k+1)
01311   //  *
01312   //  */
01313   // 
01314 
01319 df1_three_variable igam(const df1_three_variable & _a,
01320       const df1_three_variable & _x)
01321 {
01322    ADUNCONST(df1_three_variable, a)
01323       ADUNCONST(df1_three_variable, x) df1_three_variable ans, ax, c, r;
01324 
01325    if ((value(x) <= 0) || (value(a) <= 0))
01326    {
01327       df1_three_variable tmp;
01328       tmp = 0.0;
01329       return (tmp);
01330    }
01331 
01332    if ((value(x) > 1.0) && (value(x) > value(a)))
01333       return (1.0 - igamc(a, x));
01334 
01335    /* Compute  x**a * exp(-x) / gamma(a)  */
01336    ax = a * log(x) - x - lgam(a);
01337    if (value(ax) < -MAXLOG)
01338    {
01339       cerr << "igam UNDERFLOW " << endl;
01340       ad_exit(1);
01341    }
01342    ax = exp(ax);
01343 
01344    /* power series */
01345    r = a;
01346    c = 1.0;
01347    ans = 1.0;
01348 
01349    do
01350    {
01351       r += 1.0;
01352       c *= x / r;
01353       ans += c;
01354    }
01355    while (c / ans > MACHEP);
01356 
01357    return (ans * ax / a);
01358 }
01359 
01360 
01361 
01362   /*                                                  incbet.c
01363    *
01364    *  Incomplete beta integral
01365    *
01366    *
01367    * SYNOPSIS:
01368    *
01369    * double a, b, x, y, incbet();
01370    *
01371    * y = incbet( a, b, x );
01372    *
01373    *
01374    * DESCRIPTION:
01375    *
01376    * Returns incomplete beta integral of the arguments, evaluated
01377    * from zero to x.  The function is defined as
01378    *
01379    *                  x
01380    *     -            -
01381    *    | (a+b)      | |  a-1     b-1
01382    *  -----------    |   t   (1-t)   dt.
01383    *   -     -     | |
01384    *  | (a) | (b)   -
01385    *                 0
01386    *
01387    * The domain of definition is 0 <= x <= 1.  In this
01388    * implementation a and b are restricted to positive values.
01389    * The integral from x to 1 may be obtained by the symmetry
01390    * relation
01391    *
01392    *    1 - incbet( a, b, x )  =  incbet( b, a, 1-x ).
01393    *
01394    * The integral is evaluated by a continued fraction expansion
01395    * or, when b*x is small, by a power series.
01396    *
01397    * ACCURACY:
01398    *
01399    * Tested at uniformly distributed random points (a,b,x) with a and b
01400    * in "domain" and x between 0 and 1.
01401    *                                        Relative error
01402    * arithmetic   domain     # trials      peak         rms
01403    *    IEEE      0,5         10000       6.9e-15     4.5e-16
01404    *    IEEE      0,85       250000       2.2e-13     1.7e-14
01405    *    IEEE      0,1000      30000       5.3e-12     6.3e-13
01406    *    IEEE      0,10000    250000       9.3e-11     7.1e-12
01407    *    IEEE      0,100000    10000       8.7e-10     4.8e-11
01408    * Outputs smaller than the IEEE gradual underflow threshold
01409    * were excluded from these statistics.
01410    *
01411    * ERROR MESSAGES:
01412    *   message         condition      value returned
01413    * incbet domain      x<0, x>1          0.0
01414    * incbet underflow                     0.0
01415    */
01416 
01417 static df1_three_variable incbcf(const df1_three_variable & _a,
01418          const df1_three_variable & _b,
01419          const df1_three_variable & _x);
01420 static df1_three_variable pseries(const df1_three_variable & _a,
01421           const df1_three_variable & _b,
01422           const df1_three_variable & _x);
01423 static df1_three_variable incbd(const df1_three_variable & _a,
01424         const df1_three_variable & _b,
01425         const df1_three_variable & _x);
01430 df1_three_variable incbet(const df1_three_variable & _aa,
01431          const df1_three_variable & _bb,
01432          const df1_three_variable & _xx)
01433 {
01434 
01435    df1_three_variable aa;
01436    df1_three_variable bb;
01437    df1_three_variable xx;
01438    aa = _aa;
01439    bb = _bb;
01440    xx = _xx;
01441 
01442    df1_three_variable a, b, t, x, xc, w, y;
01443 
01444    int flag;
01445    df1_three_variable one;
01446    df1_three_variable zero;
01447    one = 1.0;
01448    zero = 0.0;
01449 
01450    if (value(aa) <= 0.0 || value(bb) <= 0.0)
01451       goto domerr;
01452 
01453    if ((value(xx) <= 0.0) || (value(xx) >= 1.0))
01454    {
01455       if (value(xx) == 0.0)
01456    return (zero);
01457       if (value(xx) == 1.0)
01458    return (one);
01459     domerr:
01460       cerr << "incbet DOMAIN " << endl;
01461 
01462       return (zero);
01463    }
01464 
01465    flag = 0;
01466    if (value(bb * xx) <= 1.0 && value(xx) <= 0.95)
01467    {
01468       t = pseries(aa, bb, xx);
01469       goto done;
01470    }
01471 
01472    w = 1.0 - xx;
01473 
01474    /* Reverse a and b if x is greater than the mean. */
01475    if (value(xx) > (value(aa) / value(aa + bb)))
01476    {
01477       flag = 1;
01478       a = bb;
01479       b = aa;
01480       xc = xx;
01481       x = w;
01482    } else
01483    {
01484       a = aa;
01485       b = bb;
01486       xc = w;
01487       x = xx;
01488    }
01489 
01490    if (flag == 1 && value((b * x)) <= 1.0 && value(x) <= 0.95)
01491    {
01492       t = pseries(a, b, x);
01493       goto done;
01494    }
01495 
01496    /* Choose expansion for better convergence. */
01497    y = x * (a + b - 2.0) - (a - 1.0);
01498    if (value(y) < 0.0)
01499       w = incbcf(a, b, x);
01500    else
01501       w = incbd(a, b, x) / xc;
01502 
01503    /* Multiply w by the factor
01504       a      b   _             _     _
01505       x  (1-x)   | (a+b) / ( a | (a) | (b) ) .   */
01506 
01507    y = a * log(x);
01508    t = b * log(xc);
01509    if (value((a + b)) < MAXGAM && fabs(value(y)) < MAXLOG
01510        && fabs(value(t)) < MAXLOG)
01511    {
01512       t = pow(xc, b);
01513       t *= pow(x, a);
01514       t /= a;
01515       t *= w;
01516       t *= gamma(a + b) / (gamma(a) * gamma(b));
01517       goto done;
01518    }
01519    /* Resort to logarithms.  */
01520    y += t + lgam(a + b) - lgam(a) - lgam(b);
01521    y += log(w / a);
01522    if (value(y) < MINLOG)
01523       t = 0.0;
01524    else
01525       t = exp(y);
01526 
01527  done:
01528 
01529    if (flag == 1)
01530    {
01531       if (value(t) <= MACHEP)
01532    t = 1.0 - MACHEP;
01533       else
01534    t = 1.0 - t;
01535    }
01536    return (t);
01537 }
01538 
01544 static df1_three_variable incbcf(const df1_three_variable & _a,
01545          const df1_three_variable & _b,
01546          const df1_three_variable & _x)
01547 {
01548    ADUNCONST(df1_three_variable, a)
01549       ADUNCONST(df1_three_variable, b)
01550       ADUNCONST(df1_three_variable, x)
01551       df1_three_variable xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
01552    df1_three_variable k1, k2, k3, k4, k5, k6, k7, k8;
01553    df1_three_variable r, t, ans, thresh;
01554    int n;
01555 
01556    k1 = a;
01557    k2 = a + b;
01558    k3 = a;
01559    k4 = a + 1.0;
01560    k5 = 1.0;
01561    k6 = b - 1.0;
01562    k7 = k4;
01563    k8 = a + 2.0;
01564 
01565    pkm2 = 0.0;
01566    qkm2 = 1.0;
01567    pkm1 = 1.0;
01568    qkm1 = 1.0;
01569    ans = 1.0;
01570    r = 1.0;
01571    n = 0;
01572    thresh = 3.0 * MACHEP;
01573    do
01574    {
01575 
01576       xk = -(x * k1 * k2) / (k3 * k4);
01577       pk = pkm1 + pkm2 * xk;
01578       qk = qkm1 + qkm2 * xk;
01579       pkm2 = pkm1;
01580       pkm1 = pk;
01581       qkm2 = qkm1;
01582       qkm1 = qk;
01583 
01584       xk = (x * k5 * k6) / (k7 * k8);
01585       pk = pkm1 + pkm2 * xk;
01586       qk = qkm1 + qkm2 * xk;
01587       pkm2 = pkm1;
01588       pkm1 = pk;
01589       qkm2 = qkm1;
01590       qkm1 = qk;
01591 
01592       if (value(qk) != 0)
01593    r = pk / qk;
01594       if (value(r) != 0)
01595       {
01596    t = fabs((ans - r) / r);
01597    ans = r;
01598       } else
01599    t = 1.0;
01600 
01601       if (value(t) < value(thresh))
01602    goto cdone;
01603 
01604       k1 += 1.0;
01605       k2 += 1.0;
01606       k3 += 2.0;
01607       k4 += 2.0;
01608       k5 += 1.0;
01609       k6 -= 1.0;
01610       k7 += 2.0;
01611       k8 += 2.0;
01612 
01613       if ((value(fabs(qk)) + value(fabs(pk))) > big)
01614       {
01615    pkm2 *= biginv;
01616    pkm1 *= biginv;
01617    qkm2 *= biginv;
01618    qkm1 *= biginv;
01619       }
01620       if ((value(fabs(qk)) < biginv) || (value(fabs(pk)) < biginv))
01621       {
01622    pkm2 *= big;
01623    pkm1 *= big;
01624    qkm2 *= big;
01625    qkm1 *= big;
01626       }
01627    }
01628    while (++n < 300);
01629 
01630  cdone:
01631    return (ans);
01632 }
01633 
01639 static df1_three_variable incbd(const df1_three_variable & _a,
01640         const df1_three_variable & _b,
01641         const df1_three_variable & _x)
01642 {
01643    ADUNCONST(df1_three_variable, a)
01644       ADUNCONST(df1_three_variable, b) ADUNCONST(df1_three_variable, x)
01645       //df1_three_variable s, t, u, v, n, t1, z, ai;
01646    df1_three_variable xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
01647    df1_three_variable k1, k2, k3, k4, k5, k6, k7, k8;
01648    df1_three_variable r, t, ans, z, thresh;
01649    int n;
01650 
01651    k1 = a;
01652    k2 = b - 1.0;
01653    k3 = a;
01654    k4 = a + 1.0;
01655    k5 = 1.0;
01656    k6 = a + b;
01657    k7 = a + 1.0;;
01658    k8 = a + 2.0;
01659 
01660    pkm2 = 0.0;
01661    qkm2 = 1.0;
01662    pkm1 = 1.0;
01663    qkm1 = 1.0;
01664    z = x / (1.0 - x);
01665    ans = 1.0;
01666    r = 1.0;
01667    n = 0;
01668    thresh = 3.0 * MACHEP;
01669    do
01670    {
01671 
01672       xk = -(z * k1 * k2) / (k3 * k4);
01673       pk = pkm1 + pkm2 * xk;
01674       qk = qkm1 + qkm2 * xk;
01675       pkm2 = pkm1;
01676       pkm1 = pk;
01677       qkm2 = qkm1;
01678       qkm1 = qk;
01679 
01680       xk = (z * k5 * k6) / (k7 * k8);
01681       pk = pkm1 + pkm2 * xk;
01682       qk = qkm1 + qkm2 * xk;
01683       pkm2 = pkm1;
01684       pkm1 = pk;
01685       qkm2 = qkm1;
01686       qkm1 = qk;
01687 
01688       if (value(qk) != 0)
01689    r = pk / qk;
01690       if (value(r) != 0)
01691       {
01692    t = fabs((ans - r) / r);
01693    ans = r;
01694       } else
01695    t = 1.0;
01696 
01697       if (value(t) < value(thresh))
01698    goto cdone;
01699 
01700       k1 += 1.0;
01701       k2 -= 1.0;
01702       k3 += 2.0;
01703       k4 += 2.0;
01704       k5 += 1.0;
01705       k6 += 1.0;
01706       k7 += 2.0;
01707       k8 += 2.0;
01708 
01709       if ((value(fabs(qk)) + value(fabs(pk))) > big)
01710       {
01711    pkm2 *= biginv;
01712    pkm1 *= biginv;
01713    qkm2 *= biginv;
01714    qkm1 *= biginv;
01715       }
01716       if ((value(fabs(qk)) < biginv) || (value(fabs(pk)) < biginv))
01717       {
01718    pkm2 *= big;
01719    pkm1 *= big;
01720    qkm2 *= big;
01721    qkm1 *= big;
01722       }
01723    }
01724    while (++n < 300);
01725  cdone:
01726    return (ans);
01727 }
01728 
01734 df1_three_variable pseries(const df1_three_variable & _a,
01735          const df1_three_variable & _b,
01736          const df1_three_variable & _x)
01737 {
01738    df1_three_variable a;
01739    df1_three_variable b;
01740    df1_three_variable x;
01741    a = _a;
01742    b = _b;
01743    x = _x;
01744    df1_three_variable s, t, u, v, n, t1, z, ai;
01745 
01746    ai = 1.0 / a;
01747    u = (1.0 - b) * x;
01748    v = u / (a + 1.0);
01749    t1 = v;
01750    t = u;
01751    n = 2.0;
01752    s = 0.0;
01753    z = MACHEP * ai;
01754    while (value(fabs(v)) > value(z))
01755    {
01756       u = (n - b) * x / n;
01757       t *= u;
01758       v = t / (a + n);
01759       s += v;
01760       n += 1.0;
01761    }
01762    s += t1;
01763    s += ai;
01764 
01765    u = a * log(x);
01766    if (value((a + b)) < MAXGAM && value(fabs(u)) < MAXLOG)
01767    {
01768       gamma(a + b);
01769       gamma(a);
01770       gamma(b);
01771       gamma(a) * gamma(b);
01772       t = gamma(a + b) / (gamma(a) * gamma(b));
01773       s = s * t * pow(x, a);
01774    } else
01775    {
01776       t = lgam(a + b) - lgam(a) - lgam(b) + u + log(s);
01777       if (value(t) < MINLOG)
01778    s = 0.0;
01779       else
01780    s = exp(t);
01781    }
01782    return (s);
01783 }
01784 
01789 static double get_values(double a, double b, double x)
01790 {
01791    df1_three_variable va, vb, vx;
01792    va.initialize();
01793    vb.initialize();
01794    vx.initialize();
01795    va = a;
01796    vb = b;
01797    vx = x;
01798    *va.get_u_x() = 1.0;
01799    *vb.get_u_y() = 1.0;
01800    *vx.get_u_z() = 1.0;
01801    df1_three_variable vy = incbet(va, vb, vx);
01802    return *vy.get_u();
01803 }
01804 
01809 static df1_three_variable df3_get_values(double a, double b, double x)
01810 {
01811    df1_three_variable va, vb, vx;
01812    va.initialize();
01813    vb.initialize();
01814    vx.initialize();
01815    va = a;
01816    vb = b;
01817    vx = x;
01818    *va.get_u_x() = 1.0;
01819    *vb.get_u_y() = 1.0;
01820    *vx.get_u_z() = 1.0;
01821    df1_three_variable vy = incbet(va, vb, vx);
01822    return vy;
01823 }
01824 
01830 double pseries(const double & _a, const double & _b, const double & _x)
01831 {
01832    double a=_a;
01833    double b=_b;
01834    double x=_x;
01835 
01836    double s, t, u, v, n, t1, z, ai;
01837 
01838    ai = 1.0 / a;
01839    u = (1.0 - b) * x;
01840    v = u / (a + 1.0);
01841    t1 = v;
01842    t = u;
01843    n = 2.0;
01844    s = 0.0;
01845    z = MACHEP * ai;
01846    while (fabs(v) > z)
01847    {
01848       u = (n - b) * x / n;
01849       t *= u;
01850       v = t / (a + n);
01851       s += v;
01852       n += 1.0;
01853    }
01854    s += t1;
01855    s += ai;
01856 
01857    //u = a * log(x);
01858    if ((a + b) < MAXGAM && fabs(u) < MAXLOG)
01859    {
01860       //gamma(a + b);
01861       //gamma(a);
01862       //gamma(b);
01863       //gamma(a) * gamma(b);
01864       t = gamma(a + b) / (gamma(a) * gamma(b));
01865       s = s * t * pow(x, a);
01866    }/* else
01867    {
01868       t = lgam(a + b) - lgam(a) - lgam(b) + u + log(s);
01869       if (t < MINLOG)
01870    s = 0.0;
01871       else
01872    s = exp(t);
01873    }*/
01874    return (s);
01875 }
01876