summaryrefslogtreecommitdiff
path: root/libm/ldouble/exp2l.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/ldouble/exp2l.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/ldouble/exp2l.c')
-rw-r--r--libm/ldouble/exp2l.c166
1 files changed, 166 insertions, 0 deletions
diff --git a/libm/ldouble/exp2l.c b/libm/ldouble/exp2l.c
new file mode 100644
index 000000000..076f8bca5
--- /dev/null
+++ b/libm/ldouble/exp2l.c
@@ -0,0 +1,166 @@
+/* exp2l.c
+ *
+ * Base 2 exponential function, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, exp2l();
+ *
+ * y = exp2l( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns 2 raised to the x power.
+ *
+ * Range reduction is accomplished by separating the argument
+ * into an integer k and fraction f such that
+ * x k f
+ * 2 = 2 2.
+ *
+ * A Pade' form
+ *
+ * 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )
+ *
+ * approximates 2**x in the basic range [-0.5, 0.5].
+ *
+ *
+ * ACCURACY:
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE +-16300 300000 9.1e-20 2.6e-20
+ *
+ *
+ * See exp.c for comments on error amplification.
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * exp2l underflow x < -16382 0.0
+ * exp2l overflow x >= 16384 MAXNUM
+ *
+ */
+
+
+/*
+Cephes Math Library Release 2.7: May, 1998
+Copyright 1984, 1991, 1998 by Stephen L. Moshier
+*/
+
+
+
+#include <math.h>
+
+#ifdef UNK
+static long double P[] = {
+ 6.0614853552242266094567E1L,
+ 3.0286971917562792508623E4L,
+ 2.0803843631901852422887E6L,
+};
+static long double Q[] = {
+/* 1.0000000000000000000000E0,*/
+ 1.7492876999891839021063E3L,
+ 3.2772515434906797273099E5L,
+ 6.0027204078348487957118E6L,
+};
+#endif
+
+
+#ifdef IBMPC
+static short P[] = {
+0xffd8,0x6ad6,0x9c2b,0xf275,0x4004, XPD
+0x3426,0x2dc5,0xf19f,0xec9d,0x400d, XPD
+0x7ec0,0xd041,0x02e7,0xfdf4,0x4013, XPD
+};
+static short Q[] = {
+/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/
+0x575b,0x9b93,0x34d6,0xdaa9,0x4009, XPD
+0xe38d,0x6d74,0xa4f0,0xa005,0x4011, XPD
+0xb37e,0xcfba,0x40d0,0xb730,0x4015, XPD
+};
+#endif
+
+#ifdef MIEEE
+static long P[] = {
+0x40040000,0xf2759c2b,0x6ad6ffd8,
+0x400d0000,0xec9df19f,0x2dc53426,
+0x40130000,0xfdf402e7,0xd0417ec0,
+};
+static long Q[] = {
+/*0x3fff0000,0x80000000,0x00000000,*/
+0x40090000,0xdaa934d6,0x9b93575b,
+0x40110000,0xa005a4f0,0x6d74e38d,
+0x40150000,0xb73040d0,0xcfbab37e,
+};
+#endif
+
+#define MAXL2L 16384.0L
+#define MINL2L -16382.0L
+
+
+extern long double MAXNUML;
+#ifdef ANSIPROT
+extern long double polevll ( long double, void *, int );
+extern long double p1evll ( long double, void *, int );
+extern long double floorl ( long double );
+extern long double ldexpl ( long double, int );
+extern int isnanl ( long double );
+#else
+long double polevll(), p1evll(), floorl(), ldexpl(), isnanl();
+#endif
+#ifdef INFINITIES
+extern long double INFINITYL;
+#endif
+
+long double exp2l(x)
+long double x;
+{
+long double px, xx;
+int n;
+
+#ifdef NANS
+if( isnanl(x) )
+ return(x);
+#endif
+if( x > MAXL2L)
+ {
+#ifdef INFINITIES
+ return( INFINITYL );
+#else
+ mtherr( "exp2l", OVERFLOW );
+ return( MAXNUML );
+#endif
+ }
+
+if( x < MINL2L )
+ {
+#ifndef INFINITIES
+ mtherr( "exp2l", UNDERFLOW );
+#endif
+ return(0.0L);
+ }
+
+xx = x; /* save x */
+/* separate into integer and fractional parts */
+px = floorl(x+0.5L);
+n = px;
+x = x - px;
+
+/* rational approximation
+ * exp2(x) = 1.0 + 2xP(xx)/(Q(xx) - P(xx))
+ * where xx = x**2
+ */
+xx = x * x;
+px = x * polevll( xx, P, 2 );
+x = px / ( p1evll( xx, Q, 3 ) - px );
+x = 1.0L + ldexpl( x, 1 );
+
+/* scale by power of 2 */
+x = ldexpl( x, n );
+return(x);
+}