uClibc now has a math library. muahahahaha!

 -Erik

82 files changed, 20650 insertions, 0 deletions

diff --git a/libm/float/Makefile b/libm/float/Makefile
new file mode 100644
index 000000000..389ac1a5d
--- /dev/null
+++ b/libm/float/Makefile
@@ -0,0 +1,62 @@
+# 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= acoshf.c airyf.c asinf.c asinhf.c atanf.c \
+	atanhf.c bdtrf.c betaf.c cbrtf.c chbevlf.c chdtrf.c \
+	clogf.c cmplxf.c constf.c coshf.c dawsnf.c ellief.c \
+	ellikf.c ellpef.c ellpkf.c ellpjf.c expf.c exp2f.c \
+	exp10f.c expnf.c facf.c fdtrf.c floorf.c fresnlf.c \
+	gammaf.c gdtrf.c hypergf.c hyp2f1f.c igamf.c igamif.c \
+	incbetf.c incbif.c i0f.c i1f.c ivf.c j0f.c j1f.c \
+	jnf.c jvf.c k0f.c k1f.c knf.c logf.c log2f.c \
+	log10f.c nbdtrf.c ndtrf.c ndtrif.c pdtrf.c polynf.c \
+	powif.c powf.c psif.c rgammaf.c shichif.c sicif.c \
+	sindgf.c sinf.c sinhf.c spencef.c sqrtf.c stdtrf.c \
+	struvef.c tandgf.c tanf.c tanhf.c ynf.c zetaf.c \
+	zetacf.c polevlf.c setprec.c mtherr.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
+
diff --git a/libm/float/README.txt b/libm/float/README.txt
new file mode 100644
index 000000000..30a10b083
--- /dev/null
+++ b/libm/float/README.txt
@@ -0,0 +1,4721 @@
+/*							acoshf.c
+ *
+ *	Inverse hyperbolic cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, acoshf();
+ *
+ * y = acoshf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns inverse hyperbolic cosine of argument.
+ *
+ * If 1 <= x < 1.5, a polynomial approximation
+ *
+ *	sqrt(z) * P(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         100000      1.8e-7       3.9e-8
+ *    IEEE      1,2000      100000                   3.0e-8
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * acoshf domain      |x| < 1            0.0
+ *
+ */
+
+/*							airy.c
+ *
+ *	Airy function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, ai, aip, bi, bip;
+ * int airyf();
+ *
+ * airyf( x, _&ai, _&aip, _&bi, _&bip );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Solution of the differential equation
+ *
+ *	y"(x) = xy.
+ *
+ * The function returns the two independent solutions Ai, Bi
+ * and their first derivatives Ai'(x), Bi'(x).
+ *
+ * Evaluation is by power series summation for small x,
+ * by rational minimax approximations for large x.
+ *
+ *
+ *
+ * ACCURACY:
+ * Error criterion is absolute when function <= 1, relative
+ * when function > 1, except * denotes relative error criterion.
+ * For large negative x, the absolute error increases as x^1.5.
+ * For large positive x, the relative error increases as x^1.5.
+ *
+ * Arithmetic  domain   function  # trials      peak         rms
+ * IEEE        -10, 0     Ai        50000       7.0e-7      1.2e-7
+ * IEEE          0, 10    Ai        50000       9.9e-6*     6.8e-7*
+ * IEEE        -10, 0     Ai'       50000       2.4e-6      3.5e-7
+ * IEEE          0, 10    Ai'       50000       8.7e-6*     6.2e-7*
+ * IEEE        -10, 10    Bi       100000       2.2e-6      2.6e-7
+ * IEEE        -10, 10    Bi'       50000       2.2e-6      3.5e-7
+ *
+ */
+
+/*							asinf.c
+ *
+ *	Inverse circular sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, asinf();
+ *
+ * y = asinf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns radian angle between -pi/2 and +pi/2 whose sine is x.
+ *
+ * A polynomial of the form x + x**3 P(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       100000       2.5e-7       5.0e-8
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * asinf domain        |x| > 1           0.0
+ *
+ */
+/*							acosf()
+ *
+ *	Inverse circular cosine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, acosf();
+ *
+ * y = acosf( 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      100000       1.4e-7      4.2e-8
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * acosf domain        |x| > 1           0.0
+ */
+
+/*							asinhf.c
+ *
+ *	Inverse hyperbolic sine
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, asinhf();
+ *
+ * y = asinhf( 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        100000       2.4e-7      4.1e-8
+ *
+ */
+
+/*							atanf.c
+ *
+ *	Inverse circular tangent
+ *      (arctangent)
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, atanf();
+ *
+ * y = atanf( 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 ).  A polynomial approximates
+ * the function in this basic interval.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -10, 10     100000      1.9e-7      4.1e-8
+ *
+ */
+/*							atan2f()
+ *
+ *	Quadrant correct inverse circular tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, z, atan2f();
+ *
+ * z = atan2f( 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     100000      1.9e-7      4.1e-8
+ * See atan.c.
+ *
+ */
+
+/*							atanhf.c
+ *
+ *	Inverse hyperbolic tangent
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, atanhf();
+ *
+ * y = atanhf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns inverse hyperbolic tangent of argument in the range
+ * MINLOGF to MAXLOGF.
+ *
+ * If |x| < 0.5, a polynomial approximation is used.
+ * Otherwise,
+ *        atanh(x) = 0.5 * log( (1+x)/(1-x) ).
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      -1,1        100000      1.4e-7      3.1e-8
+ *
+ */
+
+/*							bdtrf.c
+ *
+ *	Binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrf();
+ *
+ * y = bdtrf( 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:
+ *
+ *        Relative error (p varies from 0 to 1):
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,100       2000       6.9e-5      1.1e-5
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtrf domain        k < 0            0.0
+ *                     n < k
+ *                     x < 0, x > 1
+ *
+ */
+/*							bdtrcf()
+ *
+ *	Complemented binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrcf();
+ *
+ * y = bdtrcf( 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:
+ *
+ *        Relative error (p varies from 0 to 1):
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,100       2000       6.0e-5      1.2e-5
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtrcf domain     x<0, x>1, n<k       0.0
+ */
+/*							bdtrif()
+ *
+ *	Inverse binomial distribution
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * int k, n;
+ * float p, y, bdtrif();
+ *
+ * p = bdtrf( 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:
+ *
+ *        Relative error (p varies from 0 to 1):
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,100       2000       3.5e-5      3.3e-6
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * bdtrif domain    k < 0, n <= k         0.0
+ *                  x < 0, x > 1
+ *
+ */
+
+/*							betaf.c
+ *
+ *	Beta function
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float a, b, y, betaf();
+ *
+ * y = betaf( a, b );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ *                   -     -
+ *                  | (a) | (b)
+ * beta( a, b )  =  -----------.
+ *                     -
+ *                    | (a+b)
+ *
+ * For large arguments the logarithm of the function is
+ * evaluated using lgam(), then exponentiated.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE       0,30       10000       4.0e-5      6.0e-6
+ *    IEEE       -20,0      10000       4.9e-3      5.4e-5
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition          value returned
+ * betaf overflow   log(beta) > MAXLOG       0