[src/trunk]: src/lib/libm Add tgamma{, f} from FreeBSD via rudolf, netbsd at ...

[christos <christos%NetBSD.org@localhost>]

details:   https://anonhg.NetBSD.org/src/rev/7fb0a0340bf9
branches:  trunk
changeset: 779127:7fb0a0340bf9
user:      christos <christos%NetBSD.org@localhost>
date:      Sat May 05 17:54:13 2012 +0000

description:
Add tgamma{,f} from FreeBSD via  rudolf, netbsd at eq dot cz

diffstat:

 lib/libm/Makefile           |    9 +-
 lib/libm/man/lgamma.3       |   39 +++-
 lib/libm/src/b_exp.c        |  137 ++++++++++++++
 lib/libm/src/b_log.c        |  406 ++++++++++++++++++++++++++++++++++++++++++++
 lib/libm/src/b_tgamma.c     |  317 ++++++++++++++++++++++++++++++++++
 lib/libm/src/math_private.h |   31 +++-
 lib/libm/src/s_tgammaf.c    |   47 +++++
 7 files changed, 978 insertions(+), 8 deletions(-)

diffs (truncated from 1088 to 300 lines):

diff -r 46e6b92e7295 -r 7fb0a0340bf9 lib/libm/Makefile
--- a/lib/libm/Makefile Sat May 05 17:52:27 2012 +0000
+++ b/lib/libm/Makefile Sat May 05 17:54:13 2012 +0000
@@ -1,4 +1,4 @@
-#  $NetBSD: Makefile,v 1.123 2012/04/04 10:59:46 joerg Exp $
+#  $NetBSD: Makefile,v 1.124 2012/05/05 17:54:13 christos Exp $
 #
 #  @(#)Makefile 5.1beta 93/09/24
 #
@@ -129,7 +129,8 @@
 CPPFLAGS+=-DLIBM_SCCS
 
 LIB=   m
-COMMON_SRCS+= e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
+COMMON_SRCS+= b_exp.c b_log.c b_tgamma.c \
+       e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
        e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
        e_expf.c e_fmod.c e_fmodf.c e_hypot.c e_hypotf.c e_j0.c e_j0f.c \
        e_j1.c e_j1f.c e_jn.c e_jnf.c e_lgamma_r.c e_lgammaf_r.c e_log.c \
@@ -148,7 +149,7 @@
        s_matherr.c s_modff.c s_nextafter.c \
        s_nextafterf.c s_remquo.c s_remquof.c s_rint.c s_rintf.c s_round.c s_roundf.c s_scalbn.c \
        s_scalbnf.c s_scalbnl.c s_signgam.c s_significand.c s_significandf.c s_sin.c \
-       s_sinf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c \
+       s_sinf.c s_tan.c s_tanf.c s_tanh.c s_tanhf.c s_tgammaf.c s_trunc.c s_truncf.c \
        w_acos.c w_acosf.c w_acosh.c w_acoshf.c w_asin.c w_asinf.c w_atan2.c \
        w_atan2f.c w_atanh.c w_atanhf.c w_cosh.c w_coshf.c \
        w_drem.c w_dremf.c w_exp.c w_expf.c w_fmod.c w_fmodf.c w_gamma.c \
@@ -269,7 +270,7 @@
        j0.3 y0.3 j0.3 y0f.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3 j0.3 ynf.3
 MLINKS+=lgamma.3 lgammaf.3 lgamma.3 lgamma_r.3 lgamma.3 lgammaf_r.3 \
        lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 gamma_r.3 \
-       lgamma.3 gammaf_r.3
+       lgamma.3 gammaf_r.3 lgamma.3 tgamma.3 lgamma.3 tgammaf.3
 MLINKS+=nextafter.3 nextafterf.3 \
        nextafter.3 nextafterl.3 \
        nextafter.3 nexttoward.3
diff -r 46e6b92e7295 -r 7fb0a0340bf9 lib/libm/man/lgamma.3
--- a/lib/libm/man/lgamma.3     Sat May 05 17:52:27 2012 +0000
+++ b/lib/libm/man/lgamma.3     Sat May 05 17:54:13 2012 +0000
@@ -26,9 +26,9 @@
 .\" SUCH DAMAGE.
 .\"
 .\"     from: @(#)lgamma.3     6.6 (Berkeley) 12/3/92
-.\"    $NetBSD: lgamma.3,v 1.21 2003/08/07 16:44:48 agc Exp $
+.\"    $NetBSD: lgamma.3,v 1.22 2012/05/05 17:54:13 christos Exp $
 .\"
-.Dd December 3, 1992
+.Dd May 4, 2012
 .Dt LGAMMA 3
 .Os
 .Sh NAME
@@ -39,7 +39,9 @@
 .Nm gamma ,
 .Nm gammaf ,
 .Nm gamma_r ,
-.Nm gammaf_r
+.Nm gammaf_r ,
+.Nm tgamma ,
+.Nm tgammaf
 .Nd log gamma function
 .Sh LIBRARY
 .Lb libm
@@ -64,6 +66,10 @@
 .Fn gamma_r "double x" "int *sign"
 .Ft float
 .Fn gammaf_r "float x" "int *sign"
