1 files changed, 333 insertions, 0 deletions
diff --git a/libm/double/k0.c b/libm/double/k0.c
new file mode 100644
index 000000000..7d09cb4a1
--- /dev/null
+++ b/libm/double/k0.c
@@ -0,0 +1,333 @@
+/*							k0.c
+ *
+ *	Modified Bessel function, third kind, order zero
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, k0();
+ *
+ * y = k0( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns modified Bessel function of the third kind
+ * of order zero of the argument.
+ *
+ * The range is partitioned into the two intervals [0,8] and
+ * (8, infinity).  Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at 2000 random points between 0 and 8.  Peak absolute
+ * error (relative when K0 > 1) was 1.46e-14; rms, 4.26e-15.
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 30        3100       1.3e-16     2.1e-17
+ *    IEEE      0, 30       30000       1.2e-15     1.6e-16
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ *  K0 domain          x <= 0          MAXNUM
+ *
+ */
+/*							k0e()
+ *
+ *	Modified Bessel function, third kind, order zero,
+ *	exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, k0e();
+ *
+ * y = k0e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of the third kind of order zero of the argument.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0, 30       30000       1.4e-15     1.4e-16
+ * See k0().
+ *
+ */
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1984, 1987, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+/* Chebyshev coefficients for K0(x) + log(x/2) I0(x)
+ * in the interval [0,2].  The odd order coefficients are all
+ * zero; only the even order coefficients are listed.
+ * 
+ * lim(x->0){ K0(x) + log(x/2) I0(x) } = -EUL.
+ */
+
+#ifdef UNK
+static double A[] =
+{
+ 1.37446543561352307156E-16,
+ 4.25981614279661018399E-14,
+ 1.03496952576338420167E-11,
+ 1.90451637722020886025E-9,
+ 2.53479107902614945675E-7,
+ 2.28621210311945178607E-5,
+ 1.26461541144692592338E-3,
+ 3.59799365153615016266E-2,
+ 3.44289899924628486886E-1,
+-5.35327393233902768720E-1
+};
+#endif
+
+#ifdef DEC
+static unsigned short A[] = {
+0023036,0073417,0032477,0165673,
+0025077,0154126,0016046,0012517,
+0027066,0011342,0035211,0005041,
+0031002,0160233,0037454,0050224,
+0032610,0012747,0037712,0173741,
+0034277,0144007,0172147,0162375,
+0035645,0140563,0125431,0165626,
+0037023,0057662,0125124,0102051,
+0037660,0043304,0004411,0166707,
+0140011,0005467,0047227,0130370
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short A[] = {
+0xfd77,0xe6a7,0xcee1,0x3ca3,
+0xc2aa,0xc384,0xfb0a,0x3d27,
+0x2144,0x4751,0xc25c,0x3da6,
+0x8a13,0x67e5,0x5c13,0x3e20,
+0x5efc,0xe7f9,0x02bc,0x3e91,
+0xfca0,0xfe8c,0xf900,0x3ef7,
+0x3d73,0x7563,0xb82e,0x3f54,
+0x9085,0x554a,0x6bf6,0x3fa2,
+0x3db9,0x8121,0x08d8,0x3fd6,
+0xf61f,0xe9d2,0x2166,0xbfe1
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short A[] = {
+0x3ca3,0xcee1,0xe6a7,0xfd77,
+0x3d27,0xfb0a,0xc384,0xc2aa,
+0x3da6,0xc25c,0x4751,0x2144,
+0x3e20,0x5c13,0x67e5,0x8a13,
+0x3e91,0x02bc,0xe7f9,0x5efc,
+0x3ef7,0xf900,0xfe8c,0xfca0,
+0x3f54,0xb82e,0x7563,0x3d73,
+0x3fa2,0x6bf6,0x554a,0x9085,
+0x3fd6,0x08d8,0x8121,0x3db9,
+0xbfe1,0x2166,0xe9d2,0xf61f
+};
+#endif
+
+
+
+/* Chebyshev coefficients for exp(x) sqrt(x) K0(x)
+ * in the inverted interval [2,infinity].
+ * 
+ * lim(x->inf){ exp(x) sqrt(x) K0(x) } = sqrt(pi/2).
+ */
+
+#ifdef UNK
+static double B[] = {
+ 5.30043377268626276149E-18,
+-1.64758043015242134646E-17,
+ 5.21039150503902756861E-17,
+-1.67823109680541210385E-16,
+ 5.51205597852431940784E-16,
+-1.84859337734377901440E-15,
+ 6.34007647740507060557E-15,
+-2.22751332699166985548E-14,
+ 8.03289077536357521100E-14,
+-2.98009692317273043925E-13,
+ 1.14034058820847496303E-12,
+-4.51459788337394416547E-12,
+ 1.85594911495471785253E-11,
+-7.95748924447710747776E-11,
+ 3.57739728140030116597E-10,
+-1.69753450938905987466E-9,
+ 8.57403401741422608519E-9,
+-4.66048989768794782956E-8,
+ 2.76681363944501510342E-7,
+-1.83175552271911948767E-6,
+ 1.39498137188764993662E-5,
+-1.28495495816278026384E-4,
+ 1.56988388573005337491E-3,
+-3.14481013119645005427E-2,
+ 2.44030308206595545468E0
+};
+#endif
+
+#ifdef DEC
+static unsigned short B[] = {
+0021703,0106456,0076144,0173406,
+0122227,0173144,0116011,0030033,
+0022560,0044562,0006506,0067642,
+0123101,0076243,0123273,0131013,
+0023436,0157713,0056243,0141331,
+0124005,0032207,0063726,0164664,
+0024344,0066342,0051756,0162300,
+0124710,0121365,0154053,0077022,
+0025264,0161166,0066246,0077420,
+0125647,0141671,0006443,0103212,
+0026240,0076431,0077147,0160445,
+0126636,0153741,0174002,0105031,
+0027243,0040102,0035375,0163073,
+0127656,0176256,0113476,0044653,
+0030304,0125544,0006377,0130104,
+0130751,0047257,0110537,0127324,
+0031423,0046400,0014772,0012164,
+0132110,0025240,0155247,0112570,
+0032624,0105314,0007437,0021574,
+0133365,0155243,0174306,0116506,
+0034152,0004776,0061643,0102504,
+0135006,0136277,0036104,0175023,
+0035715,0142217,0162474,0115022,
+0137000,0147671,0065177,0134356,
+0040434,0026754,0175163,0044070
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short B[] = {
+0x9ee1,0xcf8c,0x71a5,0x3c58,
+0x2603,0x9381,0xfecc,0xbc72,
+0xcdf4,0x41a8,0x092e,0x3c8e,
+0x7641,0x74d7,0x2f94,0xbca8,
+0x785b,0x6b94,0xdbf9,0x3cc3,
+0xdd36,0xecfa,0xa690,0xbce0,
+0xdc98,0x4a7d,0x8d9c,0x3cfc,
+0x6fc2,0xbb05,0x145e,0xbd19,
+0xcfe2,0xcd94,0x9c4e,0x3d36,
+0x70d1,0x21a4,0xf877,0xbd54,
+0xfc25,0x2fcc,0x0fa3,0x3d74,
+0x5143,0x3f00,0xdafc,0xbd93,
+0xbcc7,0x475f,0x6808,0x3db4,
+0xc935,0xd2e7,0xdf95,0xbdd5,
+0xf608,0x819f,0x956c,0x3df8,
+0xf5db,0xf22b,0x29d5,0xbe1d,
+0x428e,0x033f,0x69a0,0x3e42,
+0xf2af,0x1b54,0x0554,0xbe69,
+0xe46f,0x81e3,0x9159,0x3e92,
+0xd3a9,0x7f18,0xbb54,0xbebe,
+0x70a9,0xcc74,0x413f,0x3eed,
+0x9f42,0xe788,0xd797,0xbf20,
+0x9342,0xfca7,0xb891,0x3f59,
+0xf71e,0x2d4f,0x19f7,0xbfa0,
+0x6907,0x9f4e,0x85bd,0x4003
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short B[] = {
+0x3c58,0x71a5,0xcf8c,0x9ee1,
+0xbc72,0xfecc,0x9381,0x2603,
+0x3c8e,0x092e,0x41a8,0xcdf4,
+0xbca8,0x2f94,0x74d7,0x7641,
+0x3cc3,0xdbf9,0x6b94,0x785b,
+0xbce0,0xa690,0xecfa,0xdd36,
+0x3cfc,0x8d9c,0x4a7d,0xdc98,
+0xbd19,0x145e,0xbb05,0x6fc2,
+0x3d36,0x9c4e,0xcd94,0xcfe2,
+0xbd54,0xf877,0x21a4,0x70d1,
+0x3d74,0x0fa3,0x2fcc,0xfc25,
+0xbd93,0xdafc,0x3f00,0x5143,
+0x3db4,0x6808,0x475f,0xbcc7,
+0xbdd5,0xdf95,0xd2e7,0xc935,
+0x3df8,0x956c,0x819f,0xf608,
+0xbe1d,0x29d5,0xf22b,0xf5db,
+0x3e42,0x69a0,0x033f,0x428e,
+0xbe69,0x0554,0x1b54,0xf2af,
+0x3e92,0x9159,0x81e3,0xe46f,
+0xbebe,0xbb54,0x7f18,0xd3a9,
+0x3eed,0x413f,0xcc74,0x70a9,
+0xbf20,0xd797,0xe788,0x9f42,
+0x3f59,0xb891,0xfca7,0x9342,
+0xbfa0,0x19f7,0x2d4f,0xf71e,
+0x4003,0x85bd,0x9f4e,0x6907
+};
+#endif
+
+/*							k0.c	*/
+#ifdef ANSIPROT 
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double i0 ( double );
+extern double log ( double );
+extern double sqrt ( double );
+#else
+double chbevl(), exp(), i0(), log(), sqrt();
+#endif
+extern double PI;
+extern double MAXNUM;
+
+double k0(x)
+double x;
+{
+double y, z;
+
+if( x <= 0.0 )
+	{
+	mtherr( "k0", DOMAIN );
+	return( MAXNUM );
+	}
+
+if( x <= 2.0 )
+	{
+	y = x * x - 2.0;
+	y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x);
+	return( y );
+	}
+z = 8.0/x - 2.0;
+y = exp(-x) * chbevl( z, B, 25 ) / sqrt(x);
+return(y);
+}
+
+
+
+
+double k0e( x )
+double x;
+{
+double y;
+
+if( x <= 0.0 )
+	{
+	mtherr( "k0e", DOMAIN );
+	return( MAXNUM );
+	}
+
+if( x <= 2.0 )
+	{
+	y = x * x - 2.0;
+	y = chbevl( y, A, 10 ) - log( 0.5 * x ) * i0(x);
+	return( y * exp(x) );
+	}
+
+y = chbevl( 8.0/x - 2.0, B, 25 ) / sqrt(x);
+return(y);
+}