/*************************************************************************\
* Copyright (c) 2002 The University of Chicago, as Operator of Argonne
* National Laboratory.
* Copyright (c) 2002 The Regents of the University of California, as
* Operator of Los Alamos National Laboratory.
* This file is distributed subject to a Software License Agreement found
* in the file LICENSE that is included with this distribution. 
\*************************************************************************/

/* file: betai.c
 * purpose: routines for incomplete beta function
 * M. Borland, 1995.
 * Based on Abramowitz and Stegun, 26.5.8; and Numerical Recipes
 $Log: betai.c,v $
 Revision 1.10  2003/08/28 19:56:25  soliday
 Cleaned up some of the code.

 Revision 1.9  2002/08/14 16:18:53  soliday
 Added Open License

 Revision 1.8  1999/12/15 16:29:53  soliday
 Added support for vxWorks.

 Revision 1.7  1999/07/12 20:37:31  soliday
 Added call to fdlibm library for lgamma function on WIN32

 Revision 1.6  1999/07/08 22:11:53  soliday
 Added "lgamma function not implemented on WIN32" error message

 Revision 1.5  1999/05/28 14:58:58  soliday
 Removed compiler warnings under Linux.

 Revision 1.4  1996/09/27 20:48:43  borland
 Commented-out warning statement about nonconvergence of sums.

 * Revision 1.3  1995/09/06  14:13:36  saunders
 * First test release of SDDS1.5
 *
 */
#include "mdb.h"

#if defined(_WIN32)
#include "fdlibm.h"
#endif

double betaIncSum(double x, double a, double b);
double betaInc1(double x, double a, double b);

#define BETAI_ACCURACY 1e-10

double betaComp(double a, double b)
{
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented on vxWorks\n");
  exit(1);
#else
  return exp(lgamma(a)+lgamma(b)-lgamma(a+b));
#endif
}

double lnBetaComp(double a, double b)
{
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented on vxWorks\n");
  exit(1);
#else
  return lgamma(a)+lgamma(b)-lgamma(a+b);
#endif
}

double betaInc(double a, double b, double x)
{
    double xLimit, sum, lnBab=0;
    short swapped = 0;

    if (x<0 || x>1)
        return -1.0;
    if (a+b!=2) {
        xLimit = (a + 1)/(a + b - 2);
        if (x>xLimit) {
            swapped = 1;
            x = 1 - x;
            SWAP_DOUBLE(a, b);
            }
        }
    lnBab = lnBetaComp(a, b);
/*    sum = pow(x, a)*pow(1-x, b)*betaIncSum(a, b, x)/(Bab*a); */
    sum = exp(a*log(x)+b*log(1-x)-lnBab)*betaIncSum(a, b, x)/a;
    if (!swapped) 
        return sum;
    return 1 - sum;
    }

#define MAXIMUM_ITERATIONS 10000

double betaIncSum(double a, double b, double x)
{
    double f1, f2;
    double A0, B0, A_1, B_1, A1, B1, d;
    double m, aM1, aP1, aPb, mT2;

    /* evaluate the first step (d1) outside the loop to get started */
    A_1 = B_1 = A0 = 1;
    B0 = 1 - (a+b)/(a+1)*x;
    aPb = a+b;
    aP1 = a+1;
    aM1 = a-1;
    m = 1;
    do {
        mT2 = m*2;
        d = m*(b-m)*x/( (aM1+mT2)*(a+mT2) );
        A1 = A0 + d*A_1;
        B1 = B0 + d*B_1;
        f1 = A1/B1;
        A_1 = A0;
        B_1 = B0;
        A0  = A1;
        B0  = B1;
        d = -(a+m)*(aPb+m)*x/( (a+mT2)*(aP1+mT2) );
        A1 = A0 + d*A_1;
        B1 = B0 + d*B_1;
        f2 = A1/B1;
        A_1 = A0;
        B_1 = B0;
        A0  = A1;
        B0  = B1;
        if (B1) {
            /* rescale the recurrence to avoid over/underflow */
            A_1 /= B1;
            B_1 /= B1;
            A0 /= B1;
            B0 = 1;
            }
        m++;
        } while (fabs(f1-f2)>BETAI_ACCURACY && m<MAXIMUM_ITERATIONS);
/*
    if (m>=MAXIMUM_ITERATIONS)
        fprintf(stderr, "warning: too many iterations for incomplete beta function sum (%e, %e, %e) (betaInc)\n",
                a, b, x);
*/
    return f2;
    }

/* This isn't used, as it isn't as efficient as betaInc(), but it sure is simple! */
double betaInc1(double a, double b, double x)
{
    double xp, sum, term;
    long i;
    xp = x;
    i = sum = 0;
    do {
        sum += (term=betaComp(a+1, i+1.0)/betaComp(a+b, i+1.0)*xp);
        xp = x*xp;
        i++;
        } while (term>BETAI_ACCURACY);
    return (sum+1)*pow(x, a)*pow(1-x, b)/(a*betaComp(a, b));
    }