summaryrefslogtreecommitdiff
path: root/libm/double/jn.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/jn.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/jn.c')
-rw-r--r--libm/double/jn.c133
1 files changed, 133 insertions, 0 deletions
diff --git a/libm/double/jn.c b/libm/double/jn.c
new file mode 100644
index 000000000..ee05395aa
--- /dev/null
+++ b/libm/double/jn.c
@@ -0,0 +1,133 @@
+/* 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 );
+}