1 files changed, 114 insertions, 0 deletions
diff --git a/libm/double/yn.c b/libm/double/yn.c
new file mode 100644
index 000000000..0c569a925
--- /dev/null
+++ b/libm/double/yn.c
@@ -0,0 +1,114 @@
+/*							yn.c
+ *
+ *	Bessel function of second kind of integer order
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, yn();
+ * int n;
+ *
+ * y = yn( n, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order n, where n is a
+ * (possibly negative) integer.
+ *
+ * The function is evaluated by forward recurrence on
+ * n, starting with values computed by the routines
+ * y0() and y1().
+ *
+ * If n = 0 or 1 the routine for y0 or y1 is called
+ * directly.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ *                      Absolute error, except relative
+ *                      when y > 1:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 30        2200       2.9e-16     5.3e-17
+ *    IEEE      0, 30       30000       3.4e-15     4.3e-16
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * yn singularity   x = 0              MAXNUM
+ * yn overflow                         MAXNUM
+ *
+ * Spot checked against tables for x, n between 0 and 100.
+ *
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+#ifdef ANSIPROT
+extern double y0 ( double );
+extern double y1 ( double );
+extern double log ( double );
+#else
+double y0(), y1(), log();
+#endif
+extern double MAXNUM, MAXLOG;
+
+double yn( n, x )
+int n;
+double x;
+{
+double an, anm1, anm2, r;
+int k, sign;
+
+if( n < 0 )
+	{
+	n = -n;
+	if( (n & 1) == 0 )	/* -1**n */
+		sign = 1;
+	else
+		sign = -1;
+	}
+else
+	sign = 1;
+
+
+if( n == 0 )
+	return( sign * y0(x) );
+if( n == 1 )
+	return( sign * y1(x) );
+
+/* test for overflow */
+if( x <= 0.0 )
+	{
+	mtherr( "yn", SING );
+	return( -MAXNUM );
+	}
+
+/* forward recurrence on n */
+
+anm2 = y0(x);
+anm1 = y1(x);
+k = 1;
+r = 2 * k;
+do
+	{
+	an = r * anm1 / x  -  anm2;
+	anm2 = anm1;
+	anm1 = an;
+	r += 2.0;
+	++k;
+	}
+while( k < n );
+
+
+return( sign * an );
+}