lib/libm/iv.c

   1 /*                                                      iv.c
   2  *
   3  *      Modified Bessel function of noninteger order
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double v, x, y, iv();
  10  *
  11  * y = iv( 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  *                      Relative error:
  35  * arithmetic   domain     # trials      peak         rms
  36  *    DEC       0,30          2000      3.1e-15     5.4e-16
  37  *    IEEE      0,30         10000      1.7e-14     2.7e-15
  38  *
  39  * Accuracy is diminished if v is near a negative integer.
  40  *
  41  * See also hyperg.c.
  42  *
  43  */
  44 \f/*                                                     iv.c    */
  45 /*      Modified Bessel function of noninteger order            */
  46 /* If x < 0, then v must be an integer. */
  47
  48
  49 /*
  50 Cephes Math Library Release 2.8:  June, 2000
  51 Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
  52 */
  53
  54
  55 #include "mconf.h"
  56 #ifdef ANSIPROT
  57 extern double hyperg ( double, double, double );
  58 extern double exp ( double );
  59 extern double gamma ( double );
  60 extern double log ( double );
  61 extern double fabs ( double );
  62 extern double floor ( double );
  63 #else
  64 double hyperg(), exp(), gamma(), log(), fabs(), floor();
  65 #endif
  66 extern double MACHEP, MAXNUM;
  67
  68 double iv( v, x )
  69 double v, x;
  70 {
  71 int sign;
  72 double t, ax;
  73
  74 /* If v is a negative integer, invoke symmetry */
  75 t = floor(v);
  76 if( v < 0.0 )
  77         {
  78         if( t == v )
  79                 {
  80                 v = -v; /* symmetry */
  81                 t = -t;
  82                 }
  83         }
  84 /* If x is negative, require v to be an integer */
  85 sign = 1;
  86 if( x < 0.0 )
  87         {
  88         if( t != v )
  89                 {
  90                 mtherr( "iv", DOMAIN );
  91                 return( 0.0 );
  92                 }
  93         if( v != 2.0 * floor(v/2.0) )
  94                 sign = -1;
  95         }
  96
  97 /* Avoid logarithm singularity */
  98 if( x == 0.0 )
  99         {
 100         if( v == 0.0 )
 101                 return( 1.0 );
 102         if( v < 0.0 )
 103                 {
 104                 mtherr( "iv", OVERFLOW );
 105                 return( MAXNUM );
 106                 }
 107         else
 108                 return( 0.0 );
 109         }
 110
 111 ax = fabs(x);
 112 t = v * log( 0.5 * ax )  -  x;
 113 t = sign * exp(t) / gamma( v + 1.0 );
 114 ax = v + 0.5;
 115 return( t * hyperg( ax,  2.0 * ax,  2.0 * x ) );
 116 }