| 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 | /* 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 | } |
|---|
