diff options
author | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
---|---|---|
committer | Eric Andersen <andersen@codepoet.org> | 2001-11-22 14:04:29 +0000 |
commit | 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch) | |
tree | 3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/jn.c | |
parent | c117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff) |
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD).
-Erik
Diffstat (limited to 'libm/double/jn.c')
-rw-r--r-- | libm/double/jn.c | 133 |
1 files changed, 0 insertions, 133 deletions
diff --git a/libm/double/jn.c b/libm/double/jn.c deleted file mode 100644 index ee05395aa..000000000 --- a/libm/double/jn.c +++ /dev/null @@ -1,133 +0,0 @@ -/* jn.c - * - * Bessel function of integer order - * - * - * - * SYNOPSIS: - * - * int n; - * double x, y, jn(); - * - * y = jn( n, x ); - * - * - * - * DESCRIPTION: - * - * Returns Bessel function of order n, where n is a - * (possibly negative) integer. - * - * The ratio of jn(x) to j0(x) is computed by backward - * recurrence. First the ratio jn/jn-1 is found by a - * continued fraction expansion. Then the recurrence - * relating successive orders is applied until j0 or j1 is - * reached. - * - * If n = 0 or 1 the routine for j0 or j1 is called - * directly. - * - * - * - * ACCURACY: - * - * Absolute error: - * arithmetic range # trials peak rms - * DEC 0, 30 5500 6.9e-17 9.3e-18 - * IEEE 0, 30 5000 4.4e-16 7.9e-17 - * - * - * Not suitable for large n or x. Use jv() instead. - * - */ - -/* jn.c -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 2000 by Stephen L. Moshier -*/ -#include <math.h> -#ifdef ANSIPROT -extern double fabs ( double ); -extern double j0 ( double ); -extern double j1 ( double ); -#else -double fabs(), j0(), j1(); -#endif -extern double MACHEP; - -double jn( n, x ) -int n; -double x; -{ -double pkm2, pkm1, pk, xk, r, ans; -int k, sign; - -if( n < 0 ) - { - n = -n; - if( (n & 1) == 0 ) /* -1**n */ - sign = 1; - else - sign = -1; - } -else - sign = 1; - -if( x < 0.0 ) - { - if( n & 1 ) - sign = -sign; - x = -x; - } - -if( n == 0 ) - return( sign * j0(x) ); -if( n == 1 ) - return( sign * j1(x) ); -if( n == 2 ) - return( sign * (2.0 * j1(x) / x - j0(x)) ); - -if( x < MACHEP ) - return( 0.0 ); - -/* continued fraction */ -#ifdef DEC -k = 56; -#else -k = 53; -#endif - -pk = 2 * (n + k); -ans = pk; -xk = x * x; - -do - { - pk -= 2.0; - ans = pk - (xk/ans); - } -while( --k > 0 ); -ans = x/ans; - -/* backward recurrence */ - -pk = 1.0; -pkm1 = 1.0/ans; -k = n-1; -r = 2 * k; - -do - { - pkm2 = (pkm1 * r - pk * x) / x; - pk = pkm1; - pkm1 = pkm2; - r -= 2.0; - } -while( --k > 0 ); - -if( fabs(pk) > fabs(pkm1) ) - ans = j1(x)/pk; -else - ans = j0(x)/pkm1; -return( sign * ans ); -} |