1 files changed, 174 insertions, 0 deletions
diff --git a/libm/float/k1f.c b/libm/float/k1f.c
new file mode 100644
index 000000000..d5b9bdfce
--- /dev/null
+++ b/libm/float/k1f.c
@@ -0,0 +1,174 @@
+/*							k1f.c
+ *
+ *	Modified Bessel function, third kind, order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k1f();
+ *
+ * y = k1f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Computes the modified Bessel function of the third kind
+ * of order one of the argument.
+ *
+ * The range is partitioned into the two intervals [0,2] and
+ * (2, infinity).  Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0, 30       30000       4.6e-7      7.6e-8
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * k1 domain          x <= 0          MAXNUM
+ *
+ */
+/*							k1ef.c
+ *
+ *	Modified Bessel function, third kind, order one,
+ *	exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, k1ef();
+ *
+ * y = k1ef( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order one of the argument:
+ *
+ *      k1e(x) = exp(x) * k1(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0, 30       30000       4.9e-7      6.7e-8
+ * See k1().
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+#include <math.h>
+
+/* Chebyshev coefficients for x(K1(x) - log(x/2) I1(x))
+ * in the interval [0,2].
+ * 
+ * lim(x->0){ x(K1(x) - log(x/2) I1(x)) } = 1.
+ */
+
+#define MINNUMF 6.0e-39
+static float A[] =
+{
+-2.21338763073472585583E-8f,
+-2.43340614156596823496E-6f,
+-1.73028895751305206302E-4f,
+-6.97572385963986435018E-3f,
+-1.22611180822657148235E-1f,
+-3.53155960776544875667E-1f,
+ 1.52530022733894777053E0f
+};
+
+
+
+
+/* Chebyshev coefficients for exp(x) sqrt(x) K1(x)
+ * in the interval [2,infinity].
+ *
+ * lim(x->inf){ exp(x) sqrt(x) K1(x) } = sqrt(pi/2).
+ */
+
+static float B[] =
+{
+ 2.01504975519703286596E-9f,
+-1.03457624656780970260E-8f,
+ 5.74108412545004946722E-8f,
+-3.50196060308781257119E-7f,
+ 2.40648494783721712015E-6f,
+-1.93619797416608296024E-5f,
+ 1.95215518471351631108E-4f,
+-2.85781685962277938680E-3f,
+ 1.03923736576817238437E-1f,
+ 2.72062619048444266945E0f
+};
+
+
+ 
+extern float MAXNUMF;
+#ifdef ANSIC
+float chbevlf(float, float *, int);
+float expf(float), i1f(float), logf(float), sqrtf(float);
+#else
+float chbevlf(), expf(), i1f(), logf(), sqrtf();
+#endif
+
+float k1f(float xx)
+{
+float x, y;
+
+x = xx;
+if( x <= MINNUMF )
+	{
+	mtherr( "k1f", DOMAIN );
+	return( MAXNUMF );
+	}
+
+if( x <= 2.0f )
+	{
+	y = x * x - 2.0f;
+	y =  logf( 0.5f * x ) * i1f(x)  +  chbevlf( y, A, 7 ) / x;
+	return( y );
+	}
+
+return(  expf(-x) * chbevlf( 8.0f/x - 2.0f, B, 10 ) / sqrtf(x) );
+
+}
+
+
+
+float k1ef( float xx )
+{
+float x, y;
+
+x = xx;
+if( x <= 0.0f )
+	{
+	mtherr( "k1ef", DOMAIN );
+	return( MAXNUMF );
+	}
+
+if( x <= 2.0f )
+	{
+	y = x * x - 2.0f;
+	y =  logf( 0.5f * x ) * i1f(x)  +  chbevlf( y, A, 7 ) / x;
+	return( y * expf(x) );
+	}
+
+return(  chbevlf( 8.0f/x - 2.0f, B, 10 ) / sqrtf(x) );
+
+}