uClibc now has a math library. muahahahaha!

 -Erik

69 files changed, 27634 insertions, 0 deletions

diff --git a/libm/ldouble/Makefile b/libm/ldouble/Makefile
new file mode 100644
index 000000000..43395a140
--- /dev/null
+++ b/libm/ldouble/Makefile
@@ -0,0 +1,123 @@
+# Makefile for uClibc's math library
+#
+# Copyright (C) 2001 by Lineo, inc.
+#
+# This program is free software; you can redistribute it and/or modify it under
+# the terms of the GNU Library General Public License as published by the Free
+# Software Foundation; either version 2 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Library General Public License
+# along with this program; if not, write to the Free Software Foundation, Inc.,
+# 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# Derived in part from the Linux-8086 C library, the GNU C Library, and several
+# other sundry sources.  Files within this library are copyright by their
+# respective copyright holders.
+
+TOPDIR=../../
+include $(TOPDIR)Rules.mak
+
+LIBM=../libm.a
+TARGET_CC= $(TOPDIR)/extra/gcc-uClibc/$(TARGET_ARCH)-uclibc-gcc
+
+CSRC=acoshl.c asinhl.c asinl.c atanhl.c atanl.c bdtrl.c btdtrl.c cbrtl.c \
+	chdtrl.c coshl.c ellpel.c ellpkl.c elliel.c ellikl.c ellpjl.c \
+	exp10l.c exp2l.c expl.c fdtrl.c gammal.c gdtrl.c igamil.c igaml.c \
+	incbetl.c incbil.c isnanl.c j0l.c j1l.c jnl.c ldrand.c log10l.c log2l.c \
+	logl.c nbdtrl.c ndtril.c ndtrl.c pdtrl.c powl.c powil.c sinhl.c sinl.c \
+	sqrtl.c stdtrl.c tanhl.c tanl.c unityl.c ynl.c \
+	floorl.c polevll.c mtherr.c #cmplxl.c clogl.c
+COBJS=$(patsubst %.c,%.o, $(CSRC))
+
+
+OBJS=$(COBJS)
+
+all: $(OBJS) $(LIBM)
+
+$(LIBM): ar-target
+
+ar-target: $(OBJS)
+	$(AR) $(ARFLAGS) $(LIBM) $(OBJS)
+
+$(COBJS): %.o : %.c
+	$(TARGET_CC) $(CFLAGS) -c $< -o $@
+	$(STRIPTOOL) -x -R .note -R .comment $*.o
+
+$(OBJ): Makefile
+
+clean:
+	rm -f *.[oa] *~ core
+
+
+
+#-----------------------------------------
+
+
+#all: mtstl lparanoi lcalc fltestl nantst testvect monotl libml.a
+
+mtstl: libml.a mtstl.o $(OBJS)
+	$(CC) $(CFLAGS) -o mtstl mtstl.o libml.a $(LIBS)
+
+mtstl.o: mtstl.c
+
+lparanoi: libml.a lparanoi.o setprec.o ieee.o econst.o $(OBJS)
+	$(CC) $(CFLAGS) -o lparanoi lparanoi.o setprec.o ieee.o econst.o libml.a $(LIBS)
+
+lparanoi.o: lparanoi.c
+	$(CC) $(CFLAGS) -Wno-implicit -c lparanoi.c
+
+econst.o: econst.c ehead.h
+
+lcalc: libml.a lcalc.o ieee.o econst.o $(OBJS)
+	$(CC) $(CFLAGS) -o lcalc lcalc.o ieee.o econst.o libml.a $(LIBS)
+
+lcalc.o: lcalc.c lcalc.h ehead.h
+
+ieee.o: ieee.c ehead.h
+
+# Use $(OBJS) in ar command for libml.a if possible; else *.o
+libml.a: $(OBJS) mconf.h
+	ar -rv libml.a $(OBJS)
+	ranlib libml.a
+
+
+fltestl: fltestl.c libml.a
+	$(CC) $(CFLAGS) -o fltestl fltestl.c libml.a
+
+fltestl.o: fltestl.c
+
+flrtstl: flrtstl.c libml.a
+	$(CC) $(CFLAGS) -o flrtstl flrtstl.c libml.a
+
+flrtstl.o: flrtstl.c
+
+nantst: nantst.c libml.a
+	$(CC) $(CFLAGS) -o nantst nantst.c libml.a
+
+nantst.o: nantst.c
+
+testvect: testvect.o libml.a
+	$(CC) $(CFLAGS) -o testvect testvect.o libml.a
+
+testvect.o: testvect.c
+	$(CC) -g -c -o testvect.o testvect.c
+
+monotl: monotl.o libml.a
+	$(CC) $(CFLAGS) -o monotl monotl.o libml.a
+
+monotl.o: monotl.c
+	$(CC) -g -c -o monotl.o monotl.c
+
+# Run test programs
+check: mtstl fltestl testvect monotl libml.a
+	-mtstl
+	-fltestl
+	-testvect
+	-monotl
+
diff --git a/libm/ldouble/README.txt b/libm/ldouble/README.txt
new file mode 100644
index 000000000..30fcaad36
--- /dev/null
+++ b/libm/ldouble/README.txt
@@ -0,0 +1,3502 @@
+/*							acoshl.c
+ *
+ *	Inverse hyperbolic cosine, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, acoshl();
+ *
+ * y = acoshl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns inverse hyperbolic cosine of argument.
+ *
+ * If 1 <= x < 1.5, a rational approximation
+ *
+ *	sqrt(2z) * P(z)/Q(z)
+ *
+ * where z = x-1, is used.  Otherwise,
+ *
+ * acosh(x)  =  log( x + sqrt( (x-1)(x+1) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      1,3         30000       2.0e-19     3.9e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * acoshl domain      |x| < 1            0.0
+ *
+ */
+
+/*							asinhl.c
+ *
+ *	Inverse hyperbolic sine, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, asinhl();
+ *
+ * y = asinhl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns inverse hyperbolic sine of argument.
+ *
+ * If |x| < 0.5, the function is approximated by a rational
+ * form  x + x**3 P(x)/Q(x).  Otherwise,
+ *
+ *     asinh(x) = log( x + sqrt(1 + x*x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -3,3         30000       1.7e-19     3.5e-20
+ *
+ */
+
+/*							asinl.c
+ *
+ *	Inverse circular sine, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, asinl();
+ *
+ * y = asinl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
+ *
+ * A rational function of the form x + x**3 P(x**2)/Q(x**2)
+ * is used for |x| in the interval [0, 0.5].  If |x| > 0.5 it is
+ * transformed by the identity
+ *
+ *    asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     -1, 1        30000       2.7e-19     4.8e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * asin domain        |x| > 1           0.0
+ *
+ */
+/*							acosl()
+ *
+ *	Inverse circular cosine, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * double x, y, acosl();
+ *
+ * y = acosl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose cosine
+ * is x.
+ *
+ * Analytically, acos(x) = pi/2 - asin(x).  However if |x| is
+ * near 1, there is cancellation error in subtracting asin(x)
+ * from pi/2.  Hence if x < -0.5,
+ *
+ *    acos(x) =	 pi - 2.0 * asin( sqrt((1+x)/2) );
+ *
+ * or if x > +0.5,
+ *
+ *    acos(x) =	 2.0 * asin(  sqrt((1-x)/2) ).
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -1, 1       30000       1.4e-19     3.5e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * asin domain        |x| > 1           0.0
+ */
+
+/*							atanhl.c
+ *
+ *	Inverse hyperbolic tangent, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, atanhl();
+ *
+ * y = atanhl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns inverse hyperbolic tangent of argument in the range
+ * MINLOGL to MAXLOGL.
+ *
+ * If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is
+ * employed.  Otherwise,
+ *        atanh(x) = 0.5 * log( (1+x)/(1-x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -1,1        30000       1.1e-19     3.3e-20
+ *
+ */
+
+/*							atanl.c
+ *
+ *	Inverse circular tangent, long double precision
+ *      (arctangent)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, atanl();
+ *
+ * y = atanl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose tangent
+ * is x.
+ *
+ * Range reduction is from four intervals into the interval
+ * from zero to  tan( pi/8 ).  The approximant uses a rational
+ * function of degree 3/4 of the form x + x**3 P(x)/Q(x).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10, 10    150000       1.3e-19     3.0e-20
+ *
+ */
+/*							atan2l()
+ *
+ *	Quadrant correct inverse circular tangent,
+ *	long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, z, atan2l();
+ *
+ * z = atan2l( y, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle whose tangent is y/x.
+ * Define compile time symbol ANSIC = 1 for ANSI standard,
+ * range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
+ * 0 to 2PI, args (x,y).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10, 10     60000       1.7e-19     3.2e-20
+ * See atan.c.
+ *
+ */
+
+/*							bdtrl.c
+ *
+ *	Binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtrl();
+ *
+ * y = bdtrl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms 0 through k of the Binomial
+ * probability density:
+ *
+ *   k
+ *   --  ( n )   j      n-j
+ *   >   (   )  p  (1-p)
+ *   --  ( j )
+ *  j=0
+ *
+ * The terms are not summed directly; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * Tested at random points (k,n,p) with a and b between 0
+ * and 10000 and p between 0 and 1.
+ *    Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,10000      3000       1.6e-14     2.2e-15
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtrl domain        k < 0            0.0
+ *                     n < k
+ *                     x < 0, x > 1
+ *
+ */
+/*							bdtrcl()
+ *
+ *	Complemented binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtrcl();
+ *
+ * y = bdtrcl( k, n, p );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the sum of the terms k+1 through n of the Binomial
+ * probability density:
+ *
+ *   n
+ *   --  ( n )   j      n-j
+ *   >   (   )  p  (1-p)
+ *   --  ( j )
+ *  j=k+1
+ *
+ * The terms are not summed directly; instead the incomplete
+ * beta integral is employed, according to the formula
+ *
+ * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
+ *
+ * The arguments must be positive, with p ranging from 0 to 1.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * See incbet.c.
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtrcl domain     x<0, x>1, n<k       0.0
+ */
+/*							bdtril()
+ *
+ *	Inverse binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * long double p, y, bdtril();
+ *
+ * p = bdtril( k, n, y );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Finds the event probability p such that the sum of the
+ * terms 0 through k of the Binomial probability density
+ * is equal to the given cumulative probability y.
+ *
+ * This is accomplished using the inverse beta integral
+ * function and the relation
+ *
+ * 1 - p = incbi( n-k, k+1, y ).
+ *
+ * ACCURACY:
+ *
+ * See incbi.c.
+ * Tested at random k, n between 1 and 10000.  The "domain" refers to p:
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,1        3500       2.0e-15     8.2e-17
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtril domain     k < 0, n <= k         0.0
+ *                  x < 0, x > 1
+ */
+
+
+/*							btdtrl.c
+ *
+ *	Beta distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double a, b, x, y, btdtrl();
+ *
+ * y = btdtrl( a, b, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the area from zero to x under the beta density
+ * function:
+ *
+ *
+ *                          x
+ *            -             -
+ *           | (a+b)       | |  a-1      b-1
+ * P(x)  =  ----------     |   t    (1-t)    dt
+ *           -     -     | |
+ *          | (a) | (b)   -
+ *                         0
+ *
+ *
+ * The mean value of this distribution is a/(a+b).  The variance
+ * is ab/[(a+b)^2 (a+b+1)].
+ *
+ * This function is identical to the incomplete beta integral
+ * function, incbetl(a, b, x).
+ *
+ * The complemented function is
+ *
+ * 1 - P(1-x)  =  incbetl( b, a, x );
+ *
+ *
+ * ACCURACY:
+ *
+ * See incbetl.c.
+ *
+ */
+
+/*							cbrtl.c
+ *
+ *	Cube root, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, cbrtl();
+ *
+ * y = cbrtl( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the cube root of the argument, which may be negative.
+ *
+ * Range reduction involves determining the power of 2 of
+ * the argument.  A polynomial of degree 2 applied to the
+ * mantissa, and multiplication by the cube root of 1, 2, or 4
+ * approximates the root to within about 0.1%.  Then Newton's
+ * iteration is used three times to converge to an accurate
+ * result.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE     .125,8        80000      7.0e-20     2.2e-20
+ *    IEEE    exp(+-707)    100000      7.0e-20     2.4e-20
+ *
+ */
+
+/*							chdtrl.c
+ *
+ *	Chi-square distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double df, x, y, chdtrl();
+ *
+ * y = chdtrl( df, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the area under the left hand tail (from 0 to x)
+ * of the Chi square probability density function with
+ * v degrees of freedom.
+ *
+ *
+ *                                  inf.
+ *                                    -
+ *                        1          | |  v/2-1  -t/2
+ *  P( x | v )   =   -----------     |   t      e     dt
+ *                    v/2  -       | |
+ *                   2    | (v/2)   -
+ *                                   x
+ *
+ * where x is the Chi-square variable.
+ *
+ * The incomplete gamma integral is used, according to the
+ * formula
+ *
+ *	y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
+ *
+ *
+ * The arguments must both be positive.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ * See igam().
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * chdtr domain   x < 0 or v < 1        0.0
+ */
+/*							chdtrcl()
+ *
+ *	Complemented Chi-square distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double v, x, y, chdtrcl();
+ *
+ * y = chdtrcl( v, x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the area under the right hand tail (from x to
+ * infinity) of the Chi square probability dens