Totally rework the math library, this time based on the MacOs X

math library (which is itself based on the math lib from FreeBSD).
 -Erik

69 files changed, 0 insertions, 27633 deletions

diff --git a/libm/ldouble/Makefile b/libm/ldouble/Makefile
deleted file mode 100644
index dad448840..000000000
--- a/libm/ldouble/Makefile
+++ /dev/null
@@ -1,122 +0,0 @@
-# Makefile for uClibc's math library
-# Copyright (C) 2001 by Lineo, inc.
-#
-# This math library is derived primarily from the Cephes Math Library,
-# copyright by Stephen L. Moshier <moshier@world.std.com>
-#
-# 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
-#
-
-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) $(TARGET_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)
-	$(TARGET_CC) $(TARGET_CFLAGS) -o mtstl mtstl.o libml.a $(LIBS)
-
-mtstl.o: mtstl.c
-
-lparanoi: libml.a lparanoi.o setprec.o ieee.o econst.o $(OBJS)
-	$(TARGET_CC) $(TARGET_CFLAGS) -o lparanoi lparanoi.o setprec.o ieee.o econst.o libml.a $(LIBS)
-
-lparanoi.o: lparanoi.c
-	$(TARGET_CC) $(TARGET_CFLAGS) -Wno-implicit -c lparanoi.c
-
-econst.o: econst.c ehead.h
-
-lcalc: libml.a lcalc.o ieee.o econst.o $(OBJS)
-	$(TARGET_CC) $(TARGET_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
-	$(TARGET_CC) $(TARGET_CFLAGS) -o fltestl fltestl.c libml.a
-
-fltestl.o: fltestl.c
-
-flrtstl: flrtstl.c libml.a
-	$(TARGET_CC) $(TARGET_CFLAGS) -o flrtstl flrtstl.c libml.a
-
-flrtstl.o: flrtstl.c
-
-nantst: nantst.c libml.a
-	$(TARGET_CC) $(TARGET_CFLAGS) -o nantst nantst.c libml.a
-
-nantst.o: nantst.c
-
-testvect: testvect.o libml.a
-	$(TARGET_CC) $(TARGET_CFLAGS) -o testvect testvect.o libml.a
-
-testvect.o: testvect.c
-	$(TARGET_CC) -g -c -o testvect.o testvect.c
-
-monotl: monotl.o libml.a
-	$(TARGET_CC) $(TARGET_CFLAGS) -o monotl monotl.o libml.a
-
-monotl.o: monotl.c
-	$(TARGET_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
deleted file mode 100644
index 30fcaad36..000000000
--- a/libm/ldouble/README.txt
+++ /dev/null
@@ -1,3502 +0,0 @@
-/*							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: