summaryrefslogtreecommitdiff
path: root/libm/float/j1f.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/float/j1f.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/j1f.c')
-rw-r--r--libm/float/j1f.c211
1 files changed, 211 insertions, 0 deletions
diff --git a/libm/float/j1f.c b/libm/float/j1f.c
new file mode 100644
index 000000000..4306e9747
--- /dev/null
+++ b/libm/float/j1f.c
@@ -0,0 +1,211 @@
+/* j1f.c
+ *
+ * Bessel function of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, j1f();
+ *
+ * y = j1f( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of order one of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval a polynomial approximation
+ * 2
+ * (w - r ) x P(w)
+ * 1
+ * 2
+ * is used, where w = x and r is the first zero of the function.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x R(1/x^2) - 3pi/4. The function is
+ *
+ * j0(x) = Modulus(x) cos( Phase(x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 1.2e-7 2.5e-8
+ * IEEE 2, 32 100000 2.0e-7 5.3e-8
+ *
+ *
+ */
+ /* y1.c
+ *
+ * Bessel function of second kind of order one
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, y1();
+ *
+ * y = y1( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns Bessel function of the second kind of order one
+ * of the argument.
+ *
+ * The domain is divided into the intervals [0, 2] and
+ * (2, infinity). In the first interval a rational approximation
+ * R(x) is employed to compute
+ *
+ * 2
+ * y0(x) = (w - r ) x R(x^2) + 2/pi (ln(x) j1(x) - 1/x) .
+ * 1
+ *
+ * Thus a call to j1() is required.
+ *
+ * In the second interval, the modulus and phase are approximated
+ * by polynomials of the form Modulus(x) = sqrt(1/x) Q(1/x)
+ * and Phase(x) = x + 1/x S(1/x^2) - 3pi/4. Then the function is
+ *
+ * y0(x) = Modulus(x) sin( Phase(x) ).
+ *
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Absolute error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0, 2 100000 2.2e-7 4.6e-8
+ * IEEE 2, 32 100000 1.9e-7 5.3e-8
+ *
+ * (error criterion relative when |y1| > 1).
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+#include <math.h>
+
+
+static float JP[5] = {
+-4.878788132172128E-009f,
+ 6.009061827883699E-007f,
+-4.541343896997497E-005f,
+ 1.937383947804541E-003f,
+-3.405537384615824E-002f
+};
+
+static float YP[5] = {
+ 8.061978323326852E-009f,
+-9.496460629917016E-007f,
+ 6.719543806674249E-005f,
+-2.641785726447862E-003f,
+ 4.202369946500099E-002f
+};
+
+static float MO1[8] = {
+ 6.913942741265801E-002f,
+-2.284801500053359E-001f,
+ 3.138238455499697E-001f,
+-2.102302420403875E-001f,
+ 5.435364690523026E-003f,
+ 1.493389585089498E-001f,
+ 4.976029650847191E-006f,
+ 7.978845453073848E-001f
+};
+
+static float PH1[8] = {
+-4.497014141919556E+001f,
+ 5.073465654089319E+001f,
+-2.485774108720340E+001f,
+ 7.222973196770240E+000f,
+-1.544842782180211E+000f,
+ 3.503787691653334E-001f,
+-1.637986776941202E-001f,
+ 3.749989509080821E-001f
+};
+
+static float YO1 = 4.66539330185668857532f;
+static float Z1 = 1.46819706421238932572E1f;
+
+static float THPIO4F = 2.35619449019234492885f; /* 3*pi/4 */
+static float TWOOPI = 0.636619772367581343075535f; /* 2/pi */
+extern float PIO4;
+
+
+float polevlf(float, float *, int);
+float logf(float), sinf(float), cosf(float), sqrtf(float);
+
+float j1f( float xx )
+{
+float x, w, z, p, q, xn;
+
+
+x = xx;
+if( x < 0 )
+ x = -xx;
+
+if( x <= 2.0f )
+ {
+ z = x * x;
+ p = (z-Z1) * x * polevlf( z, JP, 4 );
+ return( p );
+ }
+
+q = 1.0f/x;
+w = sqrtf(q);
+
+p = w * polevlf( q, MO1, 7);
+w = q*q;
+xn = q * polevlf( w, PH1, 7) - THPIO4F;
+p = p * cosf(xn + x);
+return(p);
+}
+
+
+
+
+extern float MAXNUMF;
+
+float y1f( float xx )
+{
+float x, w, z, p, q, xn;
+
+
+x = xx;
+if( x <= 2.0f )
+ {
+ if( x <= 0.0f )
+ {
+ mtherr( "y1f", DOMAIN );
+ return( -MAXNUMF );
+ }
+ z = x * x;
+ w = (z - YO1) * x * polevlf( z, YP, 4 );
+ w += TWOOPI * ( j1f(x) * logf(x) - 1.0f/x );
+ return( w );
+ }
+
+q = 1.0f/x;
+w = sqrtf(q);
+
+p = w * polevlf( q, MO1, 7);
+w = q*q;
+xn = q * polevlf( w, PH1, 7) - THPIO4F;
+p = p * sinf(xn + x);
+return(p);
+}