libm/double/beta.c

   1 /*                                                      beta.c
   2  *
   3  *      Beta function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double a, b, y, beta();
  10  *
  11  * y = beta( a, b );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  *                   -     -
  18  *                  | (a) | (b)
  19  * beta( a, b )  =  -----------.
  20  *                     -
  21  *                    | (a+b)
  22  *
  23  * For large arguments the logarithm of the function is
  24  * evaluated using lgam(), then exponentiated.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    DEC        0,30        1700       7.7e-15     1.5e-15
  33  *    IEEE       0,30       30000       8.1e-14     1.1e-14
  34  *
  35  * ERROR MESSAGES:
  36  *
  37  *   message         condition          value returned
  38  * beta overflow    log(beta) > MAXLOG       0.0
  39  *                  a or b <0 integer        0.0
  40  *
  41  */
  42 \f
  43 /*                                                      beta.c  */
  44
  45
  46 /*
  47 Cephes Math Library Release 2.0:  April, 1987
  48 Copyright 1984, 1987 by Stephen L. Moshier
  49 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50 */
  51
  52 #include <math.h>
  53
  54 #ifdef UNK
  55 #define MAXGAM 34.84425627277176174
  56 #endif
  57 #ifdef DEC
  58 #define MAXGAM 34.84425627277176174
  59 #endif
  60 #ifdef IBMPC
  61 #define MAXGAM 171.624376956302725
  62 #endif
  63 #ifdef MIEEE
  64 #define MAXGAM 171.624376956302725
  65 #endif
  66
  67 #ifdef ANSIPROT
  68 extern double fabs ( double );
  69 extern double gamma ( double );
  70 extern double lgam ( double );
  71 extern double exp ( double );
  72 extern double log ( double );
  73 extern double floor ( double );
  74 #else
  75 double fabs(), gamma(), lgam(), exp(), log(), floor();
  76 #endif
  77 extern double MAXLOG, MAXNUM;
  78 extern int sgngam;
  79
  80 double beta( a, b )
  81 double a, b;
  82 {
  83 double y;
  84 int sign;
  85
  86 sign = 1;
  87
  88 if( a <= 0.0 )
  89         {
  90         if( a == floor(a) )
  91                 goto over;
  92         }
  93 if( b <= 0.0 )
  94         {
  95         if( b == floor(b) )
  96                 goto over;
  97         }
  98
  99
 100 y = a + b;
 101 if( fabs(y) > MAXGAM )
 102         {
 103         y = lgam(y);
 104         sign *= sgngam; /* keep track of the sign */
 105         y = lgam(b) - y;
 106         sign *= sgngam;
 107         y = lgam(a) + y;
 108         sign *= sgngam;
 109         if( y > MAXLOG )
 110                 {
 111 over:
 112                 mtherr( "beta", OVERFLOW );
 113                 return( sign * MAXNUM );
 114                 }
 115         return( sign * exp(y) );
 116         }
 117
 118 y = gamma(y);
 119 if( y == 0.0 )
 120         goto over;
 121
 122 if( a > b )
 123         {
 124         y = gamma(a)/y;
 125         y *= gamma(b);
 126         }
 127 else
 128         {
 129         y = gamma(b)/y;
 130         y *= gamma(a);
 131         }
 132
 133 return(y);
 134 }
 135
 136
 137
 138 /* Natural log of |beta|.  Return the sign of beta in sgngam.  */
 139
 140 double lbeta( a, b )
 141 double a, b;
 142 {
 143 double y;
 144 int sign;
 145
 146 sign = 1;
 147
 148 if( a <= 0.0 )
 149         {
 150         if( a == floor(a) )
 151                 goto over;
 152         }
 153 if( b <= 0.0 )
 154         {
 155         if( b == floor(b) )
 156                 goto over;
 157         }
 158
 159
 160 y = a + b;
 161 if( fabs(y) > MAXGAM )
 162         {
 163         y = lgam(y);
 164         sign *= sgngam; /* keep track of the sign */
 165         y = lgam(b) - y;
 166         sign *= sgngam;
 167         y = lgam(a) + y;
 168         sign *= sgngam;
 169         sgngam = sign;
 170         return( y );
 171         }
 172
 173 y = gamma(y);
 174 if( y == 0.0 )
 175         {
 176 over:
 177         mtherr( "lbeta", OVERFLOW );
 178         return( sign * MAXNUM );
 179         }
 180
 181 if( a > b )
 182         {
 183         y = gamma(a)/y;
 184         y *= gamma(b);
 185         }
 186 else
 187         {
 188         y = gamma(b)/y;
 189         y *= gamma(a);
 190         }
 191
 192 if( y < 0 )
 193   {
 194     sgngam = -1;
 195     y = -y;
 196   }
 197 else
 198   sgngam = 1;
 199
 200 return( log(y) );
 201 }