lib/libm/ivf.c

   1 /*                                                      ivf.c
   2  *
   3  *      Modified Bessel function of noninteger order
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float v, x, y, ivf();
  10  *
  11  * y = ivf( v, x );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Returns modified Bessel function of order v of the
  18  * argument.  If x is negative, v must be integer valued.
  19  *
  20  * The function is defined as Iv(x) = Jv( ix ).  It is
  21  * here computed in terms of the confluent hypergeometric
  22  * function, according to the formula
  23  *
  24  *              v  -x
  25  * Iv(x) = (x/2)  e   hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)
  26  *
  27  * If v is a negative integer, then v is replaced by -v.
  28  *
  29  *
  30  * ACCURACY:
  31  *
  32  * Tested at random points (v, x), with v between 0 and
  33  * 30, x between 0 and 28.
  34  * arithmetic   domain     # trials      peak         rms
  35  *                      Relative error:
  36  *    IEEE      0,15          3000      4.7e-6      5.4e-7
  37  *          Absolute error (relative when function > 1)
  38  *    IEEE      0,30          5000      8.5e-6      1.3e-6
  39  *
  40  * Accuracy is diminished if v is near a negative integer.
  41  * The useful domain for relative error is limited by overflow
  42  * of the single precision exponential function.
  43  *
  44  * See also hyperg.c.
  45  *
  46  */
  47 \f/*                                                     iv.c    */
  48 /*      Modified Bessel function of noninteger order            */
  49 /* If x < 0, then v must be an integer. */
  50
  51
  52 /*
  53 Cephes Math Library Release 2.2: June, 1992
  54 Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  55 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  56 */
  57
  58
  59 #include "mconf.h"
  60
  61 extern double MAXNUMF;
  62 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  63
  64 #ifdef ANSIC
  65 float hypergf(float, float, float);
  66 float expf(float), gammaf(float), logf(float), floorf(float);
  67
  68 float ivf( float v, float x )
  69 #else
  70 float hypergf(), expf(), gammaf(), logf(), floorf();
  71
  72 float ivf( v, x )
  73 double v, x;
  74 #endif
  75 {
  76 int sign;
  77 float t, ax;
  78
  79 /* If v is a negative integer, invoke symmetry */
  80 t = floorf(v);
  81 if( v < 0.0 )
  82         {
  83         if( t == v )
  84                 {
  85                 v = -v; /* symmetry */
  86                 t = -t;
  87                 }
  88         }
  89 /* If x is negative, require v to be an integer */
  90 sign = 1;
  91 if( x < 0.0 )
  92         {
  93         if( t != v )
  94                 {
  95                 mtherr( "ivf", DOMAIN );
  96                 return( 0.0 );
  97                 }
  98         if( v != 2.0 * floorf(v/2.0) )
  99                 sign = -1;
 100         }
 101
 102 /* Avoid logarithm singularity */
 103 if( x == 0.0 )
 104         {
 105         if( v == 0.0 )
 106                 return( 1.0 );
 107         if( v < 0.0 )
 108                 {
 109                 mtherr( "ivf", OVERFLOW );
 110                 return( MAXNUMF );
 111                 }
 112         else
 113                 return( 0.0 );
 114         }
 115
 116 ax = fabsf(x);
 117 t = v * logf( 0.5 * ax )  -  x;
 118 t = sign * expf(t) / gammaf( v + 1.0 );
 119 ax = v + 0.5;
 120 return( t * hypergf( ax,  2.0 * ax,  2.0 * x ) );
 121 }