lib/libm/betaf.c

   1 /*                                                      betaf.c
   2  *
   3  *      Beta function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float a, b, y, betaf();
  10  *
  11  * y = betaf( 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  *    IEEE       0,30       10000       4.0e-5      6.0e-6
  33  *    IEEE       -20,0      10000       4.9e-3      5.4e-5
  34  *
  35  * ERROR MESSAGES:
  36  *
  37  *   message         condition          value returned
  38  * betaf 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.2:  July, 1992
  48 Copyright 1984, 1987 by Stephen L. Moshier
  49 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  50 */
  51
  52 #include "mconf.h"
  53
  54 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  55
  56 #define MAXGAM 34.84425627277176174
  57
  58
  59 extern float MAXLOGF, MAXNUMF;
  60 extern int sgngamf;
  61
  62 #ifdef ANSIC
  63 float gammaf(float), lgamf(float), expf(float), floorf(float);
  64 #else
  65 float gammaf(), lgamf(), expf(), floorf();
  66 #endif
  67
  68 #ifdef ANSIC
  69 float betaf( float aa, float bb )
  70 #else
  71 float betaf( aa, bb )
  72 double aa, bb;
  73 #endif
  74 {
  75 float a, b, y;
  76 int sign;
  77
  78 sign = 1;
  79 a = aa;
  80 b = bb;
  81 if( a <= 0.0 )
  82         {
  83         if( a == floorf(a) )
  84                 goto over;
  85         }
  86 if( b <= 0.0 )
  87         {
  88         if( b == floorf(b) )
  89                 goto over;
  90         }
  91
  92
  93 y = a + b;
  94 if( fabsf(y) > MAXGAM )
  95         {
  96         y = lgamf(y);
  97         sign *= sgngamf; /* keep track of the sign */
  98         y = lgamf(b) - y;
  99         sign *= sgngamf;
 100         y = lgamf(a) + y;
 101         sign *= sgngamf;
 102         if( y > MAXLOGF )
 103                 {
 104 over:
 105                 mtherr( "betaf", OVERFLOW );
 106                 return( sign * MAXNUMF );
 107                 }
 108         return( sign * expf(y) );
 109         }
 110
 111 y = gammaf(y);
 112 if( y == 0.0 )
 113         goto over;
 114
 115 if( a > b )
 116         {
 117         y = gammaf(a)/y;
 118         y *= gammaf(b);
 119         }
 120 else
 121         {
 122         y = gammaf(b)/y;
 123         y *= gammaf(a);
 124         }
 125
 126 return(y);
 127 }