summaryrefslogtreecommitdiff
path: root/libm/double/jn.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/jn.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (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.c133
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 );
-}