1 files changed, 402 insertions, 0 deletions
diff --git a/libm/double/i1.c b/libm/double/i1.c
new file mode 100644
index 000000000..dfde216dc
--- /dev/null
+++ b/libm/double/i1.c
@@ -0,0 +1,402 @@
+/*							i1.c
+ *
+ *	Modified Bessel function of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i1();
+ *
+ * y = i1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns modified Bessel function of order one of the
+ * argument.
+ *
+ * The function is defined as i1(x) = -i j1( ix ).
+ *
+ * The range is partitioned into the two intervals [0,8] and
+ * (8, infinity).  Chebyshev polynomial expansions are employed
+ * in each interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    DEC       0, 30        3400       1.2e-16     2.3e-17
+ *    IEEE      0, 30       30000       1.9e-15     2.1e-16
+ *
+ *
+ */
+/*							i1e.c
+ *
+ *	Modified Bessel function of order one,
+ *	exponentially scaled
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, i1e();
+ *
+ * y = i1e( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns exponentially scaled modified Bessel function
+ * of order one of the argument.
+ *
+ * The function is defined as i1(x) = -i exp(-|x|) j1( ix ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0, 30       30000       2.0e-15     2.0e-16
+ * See i1().
+ *
+ */
+
+/*							i1.c 2		*/
+
+
+/*
+Cephes Math Library Release 2.8:  June, 2000
+Copyright 1985, 1987, 2000 by Stephen L. Moshier
+*/
+
+#include <math.h>
+
+/* Chebyshev coefficients for exp(-x) I1(x) / x
+ * in the interval [0,8].
+ *
+ * lim(x->0){ exp(-x) I1(x) / x } = 1/2.
+ */
+
+#ifdef UNK
+static double A[] =
+{
+ 2.77791411276104639959E-18,
+-2.11142121435816608115E-17,
+ 1.55363195773620046921E-16,
+-1.10559694773538630805E-15,
+ 7.60068429473540693410E-15,
+-5.04218550472791168711E-14,
+ 3.22379336594557470981E-13,
+-1.98397439776494371520E-12,
+ 1.17361862988909016308E-11,
+-6.66348972350202774223E-11,
+ 3.62559028155211703701E-10,
+-1.88724975172282928790E-9,
+ 9.38153738649577178388E-9,
+-4.44505912879632808065E-8,
+ 2.00329475355213526229E-7,
+-8.56872026469545474066E-7,
+ 3.47025130813767847674E-6,
+-1.32731636560394358279E-5,
+ 4.78156510755005422638E-5,
+-1.61760815825896745588E-4,
+ 5.12285956168575772895E-4,
+-1.51357245063125314899E-3,
+ 4.15642294431288815669E-3,
+-1.05640848946261981558E-2,
+ 2.47264490306265168283E-2,
+-5.29459812080949914269E-2,
+ 1.02643658689847095384E-1,
+-1.76416518357834055153E-1,
+ 2.52587186443633654823E-1
+};
+#endif
+
+#ifdef DEC
+static unsigned short A[] = {
+0021514,0174520,0060742,0000241,
+0122302,0137206,0016120,0025663,
+0023063,0017437,0026235,0176536,
+0123637,0052523,0170150,0125632,
+0024410,0165770,0030251,0044134,
+0125143,0012160,0162170,0054727,
+0025665,0075702,0035716,0145247,
+0126413,0116032,0176670,0015462,
+0027116,0073425,0110351,0105242,
+0127622,0104034,0137530,0037364,
+0030307,0050645,0120776,0175535,
+0131001,0130331,0043523,0037455,
+0031441,0026160,0010712,0100174,
+0132076,0164761,0022706,0017500,
+0032527,0015045,0115076,0104076,
+0133146,0001714,0015434,0144520,
+0033550,0161166,0124215,0077050,
+0134136,0127715,0143365,0157170,
+0034510,0106652,0013070,0064130,
+0135051,0117126,0117264,0123761,
+0035406,0045355,0133066,0175751,
+0135706,0061420,0054746,0122440,
+0036210,0031232,0047235,0006640,
+0136455,0012373,0144235,0011523,
+0036712,0107437,0036731,0015111,
+0137130,0156742,0115744,0172743,
+0037322,0033326,0124667,0124740,
+0137464,0123210,0021510,0144556,
+0037601,0051433,0111123,0177721
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short A[] = {
+0x4014,0x0c3c,0x9f2a,0x3c49,
+0x0576,0xc38a,0x57d0,0xbc78,
+0xbfac,0xe593,0x63e3,0x3ca6,
+0x1573,0x7e0d,0xeaaa,0xbcd3,
+0x290c,0x0615,0x1d7f,0x3d01,
+0x0b3b,0x1c8f,0x628e,0xbd2c,
+0xd955,0x4779,0xaf78,0x3d56,
+0x0366,0x5fb7,0x7383,0xbd81,
+0x3154,0xb21d,0xcee2,0x3da9,
+0x07de,0x97eb,0x5103,0xbdd2,
+0xdf6c,0xb43f,0xea34,0x3df8,
+0x67e6,0x28ea,0x361b,0xbe20,
+0x5010,0x0239,0x258e,0x3e44,
+0xc3e8,0x24b8,0xdd3e,0xbe67,
+0xd108,0xb347,0xe344,0x3e8a,
+0x992a,0x8363,0xc079,0xbeac,
+0xafc5,0xd511,0x1c4e,0x3ecd,
+0xbbcf,0xb8de,0xd5f9,0xbeeb,
+0x0d0b,0x42c7,0x11b5,0x3f09,
+0x94fe,0xd3d6,0x33ca,0xbf25,
+0xdf7d,0xb6c6,0xc95d,0x3f40,
+0xd4a4,0x0b3c,0xcc62,0xbf58,
+0xa1b4,0x49d3,0x0653,0x3f71,
+0xa26a,0x7913,0xa29f,0xbf85,
+0x2349,0xe7bb,0x51e3,0x3f99,
+0x9ebc,0x537c,0x1bbc,0xbfab,
+0xf53c,0xd536,0x46da,0x3fba,
+0x192e,0x0469,0x94d1,0xbfc6,
+0x7ffa,0x724a,0x2a63,0x3fd0
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short A[] = {
+0x3c49,0x9f2a,0x0c3c,0x4014,
+0xbc78,0x57d0,0xc38a,0x0576,
+0x3ca6,0x63e3,0xe593,0xbfac,
+0xbcd3,0xeaaa,0x7e0d,0x1573,
+0x3d01,0x1d7f,0x0615,0x290c,
+0xbd2c,0x628e,0x1c8f,0x0b3b,
+0x3d56,0xaf78,0x4779,0xd955,
+0xbd81,0x7383,0x5fb7,0x0366,
+0x3da9,0xcee2,0xb21d,0x3154,
+0xbdd2,0x5103,0x97eb,0x07de,
+0x3df8,0xea34,0xb43f,0xdf6c,
+0xbe20,0x361b,0x28ea,0x67e6,
+0x3e44,0x258e,0x0239,0x5010,
+0xbe67,0xdd3e,0x24b8,0xc3e8,
+0x3e8a,0xe344,0xb347,0xd108,
+0xbeac,0xc079,0x8363,0x992a,
+0x3ecd,0x1c4e,0xd511,0xafc5,
+0xbeeb,0xd5f9,0xb8de,0xbbcf,
+0x3f09,0x11b5,0x42c7,0x0d0b,
+0xbf25,0x33ca,0xd3d6,0x94fe,
+0x3f40,0xc95d,0xb6c6,0xdf7d,
+0xbf58,0xcc62,0x0b3c,0xd4a4,
+0x3f71,0x0653,0x49d3,0xa1b4,
+0xbf85,0xa29f,0x7913,0xa26a,
+0x3f99,0x51e3,0xe7bb,0x2349,
+0xbfab,0x1bbc,0x537c,0x9ebc,
+0x3fba,0x46da,0xd536,0xf53c,
+0xbfc6,0x94d1,0x0469,0x192e,
+0x3fd0,0x2a63,0x724a,0x7ffa
+};
+#endif
+
+/*							i1.c	*/
+
+/* Chebyshev coefficients for exp(-x) sqrt(x) I1(x)
+ * in the inverted interval [8,infinity].
+ *
+ * lim(x->inf){ exp(-x) sqrt(x) I1(x) } = 1/sqrt(2pi).
+ */
+
+#ifdef UNK
+static double B[] =
+{
+ 7.51729631084210481353E-18,
+ 4.41434832307170791151E-18,
+-4.65030536848935832153E-17,
+-3.20952592199342395980E-17,
+ 2.96262899764595013876E-16,
+ 3.30820231092092828324E-16,
+-1.88035477551078244854E-15,
+-3.81440307243700780478E-15,
+ 1.04202769841288027642E-14,
+ 4.27244001671195135429E-14,
+-2.10154184277266431302E-14,
+-4.08355111109219731823E-13,
+-7.19855177624590851209E-13,
+ 2.03562854414708950722E-12,
+ 1.41258074366137813316E-11,
+ 3.25260358301548823856E-11,
+-1.89749581235054123450E-11,
+-5.58974346219658380687E-10,
+-3.83538038596423702205E-9,
+-2.63146884688951950684E-8,
+-2.51223623787020892529E-7,
+-3.88256480887769039346E-6,
+-1.10588938762623716291E-4,
+-9.76109749136146840777E-3,
+ 7.78576235018280120474E-1
+};
+#endif
+
+#ifdef DEC
+static unsigned short B[] = {
+0022012,0125555,0115227,0043456,
+0021642,0156127,0052075,0145203,
+0122526,0072435,0111231,0011664,
+0122424,0001544,0161671,0114403,
+0023252,0144257,0163532,0142121,
+0023276,0132162,0174045,0013204,
+0124007,0077154,0057046,0110517,
+0124211,0066650,0116127,0157073,
+0024473,0133413,0130551,0107504,
+0025100,0064741,0032631,0040364,
+0124675,0045101,0071551,0012400,
+0125745,0161054,0071637,0011247,
+0126112,0117410,0035525,0122231,
+0026417,0037237,0131034,0176427,
+0027170,0100373,0024742,0025725,
+0027417,0006417,0105303,0141446,
+0127246,0163716,0121202,0060137,
+0130431,0123122,0120436,0166000,
+0131203,0144134,0153251,0124500,
+0131742,0005234,0122732,0033006,
+0132606,0157751,0072362,0121031,
+0133602,0043372,0047120,0015626,
+0134747,0165774,0001125,0046462,
+0136437,0166402,0117746,0155137,
+0040107,0050305,0125330,0124241
+};
+#endif
+
+#ifdef IBMPC
+static unsigned short B[] = {
+0xe8e6,0xb352,0x556d,0x3c61,
+0xb950,0xea87,0x5b8a,0x3c54,
+0x2277,0xb253,0xcea3,0xbc8a,
+0x3320,0x9c77,0x806c,0xbc82,
+0x588a,0xfceb,0x5915,0x3cb5,
+0xa2d1,0x5f04,0xd68e,0x3cb7,
+0xd22a,0x8bc4,0xefcd,0xbce0,
+0xfbc7,0x138a,0x2db5,0xbcf1,
+0x31e8,0x762d,0x76e1,0x3d07,
+0x281e,0x26b3,0x0d3c,0x3d28,
+0x22a0,0x2e6d,0xa948,0xbd17,
+0xe255,0x8e73,0xbc45,0xbd5c,
+0xb493,0x076a,0x53e1,0xbd69,
+0x9fa3,0xf643,0xe7d3,0x3d81,
+0x457b,0x653c,0x101f,0x3daf,
+0x7865,0xf158,0xe1a1,0x3dc1,
+0x4c0c,0xd450,0xdcf9,0xbdb4,
+0xdd80,0x5423,0x34ca,0xbe03,
+0x3528,0x9ad5,0x790b,0xbe30,
+0x46c1,0x94bb,0x4153,0xbe5c,
+0x5443,0x2e9e,0xdbfd,0xbe90,
+0x0373,0x49ca,0x48df,0xbed0,
+0xa9a6,0x804a,0xfd7f,0xbf1c,
+0xdb4c,0x53fc,0xfda0,0xbf83,
+0x1514,0xb55b,0xea18,0x3fe8
+};
+#endif
+
+#ifdef MIEEE
+static unsigned short B[] = {
+0x3c61,0x556d,0xb352,0xe8e6,
+0x3c54,0x5b8a,0xea87,0xb950,
+0xbc8a,0xcea3,0xb253,0x2277,
+0xbc82,0x806c,0x9c77,0x3320,
+0x3cb5,0x5915,0xfceb,0x588a,
+0x3cb7,0xd68e,0x5f04,0xa2d1,
+0xbce0,0xefcd,0x8bc4,0xd22a,
+0xbcf1,0x2db5,0x138a,0xfbc7,
+0x3d07,0x76e1,0x762d,0x31e8,
+0x3d28,0x0d3c,0x26b3,0x281e,
+0xbd17,0xa948,0x2e6d,0x22a0,
+0xbd5c,0xbc45,0x8e73,0xe255,
+0xbd69,0x53e1,0x076a,0xb493,
+0x3d81,0xe7d3,0xf643,0x9fa3,
+0x3daf,0x101f,0x653c,0x457b,
+0x3dc1,0xe1a1,0xf158,0x7865,
+0xbdb4,0xdcf9,0xd450,0x4c0c,
+0xbe03,0x34ca,0x5423,0xdd80,
+0xbe30,0x790b,0x9ad5,0x3528,
+0xbe5c,0x4153,0x94bb,0x46c1,
+0xbe90,0xdbfd,0x2e9e,0x5443,
+0xbed0,0x48df,0x49ca,0x0373,
+0xbf1c,0xfd7f,0x804a,0xa9a6,
+0xbf83,0xfda0,0x53fc,0xdb4c,
+0x3fe8,0xea18,0xb55b,0x1514
+};
+#endif
+
+/*							i1.c	*/
+#ifdef ANSIPROT
+extern double chbevl ( double, void *, int );
+extern double exp ( double );
+extern double sqrt ( double );
+extern double fabs ( double );
+#else
+double chbevl(), exp(), sqrt(), fabs();
+#endif
+
+double i1(x)
+double x;
+{ 
+double y, z;
+
+z = fabs(x);
+if( z <= 8.0 )
+	{
+	y = (z/2.0) - 2.0;
+	z = chbevl( y, A, 29 ) * z * exp(z);
+	}
+else
+	{
+	z = exp(z) * chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z);
+	}
+if( x < 0.0 )
+	z = -z;
+return( z );
+}
+
+/*							i1e()	*/
+
+double i1e( x )
+double x;
+{ 
+double y, z;
+
+z = fabs(x);
+if( z <= 8.0 )
+	{
+	y = (z/2.0) - 2.0;
+	z = chbevl( y, A, 29 ) * z;
+	}
+else
+	{
+	z = chbevl( 32.0/z - 2.0, B, 25 ) / sqrt(z);
+	}
+if( x < 0.0 )
+	z = -z;
+return( z );
+}