+.Ft double
+.Fn tgamma "double x"
+.Ft float
+.Fn tgammaf "float x"
 .Sh DESCRIPTION
 .Fn lgamma x
 .if t \{\
@@ -90,6 +96,26 @@
 argument and
 .Fa signgam
 is not modified.
+.Pp
+The
+.Fn tgamma x
+and
+.Fn tgammaf x
+functions return \(*G(x), with no effect on
+.Fa signgam .
+.Pp
+.Fn gamma ,
+.Fn gammaf ,
+.Fn gamma_r ,
+and
+.Fn gammaf_r
+are deprecated aliases for
+.Fn lgamma ,
+.Fn lgammaf ,
+.Fn lgamma_r ,
+and
+.Fn lgammaf_r ,
+respectively.
 .Sh IDIOSYNCRASIES
 Do not use the expression
 .Dq Li signgam\(**exp(lgamma(x))
@@ -107,6 +133,9 @@
 returns appropriate values unless an argument is out of range.
 Overflow will occur for sufficiently large positive values, and
 non-positive integers.
+For large non-integer negative values,
+.Fn tgamma
+will underflow.
 On the
 .Tn VAX ,
 the reserved operator is returned,
@@ -121,3 +150,7 @@
 .Nm lgamma
 function appeared in
 .Bx 4.3 .
+The
+.Fn tgamma
+function appeared in
+.Nx 6.0 .
diff -r 46e6b92e7295 -r 7fb0a0340bf9 lib/libm/src/b_exp.c
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/lib/libm/src/b_exp.c      Sat May 05 17:54:13 2012 +0000
@@ -0,0 +1,137 @@
+/*     $NetBSD: b_exp.c,v 1.1 2012/05/05 17:54:14 christos Exp $       */
+
+/*
+ * Copyright (c) 1985, 1993
+ *     The Regents of the University of California.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *     This product includes software developed by the University of
+ *     California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/* @(#)exp.c   8.1 (Berkeley) 6/4/93 */
+#include <sys/cdefs.h>
+#if 0
+__FBSDID("$FreeBSD: release/9.0.0/lib/msun/bsdsrc/b_exp.c 176449 2008-02-22 02:26:51Z das $");
+#else
+__RCSID("$NetBSD: b_exp.c,v 1.1 2012/05/05 17:54:14 christos Exp $");
+#endif
+
+
+/* EXP(X)
+ * RETURN THE EXPONENTIAL OF X
+ * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
+ * CODED IN C BY K.C. NG, 1/19/85;
+ * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
+ *
+ * Required system supported functions:
+ *     scalb(x,n)
+ *     copysign(x,y)
+ *     finite(x)
+ *
+ * Method:
+ *     1. Argument Reduction: given the input x, find r and integer k such
+ *        that
+ *                        x = k*ln2 + r,  |r| <= 0.5*ln2 .
+ *        r will be represented as r := z+c for better accuracy.
+ *
+ *     2. Compute exp(r) by
+ *
+ *             exp(r) = 1 + r + r*R1/(2-R1),
+ *        where
+ *             R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
+ *
+ *     3. exp(x) = 2^k * exp(r) .
+ *
+ * Special cases:
+ *     exp(INF) is INF, exp(NaN) is NaN;
+ *     exp(-INF)=  0;
+ *     for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *     exp(x) returns the exponential of x nearly rounded. In a test run
+ *     with 1,156,000 random arguments on a VAX, the maximum observed
+ *     error was 0.869 ulps (units in the last place).
+ */
+
+#include "math.h"
+#include "math_private.h"
+
+static const double p1 = 0x1.555555555553ep-3;
+static const double p2 = -0x1.6c16c16bebd93p-9;
+static const double p3 = 0x1.1566aaf25de2cp-14;
+static const double p4 = -0x1.bbd41c5d26bf1p-20;
+static const double p5 = 0x1.6376972bea4d0p-25;
+static const double ln2hi = 0x1.62e42fee00000p-1;
+static const double ln2lo = 0x1.a39ef35793c76p-33;
+static const double lnhuge = 0x1.6602b15b7ecf2p9;
+static const double lntiny = -0x1.77af8ebeae354p9;
+static const double invln2 = 0x1.71547652b82fep0;
+
+/* returns exp(r = x + c) for |c| < |x| with no overlap.  */
+
+double
+__exp__D(double x, double c)
+{
+       double  z,hi,lo;
+       int k;
+
+       if (x != x)     /* x is NaN */
+               return(x);
+       if ( x <= lnhuge ) {
+               if ( x >= lntiny ) {
+
+                   /* argument reduction : x --> x - k*ln2 */
+                       z = invln2*x;
+                       k = z + copysign(.5, x);
+
+                   /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
+
+                       hi=(x-k*ln2hi);                 /* Exact. */
+                       x= hi - (lo = k*ln2lo-c);
+                   /* return 2^k*[1+x+x*c/(2+c)]  */
+                       z=x*x;
+                       c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
+                       c = (x*c)/(2.0-c);
+
+                       return  scalb(1.+(hi-(lo - c)), k);
+               }
+               /* end of x > lntiny */
+
+               else
+                    /* exp(-big#) underflows to zero */
+                    if(finite(x))  return(scalb(1.0,-5000));
+
+                    /* exp(-INF) is zero */
+                    else return(0.0);
+       }
+       /* end of x < lnhuge */
+
+       else
+       /* exp(INF) is INF, exp(+big#) overflows to INF */
+           return( finite(x) ?  scalb(1.0,5000)  : x);
+}
diff -r 46e6b92e7295 -r 7fb0a0340bf9 lib/libm/src/b_log.c
--- /dev/null   Thu Jan 01 00:00:00 1970 +0000
+++ b/lib/libm/src/b_log.c      Sat May 05 17:54:13 2012 +0000
@@ -0,0 +1,406 @@
+/*     $NetBSD: b_log.c,v 1.1 2012/05/05 17:54:14 christos Exp $       */
+
+/*
+ * Copyright (c) 1992, 1993
+ *     The Regents of the University of California.  All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ * 3. All advertising materials mentioning features or use of this software
+ *    must display the following acknowledgement:
+ *     This product includes software developed by the University of
+ *     California, Berkeley and its contributors.
+ * 4. Neither the name of the University nor the names of its contributors
+ *    may be used to endorse or promote products derived from this software
+ *    without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/* @(#)log.c   8.2 (Berkeley) 11/30/93 */