From aa5ad1f90843fb5468a46b2a8cfac339c2bd1499 Mon Sep 17 00:00:00 2001 From: naddy Date: Wed, 18 Jul 2018 20:21:12 +0000 Subject: [PATCH] Remove the unused leftovers of the 4.4BSD libm, which was only used on non-IEEE platforms. Since the VAX port was discontinued, all our remaining architectures use IEEE floating point, as will any future ones. ok millert@ tb@ --- lib/libm/Makefile | 3 +- lib/libm/noieee_src/mathimpl.h | 45 --- lib/libm/noieee_src/n_acosh.c | 88 ------ lib/libm/noieee_src/n_asincos.c | 173 ----------- lib/libm/noieee_src/n_asinh.c | 89 ------ lib/libm/noieee_src/n_atan.c | 93 ------ lib/libm/noieee_src/n_atan2.c | 247 ---------------- lib/libm/noieee_src/n_atanh.c | 79 ----- lib/libm/noieee_src/n_cbrt.c | 115 -------- lib/libm/noieee_src/n_cosh.c | 118 -------- lib/libm/noieee_src/n_erf.c | 388 ------------------------ lib/libm/noieee_src/n_exp.c | 174 ----------- lib/libm/noieee_src/n_exp__E.c | 124 -------- lib/libm/noieee_src/n_expm1.c | 149 ---------- lib/libm/noieee_src/n_floor.c | 174 ----------- lib/libm/noieee_src/n_fmod.c | 80 ----- lib/libm/noieee_src/n_hypot.c | 195 ------------ lib/libm/noieee_src/n_j0.c | 431 --------------------------- lib/libm/noieee_src/n_j1.c | 439 --------------------------- lib/libm/noieee_src/n_jn.c | 306 ------------------- lib/libm/noieee_src/n_lgamma.c | 327 --------------------- lib/libm/noieee_src/n_log.c | 482 ------------------------------ lib/libm/noieee_src/n_log10.c | 83 ------ lib/libm/noieee_src/n_log1p.c | 154 ---------- lib/libm/noieee_src/n_log__L.c | 86 ------ lib/libm/noieee_src/n_lrint.c | 57 ---- lib/libm/noieee_src/n_nan.c | 29 -- lib/libm/noieee_src/n_pow.c | 205 ------------- lib/libm/noieee_src/n_sincos.c | 110 ------- lib/libm/noieee_src/n_sinh.c | 106 ------- lib/libm/noieee_src/n_support.c | 505 -------------------------------- lib/libm/noieee_src/n_tan.c | 74 ----- lib/libm/noieee_src/n_tanh.c | 95 ------ lib/libm/noieee_src/n_tgamma.c | 363 ----------------------- 34 files changed, 1 insertion(+), 6185 deletions(-) delete mode 100644 lib/libm/noieee_src/mathimpl.h delete mode 100644 lib/libm/noieee_src/n_acosh.c delete mode 100644 lib/libm/noieee_src/n_asincos.c delete mode 100644 lib/libm/noieee_src/n_asinh.c delete mode 100644 lib/libm/noieee_src/n_atan.c delete mode 100644 lib/libm/noieee_src/n_atan2.c delete mode 100644 lib/libm/noieee_src/n_atanh.c delete mode 100644 lib/libm/noieee_src/n_cbrt.c delete mode 100644 lib/libm/noieee_src/n_cosh.c delete mode 100644 lib/libm/noieee_src/n_erf.c delete mode 100644 lib/libm/noieee_src/n_exp.c delete mode 100644 lib/libm/noieee_src/n_exp__E.c delete mode 100644 lib/libm/noieee_src/n_expm1.c delete mode 100644 lib/libm/noieee_src/n_floor.c delete mode 100644 lib/libm/noieee_src/n_fmod.c delete mode 100644 lib/libm/noieee_src/n_hypot.c delete mode 100644 lib/libm/noieee_src/n_j0.c delete mode 100644 lib/libm/noieee_src/n_j1.c delete mode 100644 lib/libm/noieee_src/n_jn.c delete mode 100644 lib/libm/noieee_src/n_lgamma.c delete mode 100644 lib/libm/noieee_src/n_log.c delete mode 100644 lib/libm/noieee_src/n_log10.c delete mode 100644 lib/libm/noieee_src/n_log1p.c delete mode 100644 lib/libm/noieee_src/n_log__L.c delete mode 100644 lib/libm/noieee_src/n_lrint.c delete mode 100644 lib/libm/noieee_src/n_nan.c delete mode 100644 lib/libm/noieee_src/n_pow.c delete mode 100644 lib/libm/noieee_src/n_sincos.c delete mode 100644 lib/libm/noieee_src/n_sinh.c delete mode 100644 lib/libm/noieee_src/n_support.c delete mode 100644 lib/libm/noieee_src/n_tan.c delete mode 100644 lib/libm/noieee_src/n_tanh.c delete mode 100644 lib/libm/noieee_src/n_tgamma.c diff --git a/lib/libm/Makefile b/lib/libm/Makefile index 2921d8237ee..2031f8b0df3 100644 --- a/lib/libm/Makefile +++ b/lib/libm/Makefile @@ -1,4 +1,4 @@ -# $OpenBSD: Makefile,v 1.118 2018/03/10 20:52:58 kettenis Exp $ +# $OpenBSD: Makefile,v 1.119 2018/07/18 20:21:12 naddy Exp $ # $NetBSD: Makefile,v 1.28 1995/11/20 22:06:19 jtc Exp $ # # @(#)Makefile 5.1beta 93/09/24 @@ -67,7 +67,6 @@ ARCH_SRCS = e_sqrtl.c .PATH: ${.CURDIR}/man .PATH: ${.CURDIR}/src -.PATH: ${.CURDIR}/noieee_src LIB= m COMMON_SRCS = b_exp__D.c b_log__D.c b_tgamma.c \ diff --git a/lib/libm/noieee_src/mathimpl.h b/lib/libm/noieee_src/mathimpl.h deleted file mode 100644 index 163843ab186..00000000000 --- a/lib/libm/noieee_src/mathimpl.h +++ /dev/null @@ -1,45 +0,0 @@ -/* $OpenBSD: mathimpl.h,v 1.10 2009/04/11 20:03:21 martynas Exp $ */ -/* $NetBSD: mathimpl.h,v 1.1 1995/10/10 23:36:31 ragge Exp $ */ -/* - * Copyright (c) 1988, 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. 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. - * - * @(#)mathimpl.h 8.1 (Berkeley) 6/4/93 - */ - -/* - * Functions internal to the math package, yet not static. - */ -extern double __exp__E(); -extern double __log__L(); - -struct Double { - double a, b; -}; - -double __exp__D(double, double); -struct Double __log__D(double); diff --git a/lib/libm/noieee_src/n_acosh.c b/lib/libm/noieee_src/n_acosh.c deleted file mode 100644 index 46b2914339f..00000000000 --- a/lib/libm/noieee_src/n_acosh.c +++ /dev/null @@ -1,88 +0,0 @@ -/* $OpenBSD: n_acosh.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_acosh.c,v 1.1 1995/10/10 23:36:33 ragge 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. 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. - */ - -/* ACOSH(X) - * RETURN THE INVERSE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/6/85, 3/24/85, 4/16/85, 8/17/85. - * - * Required system supported functions : - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log1p(x)+ln2, if (x > 1.0E20); else - * acosh(x) := log1p( sqrt(x-1) * (sqrt(x-1) + sqrt(x+1)) ) . - * These formulae avoid the over/underflow complication. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - * - * Accuracy: - * acosh(x) returns the exact inverse hyperbolic cosine of x nearly - * rounded. In a test run with 512,000 random arguments on a VAX, the - * maximum observed error was 3.30 ulps (units of the last place) at - * x=1.0070493753568216 . - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double ln2hi = 6.9314718055829871446E-1; -static const double ln2lo = 1.6465949582897081279E-12; - -double -acosh(double x) -{ - double t,big=1.E20; /* big+1==big */ - - if (isnan(x)) - return (x); - - /* return log1p(x) + log(2) if x is large */ - if(x>big) {t=log1p(x)+ln2lo; return(t+ln2hi);} - - t=sqrt(x-1.0); - return(log1p(t*(t+sqrt(x+1.0)))); -} diff --git a/lib/libm/noieee_src/n_asincos.c b/lib/libm/noieee_src/n_asincos.c deleted file mode 100644 index a6fe371d39a..00000000000 --- a/lib/libm/noieee_src/n_asincos.c +++ /dev/null @@ -1,173 +0,0 @@ -/* $OpenBSD: n_asincos.c,v 1.15 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_asincos.c,v 1.1 1995/10/10 23:36:34 ragge 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. 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. - */ - -/* ASIN(X) - * RETURNS ARC SINE OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method: - * asin(x) = atan2(x,sqrt(1-x*x)); for better accuracy, 1-x*x is - * computed as follows - * 1-x*x if x < 0.5, - * 2*(1-|x|)-(1-|x|)*(1-|x|) if x >= 0.5. - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * asin(x) returns (PI/pi) * (the exact arc sine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.06 ulps. (comparing against (PI/pi)*(exact asin(x))); - * - * 2) If atan2() uses true pi, then - * - * asin(x) returns the exact asin(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 1.99 ulps. - */ - -#include - -#include "mathimpl.h" - -double -asin(double x) -{ - double s, t, one = 1.0; - - if (isnan(x)) - return (x); - - s=copysign(x,one); - if(s <= 0.5) - return(atan2(x,sqrt(one-x*x))); - else - { t=one-s; s=t+t; return(atan2(x,sqrt(s-t*t))); } - -} -DEF_STD(asin); -LDBL_CLONE(asin); - -/* ACOS(X) - * RETURNS ARC COS OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required system supported functions: - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * atan2(y,x) - * - * Method: - * ________ - * / 1 - x - * acos(x) = 2*atan2( / -------- , 1 ) . - * \/ 1 + x - * - * Special cases: - * if x is NaN, return x itself; - * if |x|>1, return NaN. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * acos(x) returns (PI/pi) * (the exact arc cosine of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.07 ulps. (comparing against (PI/pi)*(exact acos(x))); - * - * 2) If atan2() uses true pi, then - * - * acos(x) returns the exact acos(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 2.15 ulps. - */ - -double -acos(double x) -{ - double t, one = 1.0; - - if (isnan(x)) - return (x); - - if( x != -1.0) - t=atan2(sqrt((one-x)/(one+x)),one); - else - t=atan2(one,0.0); /* t = PI/2 */ - return(t+t); -} -DEF_STD(acos); -LDBL_UNUSED_CLONE(acos); diff --git a/lib/libm/noieee_src/n_asinh.c b/lib/libm/noieee_src/n_asinh.c deleted file mode 100644 index 6307c9356b4..00000000000 --- a/lib/libm/noieee_src/n_asinh.c +++ /dev/null @@ -1,89 +0,0 @@ -/* $OpenBSD: n_asinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_asinh.c,v 1.1 1995/10/10 23:36:35 ragge 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. 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. - */ - -/* ASINH(X) - * RETURN THE INVERSE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 2/16/85; - * REVISED BY K.C. NG on 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * sqrt(x) - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log1p(x)+ln2)) if sqrt(1+x*x)=x, else - * := sign(x)*log1p(|x| + |x|/(1/|x| + sqrt(1+(1/|x|)^2)) ) - * - * Accuracy: - * asinh(x) returns the exact inverse hyperbolic sine of x nearly rounded. - * In a test run with 52,000 random arguments on a VAX, the maximum - * observed error was 1.58 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double ln2hi = 6.9314718055829871446E-1; -static const double ln2lo = 1.6465949582897081279E-12; - -double -asinh(double x) -{ - double t,s; - static const double small=1.0E-10, /* fl(1+small*small) == 1 */ - big =1.0E20, /* fl(1+big) == big */ - one =1.0 ; - - if (isnan(x)) - return (x); - - if((t=copysign(x,one))>small) - if(t big */ - {s=log1p(t)+ln2lo; return(copysign(s+ln2hi,x));} - else /* if |x| < small */ - return(x); -} diff --git a/lib/libm/noieee_src/n_atan.c b/lib/libm/noieee_src/n_atan.c deleted file mode 100644 index 79747e2c170..00000000000 --- a/lib/libm/noieee_src/n_atan.c +++ /dev/null @@ -1,93 +0,0 @@ -/* $OpenBSD: n_atan.c,v 1.14 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_atan.c,v 1.1 1995/10/10 23:36:36 ragge 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. 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. - */ - -/* ATAN(X) - * RETURNS ARC TANGENT OF X - * DOUBLE PRECISION (IEEE DOUBLE 53 bits, VAX D FORMAT 56 bits) - * CODED IN C BY K.C. NG, 4/16/85, REVISED ON 6/10/85. - * - * Required kernel function: - * atan2(y,x) - * - * Method: - * atan(x) = atan2(x,1.0). - * - * Special case: - * if x is NaN, return x itself. - * - * Accuracy: - * 1) If atan2() uses machine PI, then - * - * atan(x) returns (PI/pi) * (the exact arc tangent of x) nearly rounded; - * and PI is the exact pi rounded to machine precision (see atan2 for - * details): - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with more than 200,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.86 ulps. (comparing against (PI/pi)*(exact atan(x))). - * - * 2) If atan2() uses true pi, then - * - * atan(x) returns the exact atan(x) with error below about 2 ulps. - * - * In a test run with more than 1,024,000 random arguments on a VAX, the - * maximum observed error in ulps (units in the last place) was - * 0.85 ulps. - */ - -#include - -double -atan(double x) -{ - double one=1.0; - return(atan2(x,one)); -} -DEF_STD(atan); -LDBL_CLONE(atan); - -float -atanf(float x) -{ - float one=1.0f; - return(atan2f(x,one)); -} -DEF_STD(atanf); diff --git a/lib/libm/noieee_src/n_atan2.c b/lib/libm/noieee_src/n_atan2.c deleted file mode 100644 index e99441349bc..00000000000 --- a/lib/libm/noieee_src/n_atan2.c +++ /dev/null @@ -1,247 +0,0 @@ -/* $OpenBSD: n_atan2.c,v 1.20 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_atan2.c,v 1.1 1995/10/10 23:36:37 ragge 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. 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. - */ - -/* ATAN2(Y,X) - * RETURN ARG (X+iY) - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85, 6/29/85. - * - * Required system supported functions : - * copysign(x,y) - * scalbn(x,y) - * logb(x) - * - * Method : - * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): - * ARG (x+iy) = arctan(y/x) ... if x > 0, - * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, - * 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument - * is further reduced to one of the following intervals and the - * arctangent of y/x is evaluated by the corresponding formula: - * - * [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) - * [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) - * [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) - * [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) - * [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) - * - * Special cases: - * Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). - * - * ARG( NAN , (anything) ) is NaN; - * ARG( (anything), NaN ) is NaN; - * ARG(+(anything but NaN), +-0) is +-0 ; - * ARG(-(anything but NaN), +-0) is +-PI ; - * ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; - * ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; - * ARG( -INF,+-(anything but INF and NaN) ) is +-PI; - * ARG( +INF,+-INF ) is +-PI/4 ; - * ARG( -INF,+-INF ) is +-3PI/4; - * ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; - * - * Accuracy: - * atan2(y,x) returns (PI/pi) * the exact ARG (x+iy) nearly rounded, - * where - * - * in decimal: - * pi = 3.141592653589793 23846264338327 ..... - * 53 bits PI = 3.141592653589793 115997963 ..... , - * 56 bits PI = 3.141592653589793 227020265 ..... , - * - * in hexadecimal: - * pi = 3.243F6A8885A308D313198A2E.... - * 53 bits PI = 3.243F6A8885A30 = 2 * 1.921FB54442D18 error=.276ulps - * 56 bits PI = 3.243F6A8885A308 = 4 * .C90FDAA22168C2 error=.206ulps - * - * In a test run with 356,000 random argument on [-1,1] * [-1,1] on a - * VAX, the maximum observed error was 1.41 ulps (units of the last place) - * compared with (PI/pi)*(the exact ARG(x+iy)). - * - * Note: - * We use machine PI (the true pi rounded) in place of the actual - * value of pi for all the trig and inverse trig functions. In general, - * if trig is one of sin, cos, tan, then computed trig(y) returns the - * exact trig(y*pi/PI) nearly rounded; correspondingly, computed arctrig - * returns the exact arctrig(y)*PI/pi nearly rounded. These guarantee the - * trig functions have period PI, and trig(arctrig(x)) returns x for - * all critical values x. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include - -#include "mathimpl.h" - -static const double athfhi = 4.6364760900080611433E-1; -static const double athflo = 1.9338828231967579916E-19; -static const double PIo4 = 7.8539816339744830676E-1; -static const double at1fhi = 9.8279372324732906796E-1; -static const double at1flo = -3.5540295636764633916E-18; -static const double PIo2 = 1.5707963267948966135E0; -static const double PI = 3.1415926535897932270E0; -static const double a1 = 3.3333333333333473730E-1; -static const double a2 = -2.0000000000017730678E-1; -static const double a3 = 1.4285714286694640301E-1; -static const double a4 = -1.1111111135032672795E-1; -static const double a5 = 9.0909091380563043783E-2; -static const double a6 = -7.6922954286089459397E-2; -static const double a7 = 6.6663180891693915586E-2; -static const double a8 = -5.8772703698290408927E-2; -static const double a9 = 5.2170707402812969804E-2; -static const double a10 = -4.4895863157820361210E-2; -static const double a11 = 3.3006147437343875094E-2; -static const double a12 = -1.4614844866464185439E-2; - -float -atan2f(float x, float y) -{ - return (float)atan2((double)x, (double)y); -} - -double -atan2(double y, double x) -{ - static const double zero=0, one=1, small=1.0E-9, big=1.0E18; - double t,z,signy,signx,hi,lo; - int k,m; - - /* if x or y is NAN */ - if (isnan(x)) - return (x); - if (isnan(y)) - return (y); - - /* copy down the sign of y and x */ - signy = copysign(one,y) ; - signx = copysign(one,x) ; - - /* if x is 1.0, goto begin */ - if(x==1) { y=copysign(y,one); t=y; if(isfinite(t)) goto begin;} - - /* when y = 0 */ - if(y==zero) return((signx==one)?y:copysign(PI,signy)); - - /* when x = 0 */ - if(x==zero) return(copysign(PIo2,signy)); - - /* when x is INF */ - if(!isfinite(x)) - if(!isfinite(y)) - return(copysign((signx==one)?PIo4:3*PIo4,signy)); - else - return(copysign((signx==one)?zero:PI,signy)); - - /* when y is INF */ - if(!isfinite(y)) return(copysign(PIo2,signy)); - - /* compute y/x */ - x=copysign(x,one); - y=copysign(y,one); - if((m=(k=logb(y))-logb(x)) > 60) t=big+big; - else if(m < -80 ) t=y/x; - else { t = y/x ; y = scalbn(y,-k); x=scalbn(x,-k); } - - /* begin argument reduction */ -begin: - if (t < 2.4375) { - - /* truncate 4(t+1/16) to integer for branching */ - k = 4 * (t+0.0625); - switch (k) { - - /* t is in [0,7/16] */ - case 0: - case 1: - if (t < small) { - if (big + small > 0.0) /* raise inexact flag */ - return (copysign((signx>zero)?t:PI-t,signy)); - } - - hi = zero; lo = zero; break; - - /* t is in [7/16,11/16] */ - case 2: - hi = athfhi; lo = athflo; - z = x+x; - t = ( (y+y) - x ) / ( z + y ); break; - - /* t is in [11/16,19/16] */ - case 3: - case 4: - hi = PIo4; lo = zero; - t = ( y - x ) / ( x + y ); break; - - /* t is in [19/16,39/16] */ - default: - hi = at1fhi; lo = at1flo; - z = y-x; y=y+y+y; t = x+x; - t = ( (z+z)-x ) / ( t + y ); break; - } - } - /* end of if (t < 2.4375) */ - - else - { - hi = PIo2; lo = zero; - - /* t is in [2.4375, big] */ - if (t <= big) t = - x / y; - - /* t is in [big, INF] */ - else { - if (big + small > 0.0) /* raise inexact flag */ - t = zero; - } - } - /* end of argument reduction */ - - /* compute atan(t) for t in [-.4375, .4375] */ - z = t*t; -#if defined(__vax__) - z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ - z*(a9+z*(a10+z*(a11+z*a12)))))))))))); -#else /* defined(__vax__) */ - z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ - z*(a9+z*(a10+z*a11))))))))))); -#endif /* defined(__vax__) */ - z = lo - z; z += t; z += hi; - - return(copysign((signx>zero)?z:PI-z,signy)); -} -DEF_STD(atan2); -LDBL_CLONE(atan2); diff --git a/lib/libm/noieee_src/n_atanh.c b/lib/libm/noieee_src/n_atanh.c deleted file mode 100644 index 513065aca5e..00000000000 --- a/lib/libm/noieee_src/n_atanh.c +++ /dev/null @@ -1,79 +0,0 @@ -/* $OpenBSD: n_atanh.c,v 1.7 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_atanh.c,v 1.1 1995/10/10 23:36:38 ragge 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. 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. - */ - -/* ATANH(X) - * RETURN THE HYPERBOLIC ARC TANGENT OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/7/85, 3/7/85, 8/18/85. - * - * Required kernel function: - * log1p(x) ...return log(1+x) - * - * Method: - * Return - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - * Accuracy: - * atanh(x) returns the exact hyperbolic arc tangent of x nearly rounded. - * In a test run with 512,000 random arguments on a VAX, the maximum - * observed error was 1.87 ulps (units in the last place) at - * x= -3.8962076028810414000e-03. - */ - -#include "math.h" -#include "mathimpl.h" - -#if defined(__vax__) -#include -#endif /* defined(__vax__) */ - -double -atanh(double x) -{ - double z; - z = copysign(0.5,x); - x = copysign(x,1.0); -#if defined(__vax__) - if (x == 1.0) { - return(copysign(1.0,z)*infnan(ERANGE)); /* sign(x)*INF */ - } -#endif /* defined(__vax__) */ - x = x/(1.0-x); - return( z*log1p(x+x) ); -} diff --git a/lib/libm/noieee_src/n_cbrt.c b/lib/libm/noieee_src/n_cbrt.c deleted file mode 100644 index eb2c1e1eaac..00000000000 --- a/lib/libm/noieee_src/n_cbrt.c +++ /dev/null @@ -1,115 +0,0 @@ -/* $OpenBSD: n_cbrt.c,v 1.8 2013/07/15 04:08:26 espie Exp $ */ -/* $NetBSD: n_cbrt.c,v 1.1 1995/10/10 23:36:40 ragge 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. 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. - */ - -#include - -/* kahan's cube root (53 bits IEEE double precision) - * for IEEE machines only - * coded in C by K.C. Ng, 4/30/85 - * - * Accuracy: - * better than 0.667 ulps according to an error analysis. Maximum - * error observed was 0.666 ulps in an 1,000,000 random arguments test. - * - * Warning: this code is semi machine dependent; the ordering of words in - * a floating point number must be known in advance. I assume that the - * long interger at the address of a floating point number will be the - * leading 32 bits of that floating point number (i.e., sign, exponent, - * and the 20 most significant bits). - * On a National machine, it has different ordering; therefore, this code - * must be compiled with flag -DNATIONAL. - */ -#if !defined(__vax__) - -static const unsigned long - B1 = 715094163, /* B1 = (682-0.03306235651)*2**20 */ - B2 = 696219795; /* B2 = (664-0.03306235651)*2**20 */ -static const double - C= 19./35., - D= -864./1225., - E= 99./70., - F= 45./28., - G= 5./14.; - -float -cbrtf(float x) -{ - return (float)cbrt((double) x); -} - -double -cbrt(double x) -{ - double r,s,t=0.0,w; - unsigned long *px = (unsigned long *) &x, - *pt = (unsigned long *) &t, - mexp,sign; - const int n0=0,n1=1; - - mexp=px[n0]&0x7ff00000; - if(mexp==0x7ff00000) return(x); /* cbrt(NaN,INF) is itself */ - if(x==0.0) return(x); /* cbrt(0) is itself */ - - sign=px[n0]&0x80000000; /* sign= sign(x) */ - px[n0] ^= sign; /* x=|x| */ - - - /* rough cbrt to 5 bits */ - if(mexp==0) /* subnormal number */ - {pt[n0]=0x43500000; /* set t= 2**54 */ - t*=x; pt[n0]=pt[n0]/3+B2; - } - else - pt[n0]=px[n0]/3+B1; - - - /* new cbrt to 23 bits, may be implemented in single precision */ - r=t*t/x; - s=C+r*t; - t*=G+F/(s+E+D/s); - - /* chopped to 20 bits and make it larger than cbrt(x) */ - pt[n1]=0; pt[n0]+=0x00000001; - - - /* one step newton iteration to 53 bits with error less than 0.667 ulps */ - s=t*t; /* t*t is exact */ - r=x/s; - w=t+t; - r=(r-t)/(w+r); /* r-t is exact */ - t=t+t*r; - - - /* retore the sign bit */ - pt[n0] |= sign; - return(t); -} -#endif diff --git a/lib/libm/noieee_src/n_cosh.c b/lib/libm/noieee_src/n_cosh.c deleted file mode 100644 index 72fe1853972..00000000000 --- a/lib/libm/noieee_src/n_cosh.c +++ /dev/null @@ -1,118 +0,0 @@ -/* $OpenBSD: n_cosh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_cosh.c,v 1.1 1995/10/10 23:36:42 ragge 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. 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. - */ - -/* COSH(X) - * RETURN THE HYPERBOLIC COSINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/23/85, 3/7/85, 3/29/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalbn(x,N) - * - * Required kernel function: - * exp(x) - * exp__E(x,c) ...return exp(x+c)-1-x for |x|<0.3465 - * - * Method : - * 1. Replace x by |x|. - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= 0.3465 : cosh(x) := 1 + ------------------- - * 2*exp(x) - * - * exp(x) + 1/exp(x) - * 0.3465 <= x <= 22 : cosh(x) := ------------------- - * 2 - * 22 <= x <= lnovfl : cosh(x) := exp(x)/2 - * lnovfl <= x <= lnovfl+log(2) - * : cosh(x) := exp(x)/2 (avoid overflow) - * log(2)+lnovfl < x < INF: overflow to INF - * - * Note: .3465 is a number near one half of ln2. - * - * Special cases: - * cosh(x) is x if x is +INF, -INF, or NaN. - * only cosh(0)=1 is exact for finite x. - * - * Accuracy: - * cosh(x) returns the exact hyperbolic cosine of x nearly rounded. - * In a test run with 768,000 random arguments on a VAX, the maximum - * observed error was 1.23 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double mln2hi = 8.8029691931113054792E1; -static const double mln2lo = -4.9650192275318476525E-16; -static const double lnovfl = 8.8029691931113053016E1; - -#if defined(__vax__) -static max = 126 ; -#else /* defined(__vax__) */ -static max = 1023 ; -#endif /* defined(__vax__) */ - -double -cosh(double x) -{ - static const double half=1.0/2.0, - one=1.0, small=1.0E-18; /* fl(1+small)==1 */ - double t; - - if (isnan(x)) - return (x); - - if((x=copysign(x,one)) <= 22) - if(x<0.3465) - if(x= 1.25 (error 1.7ulp vs ~6ulp) - * Replaced even+odd with direct calculation for x < .84375, - * to avoid destructive cancellation. - * - * Performance of erfc(x): - * In 300000 trials in the range [.83, .84375] the - * maximum observed error was 3.6ulp. - * - * In [.84735,1.25] the maximum observed error was <2.5ulp in - * 100000 runs in the range [1.2, 1.25]. - * - * In [1.25,26] (Not including subnormal results) - * the error is < 1.7ulp. - */ - -/* double erf(double x) - * double erfc(double x) - * x - * 2 |\ - * erf(x) = --------- | exp(-t*t)dt - * sqrt(pi) \| - * 0 - * - * erfc(x) = 1-erf(x) - * - * Method: - * 1. Reduce x to |x| by erf(-x) = -erf(x) - * 2. For x in [0, 0.84375] - * erf(x) = x + x*P(x^2) - * erfc(x) = 1 - erf(x) if x<=0.25 - * = 0.5 + ((0.5-x)-x*P) if x in [0.25,0.84375] - * where - * 2 2 4 20 - * P = P(x ) = (p0 + p1 * x + p2 * x + ... + p10 * x ) - * is an approximation to (erf(x)-x)/x with precision - * - * -56.45 - * | P - (erf(x)-x)/x | <= 2 - * - * - * Remark. The formula is derived by noting - * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) - * and that - * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 - * is close to one. The interval is chosen because the fixed - * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is - * near 0.6174), and by some experiment, 0.84375 is chosen to - * guarantee the error is less than one ulp for erf. - * - * 3. For x in [0.84375,1.25], let s = x - 1, and - * c = 0.84506291151 rounded to single (24 bits) - * erf(x) = c + P1(s)/Q1(s) - * erfc(x) = (1-c) - P1(s)/Q1(s) - * |P1/Q1 - (erf(x)-c)| <= 2**-59.06 - * Remark: here we use the taylor series expansion at x=1. - * erf(1+s) = erf(1) + s*Poly(s) - * = 0.845.. + P1(s)/Q1(s) - * That is, we use rational approximation to approximate - * erf(1+s) - (c = (single)0.84506291151) - * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] - * where - * P1(s) = degree 6 poly in s - * Q1(s) = degree 6 poly in s - * - * 4. For x in [1.25, 2]; [2, 4] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi) -.5z + R(z)/S(z)) - * - * Where z = 1/(x*x), R is degree 9, and S is degree 3; - * - * 5. For x in [4,28] - * erf(x) = 1.0 - tiny - * erfc(x) = (1/x)exp(-x*x-(.5*log(pi)+eps + zP(z)) - * - * Where P is degree 14 polynomial in 1/(x*x). - * - * Notes: - * Here 4 and 5 make use of the asymptotic series - * exp(-x*x) - * erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) ); - * x*sqrt(pi) - * - * where for z = 1/(x*x) - * P(z) ~ z/2*(-1 + z*3/2*(1 + z*5/2*(-1 + z*7/2*(1 +...)))) - * - * Thus we use rational approximation to approximate - * erfc*x*exp(x*x) ~ 1/sqrt(pi); - * - * The error bound for the target function, G(z) for - * the interval - * [4, 28]: - * |eps + 1/(z)P(z) - G(z)| < 2**(-56.61) - * for [2, 4]: - * |R(z)/S(z) - G(z)| < 2**(-58.24) - * for [1.25, 2]: - * |R(z)/S(z) - G(z)| < 2**(-58.12) - * - * 6. For inf > x >= 28 - * erf(x) = 1 - tiny (raise inexact) - * erfc(x) = tiny*tiny (raise underflow) - * - * 7. Special cases: - * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, - * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, - * erfc/erf(NaN) is NaN - */ - -#if defined(__vax__) -#define _IEEE 0 -#define TRUNC(x) (double) (float) (x) -#else -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -static double -tiny = 1e-300, -half = 0.5, -one = 1.0, -two = 2.0, -c = 8.45062911510467529297e-01, /* (float)0.84506291151 */ -/* - * Coefficients for approximation to erf in [0,0.84375] - */ -p0t8 = 1.02703333676410051049867154944018394163280, -p0 = 1.283791670955125638123339436800229927041e-0001, -p1 = -3.761263890318340796574473028946097022260e-0001, -p2 = 1.128379167093567004871858633779992337238e-0001, -p3 = -2.686617064084433642889526516177508374437e-0002, -p4 = 5.223977576966219409445780927846432273191e-0003, -p5 = -8.548323822001639515038738961618255438422e-0004, -p6 = 1.205520092530505090384383082516403772317e-0004, -p7 = -1.492214100762529635365672665955239554276e-0005, -p8 = 1.640186161764254363152286358441771740838e-0006, -p9 = -1.571599331700515057841960987689515895479e-0007, -p10= 1.073087585213621540635426191486561494058e-0008; -/* - * Coefficients for approximation to erf in [0.84375,1.25] - */ -static double -pa0 = -2.362118560752659485957248365514511540287e-0003, -pa1 = 4.148561186837483359654781492060070469522e-0001, -pa2 = -3.722078760357013107593507594535478633044e-0001, -pa3 = 3.183466199011617316853636418691420262160e-0001, -pa4 = -1.108946942823966771253985510891237782544e-0001, -pa5 = 3.547830432561823343969797140537411825179e-0002, -pa6 = -2.166375594868790886906539848893221184820e-0003, -qa1 = 1.064208804008442270765369280952419863524e-0001, -qa2 = 5.403979177021710663441167681878575087235e-0001, -qa3 = 7.182865441419627066207655332170665812023e-0002, -qa4 = 1.261712198087616469108438860983447773726e-0001, -qa5 = 1.363708391202905087876983523620537833157e-0002, -qa6 = 1.198449984679910764099772682882189711364e-0002; -/* - * log(sqrt(pi)) for large x expansions. - * The tail (lsqrtPI_lo) is included in the rational - * approximations. -*/ -static double - lsqrtPI_hi = .5723649429247000819387380943226; -/* - * lsqrtPI_lo = .000000000000000005132975581353913; - * - * Coefficients for approximation to erfc in [2, 4] -*/ -static double -rb0 = -1.5306508387410807582e-010, /* includes lsqrtPI_lo */ -rb1 = 2.15592846101742183841910806188e-008, -rb2 = 6.24998557732436510470108714799e-001, -rb3 = 8.24849222231141787631258921465e+000, -rb4 = 2.63974967372233173534823436057e+001, -rb5 = 9.86383092541570505318304640241e+000, -rb6 = -7.28024154841991322228977878694e+000, -rb7 = 5.96303287280680116566600190708e+000, -rb8 = -4.40070358507372993983608466806e+000, -rb9 = 2.39923700182518073731330332521e+000, -rb10 = -6.89257464785841156285073338950e-001, -sb1 = 1.56641558965626774835300238919e+001, -sb2 = 7.20522741000949622502957936376e+001, -sb3 = 9.60121069770492994166488642804e+001; -/* - * Coefficients for approximation to erfc in [1.25, 2] -*/ -static double -rc0 = -2.47925334685189288817e-007, /* includes lsqrtPI_lo */ -rc1 = 1.28735722546372485255126993930e-005, -rc2 = 6.24664954087883916855616917019e-001, -rc3 = 4.69798884785807402408863708843e+000, -rc4 = 7.61618295853929705430118701770e+000, -rc5 = 9.15640208659364240872946538730e-001, -rc6 = -3.59753040425048631334448145935e-001, -rc7 = 1.42862267989304403403849619281e-001, -rc8 = -4.74392758811439801958087514322e-002, -rc9 = 1.09964787987580810135757047874e-002, -rc10 = -1.28856240494889325194638463046e-003, -sc1 = 9.97395106984001955652274773456e+000, -sc2 = 2.80952153365721279953959310660e+001, -sc3 = 2.19826478142545234106819407316e+001; -/* - * Coefficients for approximation to erfc in [4,28] - */ -static double -rd0 = -2.1491361969012978677e-016, /* includes lsqrtPI_lo */ -rd1 = -4.99999999999640086151350330820e-001, -rd2 = 6.24999999772906433825880867516e-001, -rd3 = -1.54166659428052432723177389562e+000, -rd4 = 5.51561147405411844601985649206e+000, -rd5 = -2.55046307982949826964613748714e+001, -rd6 = 1.43631424382843846387913799845e+002, -rd7 = -9.45789244999420134263345971704e+002, -rd8 = 6.94834146607051206956384703517e+003, -rd9 = -5.27176414235983393155038356781e+004, -rd10 = 3.68530281128672766499221324921e+005, -rd11 = -2.06466642800404317677021026611e+006, -rd12 = 7.78293889471135381609201431274e+006, -rd13 = -1.42821001129434127360582351685e+007; - -double -erf(double x) -{ - double R, S, P, Q, ax, s, y, z, r; - if(!isfinite(x)) { /* erf(nan)=nan */ - if (isnan(x)) - return(x); - return (x > 0 ? one : -one); /* erf(+/-inf)= +/-1 */ - } - if ((ax = x) < 0) - ax = - ax; - if (ax < .84375) { - if (ax < 3.7e-09) { - if (ax < 1.0e-308) - return 0.125*(8.0*x+p0t8*x); /*avoid underflow */ - return x + p0*x; - } - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - return x + x*(p0+r); - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = fabs(x)-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) - return (c + P/Q); - else - return (-c - P/Q); - } - if (ax >= 6.0) { /* inf>|x|>=6 */ - if (x >= 0.0) - return (one-tiny); - else - return (tiny-one); - } - /* 1.25 <= |x| < 6 */ - z = -ax*ax; - s = -one/z; - if (ax < 2.0) { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - } else { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - } - y = (R/S -.5*s) - lsqrtPI_hi; - z += y; - z = exp(z)/ax; - if (x >= 0) - return (one-z); - else - return (z-one); -} -DEF_STD(erf); - -double -erfc(double x) -{ - double R, S, P, Q, s, ax, y, z, r; - if (!isfinite(x)) { - if (isnan(x)) /* erfc(NaN) = NaN */ - return(x); - else if (x > 0) /* erfc(+-inf)=0,2 */ - return 0.0; - else - return 2.0; - } - if ((ax = x) < 0) - ax = -ax; - if (ax < .84375) { /* |x|<0.84375 */ - if (ax < 1.38777878078144568e-17) /* |x|<2**-56 */ - return one-x; - y = x*x; - r = y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+ - y*(p6+y*(p7+y*(p8+y*(p9+y*p10))))))))); - if (ax < .0625) { /* |x|<2**-4 */ - return (one-(x+x*(p0+r))); - } else { - r = x*(p0+r); - r += (x-half); - return (half - r); - } - } - if (ax < 1.25) { /* 0.84375 <= |x| < 1.25 */ - s = ax-one; - P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); - Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); - if (x>=0) { - z = one-c; return z - P/Q; - } else { - z = c+P/Q; return one+z; - } - } - if (ax >= 28) /* Out of range */ - if (x>0) - return (tiny*tiny); - else - return (two-tiny); - z = ax; - TRUNC(z); - y = z - ax; y *= (ax+z); - z *= -z; /* Here z + y = -x^2 */ - s = one/(-z-y); /* 1/(x*x) */ - if (ax >= 4) { /* 6 <= ax */ - R = s*(rd1+s*(rd2+s*(rd3+s*(rd4+s*(rd5+ - s*(rd6+s*(rd7+s*(rd8+s*(rd9+s*(rd10 - +s*(rd11+s*(rd12+s*rd13)))))))))))); - y += rd0; - } else if (ax >= 2) { - R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+ - s*(rb6+s*(rb7+s*(rb8+s*(rb9+s*rb10))))))))); - S = one+s*(sb1+s*(sb2+s*sb3)); - y += R/S; - R = -.5*s; - } else { - R = rc0+s*(rc1+s*(rc2+s*(rc3+s*(rc4+s*(rc5+ - s*(rc6+s*(rc7+s*(rc8+s*(rc9+s*rc10))))))))); - S = one+s*(sc1+s*(sc2+s*sc3)); - y += R/S; - R = -.5*s; - } - /* return exp(-x^2 - lsqrtPI_hi + R + y)/x; */ - s = ((R + y) - lsqrtPI_hi) + z; - y = (((z-s) - lsqrtPI_hi) + R) + y; - r = __exp__D(s, y)/x; - if (x>0) - return r; - else - return two-r; -} -DEF_STD(erfc); diff --git a/lib/libm/noieee_src/n_exp.c b/lib/libm/noieee_src/n_exp.c deleted file mode 100644 index 986ff59b8f5..00000000000 --- a/lib/libm/noieee_src/n_exp.c +++ /dev/null @@ -1,174 +0,0 @@ -/* $OpenBSD: n_exp.c,v 1.12 2016/09/12 19:47:02 guenther 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. 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(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: - * scalbn(x,n) - * copysign(x,y) - * isfinite(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). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double ln2hi = 6.9314718055829871446E-1; -static const double ln2lo = 1.6465949582897081279E-12; -static const double lnhuge = 9.4961163736712506989E1; -static const double lntiny = -9.5654310917272452386E1; -static const double invln2 = 1.4426950408889634148E0; -static const double p1 = 1.6666666666666602251E-1; -static const double p2 = -2.7777777777015591216E-3; -static const double p3 = 6.6137563214379341918E-5; -static const double p4 = -1.6533902205465250480E-6; -static const double p5 = 4.1381367970572387085E-8; - -double -exp(double x) -{ - double z, hi, lo, c; - int k; - - if (isnan(x)) - return (x); - - if( x <= lnhuge ) { - if( x >= lntiny ) { - - /* argument reduction : x --> x - k*ln2 */ - - k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */ - - /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */ - - hi=x-k*ln2hi; - x=hi-(lo=k*ln2lo); - - /* return 2^k*[1+x+x*c/(2+c)] */ - z=x*x; - c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); - return scalbn(1.0+(hi-(lo-(x*c)/(2.0-c))),k); - - } - /* end of x > lntiny */ - - else - /* exp(-big#) underflows to zero */ - if(isfinite(x)) return(scalbn(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( isfinite(x) ? scalbn(1.0,5000) : x); -} -DEF_STD(exp); - -/* 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 (isnan(x)) - 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 scalbn(1.+(hi-(lo - c)), k); - } - /* end of x > lntiny */ - - else - /* exp(-big#) underflows to zero */ - if(isfinite(x)) return(scalbn(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( isfinite(x) ? scalbn(1.0,5000) : x); -} diff --git a/lib/libm/noieee_src/n_exp__E.c b/lib/libm/noieee_src/n_exp__E.c deleted file mode 100644 index ccdda138c45..00000000000 --- a/lib/libm/noieee_src/n_exp__E.c +++ /dev/null @@ -1,124 +0,0 @@ -/* $OpenBSD: n_exp__E.c,v 1.12 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_exp__E.c,v 1.1 1995/10/10 23:36:45 ragge 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. 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__E(x,c) - * ASSUMPTION: c << x SO THAT fl(x+c)=x. - * (c is the correction term for x) - * exp__E RETURNS - * - * / exp(x+c) - 1 - x , 1E-19 < |x| < .3465736 - * exp__E(x,c) = | - * \ 0 , |x| < 1E-19. - * - * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) - * KERNEL FUNCTION OF EXP, EXPM1, POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/31/85; - * REVISED BY K.C. NG on 3/16/85, 4/16/85. - * - * Required system supported function: - * copysign(x,y) - * - * Method: - * 1. Rational approximation. Let r=x+c. - * Based on - * 2 * sinh(r/2) - * exp(r) - 1 = ---------------------- , - * cosh(r/2) - sinh(r/2) - * exp__E(r) is computed using - * x*x (x/2)*W - ( Q - ( 2*P + x*P ) ) - * --- + (c + x*[---------------------------------- + c ]) - * 2 1 - W - * where P := p1*x^2 + p2*x^4, - * Q := q1*x^2 + q2*x^4 (for 56 bits precision, add q3*x^6) - * W := x/2-(Q-x*P), - * - * (See the listing below for the values of p1,p2,q1,q2,q3. The poly- - * nomials P and Q may be regarded as the approximations to sinh - * and cosh : - * sinh(r/2) = r/2 + r * P , cosh(r/2) = 1 + Q . ) - * - * The coefficients were obtained by a special Remes algorithm. - * - * Approximation error: - * - * | exp(x) - 1 | 2**(-57), (IEEE double) - * | ------------ - (exp__E(x,0)+x)/x | <= - * | x | 2**(-69). (VAX D) - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double p1 = 1.5150724356786683059E-2; -static const double p2 = 6.3112487873718332688E-5; -static const double q1 = 1.1363478204690669916E-1; -static const double q2 = 1.2624568129896839182E-3; -static const double q3 = 1.5021856115869022674E-6; - -double -__exp__E(double x, double c) -{ - static const double zero=0.0, one=1.0, half=1.0/2.0, small=1.0E-19; - double z,p,q,xp,xh,w; - if(copysign(x,one)>small) { - z = x*x ; - p = z*( p1 +z* p2 ); -#if defined(__vax__) - q = z*( q1 +z*( q2 +z* q3 )); -#else /* defined(__vax__) */ - q = z*( q1 +z* q2 ); -#endif /* defined(__vax__) */ - xp= x*p ; - xh= x*half ; - w = xh-(q-xp) ; - p = p+p; - c += x*((xh*w-(q-(p+xp)))/(one-w)+c); - return(z*half+c); - } - /* end of |x| > small */ - - else { - if(x != zero) { - if (one + small >= 1.0) /* raise the inexact flag */ - return(copysign(zero,x)); - } - else - return(copysign(zero,x)); - } - - /* NOTREACHED */ -} diff --git a/lib/libm/noieee_src/n_expm1.c b/lib/libm/noieee_src/n_expm1.c deleted file mode 100644 index 9494c41f3e4..00000000000 --- a/lib/libm/noieee_src/n_expm1.c +++ /dev/null @@ -1,149 +0,0 @@ -/* $OpenBSD: n_expm1.c,v 1.13 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_expm1.c,v 1.1 1995/10/10 23:36:46 ragge 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. 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. - */ - -/* EXPM1(X) - * RETURN THE EXPONENTIAL OF X MINUS ONE - * 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, 3/7/85, 3/21/85, 4/16/85. - * - * Required system supported functions: - * scalbn(x,n) - * copysign(x,y) - * isfinite(x) - * - * Kernel function: - * exp__E(x,c) - * - * 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 EXPM1(r)=exp(r)-1 by - * - * EXPM1(r=z+c) := z + exp__E(z,c) - * - * 3. EXPM1(x) = 2^k * ( EXPM1(r) + 1-2^-k ). - * - * Remarks: - * 1. When k=1 and z < -0.25, we use the following formula for - * better accuracy: - * EXPM1(x) = 2 * ( (z+0.5) + exp__E(z,c) ) - * 2. To avoid rounding error in 1-2^-k where k is large, we use - * EXPM1(x) = 2^k * { [z+(exp__E(z,c)-2^-k )] + 1 } - * when k>56. - * - * Special cases: - * EXPM1(INF) is INF, EXPM1(NaN) is NaN; - * EXPM1(-INF)= -1; - * for finite argument, only EXPM1(0)=0 is exact. - * - * Accuracy: - * EXPM1(x) returns the exact (exp(x)-1) nearly rounded. In a test run with - * 1,166,000 random arguments on a VAX, the maximum observed error was - * .872 ulps (units of the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double ln2hi = 6.9314718055829871446E-1; -static const double ln2lo = 1.6465949582897081279E-12; -static const double lnhuge = 9.4961163736712506989E1; -static const double invln2 = 1.4426950408889634148E0; - -double -expm1(double x) -{ - static const double one=1.0, half=1.0/2.0, tiny=1e-300; - double z,hi,lo,c; - int k; -#if defined(__vax__) - static prec=56; -#else /* defined(__vax__) */ - static prec=53; -#endif /* defined(__vax__) */ - - if (isnan(x)) - return (x); - - if( x <= lnhuge ) { - if( x >= -40.0 ) { - - /* argument reduction : x - k*ln2 */ - k= invln2 *x+copysign(0.5,x); /* k=NINT(x/ln2) */ - hi=x-k*ln2hi ; - z=hi-(lo=k*ln2lo); - c=(hi-z)-lo; - - if(k==0) return(z+__exp__E(z,c)); - if(k==1) - if(z< -0.25) - {x=z+half;x +=__exp__E(z,c); return(x+x);} - else - {z+=__exp__E(z,c); x=half+z; return(x+x);} - /* end of k=1 */ - - else { - if(k<=prec) - { x=one-scalbn(one,-k); z += __exp__E(z,c);} - else if(k<100) - { x = __exp__E(z,c)-scalbn(one,-k); x+=z; z=one;} - else - { x = __exp__E(z,c)+z; z=one;} - - return (scalbn(x+z,k)); - } - } - /* end of x > lnunfl */ - - else - /* expm1(-big#) rounded to -1 (inexact) */ - if(isfinite(x)) - return(tiny-one); - - /* expm1(-INF) is -1 */ - else return(-one); - } - /* end of x < lnhuge */ - - else - /* expm1(INF) is INF, expm1(+big#) overflows to INF */ - return( isfinite(x) ? scalbn(one,5000) : x); -} diff --git a/lib/libm/noieee_src/n_floor.c b/lib/libm/noieee_src/n_floor.c deleted file mode 100644 index 6c3888fee02..00000000000 --- a/lib/libm/noieee_src/n_floor.c +++ /dev/null @@ -1,174 +0,0 @@ -/* $OpenBSD: n_floor.c,v 1.21 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_floor.c,v 1.1 1995/10/10 23:36:48 ragge 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. 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. - */ - -#include - -#include "mathimpl.h" - -static const double L = 36028797018963968.0E0; /* 2**55 */ -static const float F = 8388608E0f; /* 2**23 */ - -/* - * floor(x) := the largest integer no larger than x; - * ceil(x) := -floor(-x), for all real x. - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -floor(double x) -{ - volatile double y; - - if (isnan(x) || x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -ceil(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x < y ? y-(double)1 : y; - } -} -DEF_STD(floor); -LDBL_CLONE(floor); - -double -ceil(double x) -{ - volatile double y; - - if (isnan(x) || x >= L) /* already an even integer */ - return x; - else if (x < (double)0) - return -floor(-x); - else { /* now 0 <= x < L */ - y = L+x; /* destructive store must be forced */ - y -= L; /* an integer, and |x-y| < 1 */ - return x > y ? y+(double)1 : y; - } -} -DEF_STD(ceil); -LDBL_UNUSED_CLONE(ceil); - -float -floorf(float x) -{ - volatile float y; - - if (isnan(x) || x >= F) /* already an even integer */ - return x; - else if (x < (float)0) - return -ceilf(-x); - else { /* now 0 <= x < F */ - y = F+x; /* destructive store must be forced */ - y -= F; /* an integer, and |x-y| < 1 */ - return x < y ? y-(float)1 : y; - } -} -DEF_STD(floorf); - -float -ceilf(float x) -{ - volatile float y; - - if (isnan(x) || x >= F) /* already an even integer */ - return x; - else if (x < (float)0) - return -floorf(-x); - else { /* now 0 <= x < F */ - y = F+x; /* destructive store must be forced */ - y -= F; /* an integer, and |x-y| < 1 */ - return x > y ? y+(float)1 : y; - } -} -DEF_STD(ceilf); - -/* - * algorithm for rint(x) in pseudo-pascal form ... - * - * real rint(x): real x; - * ... delivers integer nearest x in direction of prevailing rounding - * ... mode - * const L = (last consecutive integer)/2 - * = 2**55; for VAX D - * = 2**52; for IEEE 754 Double - * real s,t; - * begin - * if isnan(x) then return x; ... NaN - * if |x| >= L then return x; ... already an integer - * s := copysign(L,x); - * t := x + s; ... = (x+s) rounded to integer - * return t - s - * end; - * - * Note: Inexact will be signaled if x is not an integer, as is - * customary for IEEE 754. No other signal can be emitted. - */ -double -rint(double x) -{ - double s; - volatile double t; - const double one = 1.0; - - if (isnan(x)) - return (x); - - if (copysign(x, one) >= L) /* already an integer */ - return (x); - - s = copysign(L,x); - t = x + s; /* x+s rounded to integer */ - return (t - s); -} -DEF_STD(rint); -LDBL_CLONE(rint); - -float -rintf(float x) -{ - float s; - volatile float t; - const float one = 1.0f; - - if (isnan(x)) - return (x); - - if (copysignf(x, one) >= F) /* already an integer */ - return (x); - - s = copysignf(F,x); - t = x + s; /* x+s rounded to integer */ - return (t - s); -} -DEF_STD(rintf); diff --git a/lib/libm/noieee_src/n_fmod.c b/lib/libm/noieee_src/n_fmod.c deleted file mode 100644 index f43d15dce0c..00000000000 --- a/lib/libm/noieee_src/n_fmod.c +++ /dev/null @@ -1,80 +0,0 @@ -/* $OpenBSD: n_fmod.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_fmod.c,v 1.1 1995/10/10 23:36:49 ragge Exp $ */ -/* - * Copyright (c) 1989, 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. 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. - */ - -#include "math.h" -#include "mathimpl.h" - -/* fmod.c - * - * SYNOPSIS - * - * #include - * double fmod(double x, double y) - * - * DESCRIPTION - * - * The fmod function computes the floating-point remainder of x/y. - * - * RETURNS - * - * The fmod function returns the value x-i*y, for some integer i - * such that, if y is nonzero, the result has the same sign as x and - * magnitude less than the magnitude of y. - * - * On a VAX or CCI, - * - * fmod(x,0) traps/faults on floating-point divided-by-zero. - * - * On IEEE-754 conforming machines with "isnan()" primitive, - * - * fmod(x,0), fmod(INF,y) are invalid operations and NaN is returned. - * - */ - -double -fmod(double x, double y) -{ - int ir,iy; - double r,w; - - if (y == (double)0 || isnan(y) || !isfinite(x)) - return (x*y)/(x*y); - - r = fabs(x); - y = fabs(y); - (void)frexp(y,&iy); - while (r >= y) { - (void)frexp(r,&ir); - w = ldexp(y,ir-iy); - r -= w <= r ? w : w*(double)0.5; - } - return x >= (double)0 ? r : -r; -} diff --git a/lib/libm/noieee_src/n_hypot.c b/lib/libm/noieee_src/n_hypot.c deleted file mode 100644 index 6653e956844..00000000000 --- a/lib/libm/noieee_src/n_hypot.c +++ /dev/null @@ -1,195 +0,0 @@ -/* $OpenBSD: n_hypot.c,v 1.6 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_cabs.c,v 1.1 1995/10/10 23:36:39 ragge 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. 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. - */ - -/* HYPOT(X,Y) - * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 11/28/84; - * REVISED BY K.C. NG, 7/12/85. - * - * Required system supported functions : - * copysign(x,y) - * isfinite(x) - * scalbn(x,N) - * sqrt(x) - * - * Method : - * 1. replace x by |x| and y by |y|, and swap x and - * y if y > x (hence x is never smaller than y). - * 2. Hypot(x,y) is computed by: - * Case I, x/y > 2 - * - * y - * hypot = x + ----------------------------- - * 2 - * sqrt ( 1 + [x/y] ) + x/y - * - * Case II, x/y <= 2 - * y - * hypot = x + -------------------------------------------------- - * 2 - * [x/y] - 2 - * (sqrt(2)+1) + (x-y)/y + ----------------------------- - * 2 - * sqrt ( 1 + [x/y] ) + sqrt(2) - * - * - * - * Special cases: - * hypot(x,y) is INF if x or y is +INF or -INF; else - * hypot(x,y) is NAN if x or y is NAN. - * - * Accuracy: - * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units - * in the last place). See Kahan's "Interval Arithmetic Options in the - * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics - * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate - * code follows in comments.) In a test run with 500,000 random arguments - * on a VAX, the maximum observed error was .959 ulps. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double r2p1hi = 2.4142135623730950345E0; -static const double r2p1lo = 1.4349369327986523769E-17; -static const double sqrt2 = 1.4142135623730950622E0; - -float -hypotf(float x, float y) -{ - return (float)hypot((double) x, (double) y); -} - -double -hypot(double x, double y) -{ - static const double zero=0, one=1, - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double t,r; - int exp; - - if(isfinite(x)) - if(isfinite(y)) - { - x=copysign(x,one); - y=copysign(y,one); - if(y > x) - { t=x; x=y; y=t; } - if(x == zero) return(zero); - if(y == zero) return(x); - exp= logb(x); - if (exp - (int)logb(y) > ibig) { - if (one + small >= 1.0) /* raise inexact flag */ - return(x); /* return |x| */ - } - - /* start computing sqrt(x^2 + y^2) */ - r=x-y; - if(r>y) { /* x/y > 2 */ - r=x/y; - r=r+sqrt(one+r*r); } - else { /* 1 <= x/y <= 2 */ - r/=y; t=r*(r+2.0); - r+=t/(sqrt2+sqrt(2.0+t)); - r+=r2p1lo; r+=r2p1hi; } - - r=y/r; - return(x+r); - - } - - else if(isinf(y)) /* y is +-INF */ - return(copysign(y,one)); - else - return(y); /* y is NaN and x is finite */ - - else if(isinf(x)) /* x is +-INF */ - return (copysign(x,one)); - else if(isfinite(y)) - return(x); /* x is NaN, y is finite */ - else if (isnan(y)) - return (y); - else return(copysign(y,one)); /* y is INF */ -} -DEF_STD(hypot); -LDBL_CLONE(hypot); - -/* A faster but less accurate version of cabs(x,y) */ -#if 0 -double -hypot(double x, double y) -{ - static const double zero=0, one=1; - small=1.0E-18; /* fl(1+small)==1 */ - static const ibig=30; /* fl(1+2**(2*ibig))==1 */ - double temp; - int exp; - - if(isfinite(x)) - if(isfinite(y)) - { - x=copysign(x,one); - y=copysign(y,one); - if(y > x) - { temp=x; x=y; y=temp; } - if(x == zero) return(zero); - if(y == zero) return(x); - exp= logb(x); - x=scalbn(x,-exp); - if (exp - (int)logb(y) > ibig) { - if (one + small >= 1.0) /* raise inexact flag */ - return(scalbn(x,exp)); /* return |x| */ - } - else y=scalbn(y,-exp); - return(scalbn(sqrt(x*x+y*y),exp)); - } - - else if(isinf(y)) /* y is +-INF */ - return(copysign(y,one)); - else - return(y); /* y is NaN and x is finite */ - - else if(isinf(x)) /* x is +-INF */ - return (copysign(x,one)); - else if(isfinite(y)) - return(x); /* x is NaN, y is finite */ - else if(isnan(y)) return(y); /* x and y is NaN */ - else return(copysign(y,one)); /* y is INF */ -} -#endif diff --git a/lib/libm/noieee_src/n_j0.c b/lib/libm/noieee_src/n_j0.c deleted file mode 100644 index 6274b17cde2..00000000000 --- a/lib/libm/noieee_src/n_j0.c +++ /dev/null @@ -1,431 +0,0 @@ -/* $OpenBSD: n_j0.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_j0.c,v 1.1 1995/10/10 23:36:52 ragge 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. 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. - */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j0(double x), y0(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j0(x): - * 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ... - * 2. Reduce x to |x| since j0(x)=j0(-x), and - * for x in (0,2) - * j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x; - * (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 ) - * for x in (2,inf) - * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * as follow: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (cos(x) + sin(x)) - * sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j0(nan)= nan - * j0(0) = 1 - * j0(inf) = 0 - * - * Method -- y0(x): - * 1. For x<2. - * Since - * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) - * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. - * We use the following function to approximate y0, - * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - * where - * U(z) = u0 + u1*z + ... + u6*z^6 - * V(z) = 1 + v1*z + ... + v4*z^4 - * with absolute approximation error bounded by 2**-72. - * Note: For tiny x, U/V = u0 and j0(x)~1, hence - * y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27) - * 2. For x>=2. - * y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0)) - * where x0 = x-pi/4. It is better to compute sin(x0),cos(x0) - * by the method mentioned above. - * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. - */ - -#include -#include -#include - -#if defined(__vax__) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pzero(double), qzero(double); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - /* R0/S0 on [0, 2.00] */ -r02 = 1.562499999999999408594634421055018003102e-0002, -r03 = -1.899792942388547334476601771991800712355e-0004, -r04 = 1.829540495327006565964161150603950916854e-0006, -r05 = -4.618326885321032060803075217804816988758e-0009, -s01 = 1.561910294648900170180789369288114642057e-0002, -s02 = 1.169267846633374484918570613449245536323e-0004, -s03 = 5.135465502073181376284426245689510134134e-0007, -s04 = 1.166140033337900097836930825478674320464e-0009; - -double -j0(double x) -{ - double z, s,c,ss,cc,r,u,v; - - if (!isfinite(x)) - if (_IEEE) return one/(x*x); - else return (0); - x = fabs(x); - if (x >= 2.0) { /* |x| >= 2.0 */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c) 6.80564733841876927e+38) /* 2^129 */ - z = (invsqrtpi*cc)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*cc-v*ss)/sqrt(x); - } - return z; - } - if (x < 1.220703125e-004) { /* |x| < 2**-13 */ - if (huge+x > one) { /* raise inexact if x != 0 */ - if (x < 7.450580596923828125e-009) /* |x|<2**-27 */ - return one; - else return (one - 0.25*x*x); - } - } - z = x*x; - r = z*(r02+z*(r03+z*(r04+z*r05))); - s = one+z*(s01+z*(s02+z*(s03+z*s04))); - if (x < one) { /* |x| < 1.00 */ - return (one + z*(-0.25+(r/s))); - } else { - u = 0.5*x; - return ((one+u)*(one-u)+z*(r/s)); - } -} - -static double -u00 = -7.380429510868722527422411862872999615628e-0002, -u01 = 1.766664525091811069896442906220827182707e-0001, -u02 = -1.381856719455968955440002438182885835344e-0002, -u03 = 3.474534320936836562092566861515617053954e-0004, -u04 = -3.814070537243641752631729276103284491172e-0006, -u05 = 1.955901370350229170025509706510038090009e-0008, -u06 = -3.982051941321034108350630097330144576337e-0011, -v01 = 1.273048348341237002944554656529224780561e-0002, -v02 = 7.600686273503532807462101309675806839635e-0005, -v03 = 2.591508518404578033173189144579208685163e-0007, -v04 = 4.411103113326754838596529339004302243157e-0010; - -double -y0(double x) -{ - double z, s, c, ss, cc, u, v; - /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if (!isfinite(x)) - if (_IEEE) - return (one/(x+x*x)); - else - return (0); - if (x == 0) - if (_IEEE) return (-one/zero); - else return(infnan(-ERANGE)); - if (x<0) - if (_IEEE) return (zero/zero); - else return (infnan(EDOM)); - if (x >= 2.00) { /* |x| >= 2.0 */ - /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) - * where x0 = x-pi/4 - * Better formula: - * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) - * = 1/sqrt(2) * (sin(x) + cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - s = sin(x); - c = cos(x); - ss = s-c; - cc = s+c; - /* - * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) - * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) - */ - if (x < .5 * DBL_MAX) { /* make sure x+x not overflow */ - z = -cos(x+x); - if ((s*c) 6.80564733841876927e+38) /* > 2^129 */ - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pzero(x); v = qzero(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= 7.450580596923828125e-009) { /* x < 2**-27 */ - return (u00 + tpi*log(x)); - } - z = x*x; - u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); - v = one+z*(v01+z*(v02+z*(v03+z*v04))); - return (u/v + tpi*(j0(x)*log(x))); -} - -/* The asymptotic expansions of pzero is - * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. - * For x >= 2, We approximate pzero by - * pzero(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pzero(x)-1-R/S | <= 2 ** ( -60.26) - */ -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -7.031249999999003994151563066182798210142e-0002, - -8.081670412753498508883963849859423939871e+0000, - -2.570631056797048755890526455854482662510e+0002, - -2.485216410094288379417154382189125598962e+0003, - -5.253043804907295692946647153614119665649e+0003, -}; -static double ps8[5] = { - 1.165343646196681758075176077627332052048e+0002, - 3.833744753641218451213253490882686307027e+0003, - 4.059785726484725470626341023967186966531e+0004, - 1.167529725643759169416844015694440325519e+0005, - 4.762772841467309430100106254805711722972e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -1.141254646918944974922813501362824060117e-0011, - -7.031249408735992804117367183001996028304e-0002, - -4.159610644705877925119684455252125760478e+0000, - -6.767476522651671942610538094335912346253e+0001, - -3.312312996491729755731871867397057689078e+0002, - -3.464333883656048910814187305901796723256e+0002, -}; -static double ps5[5] = { - 6.075393826923003305967637195319271932944e+0001, - 1.051252305957045869801410979087427910437e+0003, - 5.978970943338558182743915287887408780344e+0003, - 9.625445143577745335793221135208591603029e+0003, - 2.406058159229391070820491174867406875471e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -2.547046017719519317420607587742992297519e-0009, - -7.031196163814817199050629727406231152464e-0002, - -2.409032215495295917537157371488126555072e+0000, - -2.196597747348830936268718293366935843223e+0001, - -5.807917047017375458527187341817239891940e+0001, - -3.144794705948885090518775074177485744176e+0001, -}; -static double ps3[5] = { - 3.585603380552097167919946472266854507059e+0001, - 3.615139830503038919981567245265266294189e+0002, - 1.193607837921115243628631691509851364715e+0003, - 1.127996798569074250675414186814529958010e+0003, - 1.735809308133357510239737333055228118910e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -8.875343330325263874525704514800809730145e-0008, - -7.030309954836247756556445443331044338352e-0002, - -1.450738467809529910662233622603401167409e+0000, - -7.635696138235277739186371273434739292491e+0000, - -1.119316688603567398846655082201614524650e+0001, - -3.233645793513353260006821113608134669030e+0000, -}; -static double ps2[5] = { - 2.222029975320888079364901247548798910952e+0001, - 1.362067942182152109590340823043813120940e+0002, - 2.704702786580835044524562897256790293238e+0002, - 1.538753942083203315263554770476850028583e+0002, - 1.465761769482561965099880599279699314477e+0001, -}; - -static double pzero(double x) -{ - double *p,*q,z,r,s; - if (x >= 8.00) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else if (x >= 2.00) {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return one+ r/s; -} - - -/* For x >= 8, the asymptotic expansions of qzero is - * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. - * We approximate pzero by - * qzero(x) = s*(-1.25 + (R/S)) - * where R = qr0 + qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs0*s^2 + ... + qs5*s^12 - * and - * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) - */ -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 7.324218749999350414479738504551775297096e-0002, - 1.176820646822526933903301695932765232456e+0001, - 5.576733802564018422407734683549251364365e+0002, - 8.859197207564685717547076568608235802317e+0003, - 3.701462677768878501173055581933725704809e+0004, -}; -static double qs8[6] = { - 1.637760268956898345680262381842235272369e+0002, - 8.098344946564498460163123708054674227492e+0003, - 1.425382914191204905277585267143216379136e+0005, - 8.033092571195144136565231198526081387047e+0005, - 8.405015798190605130722042369969184811488e+0005, - -3.438992935378666373204500729736454421006e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.840859635945155400568380711372759921179e-0011, - 7.324217666126847411304688081129741939255e-0002, - 5.835635089620569401157245917610984757296e+0000, - 1.351115772864498375785526599119895942361e+0002, - 1.027243765961641042977177679021711341529e+0003, - 1.989977858646053872589042328678602481924e+0003, -}; -static double qs5[6] = { - 8.277661022365377058749454444343415524509e+0001, - 2.077814164213929827140178285401017305309e+0003, - 1.884728877857180787101956800212453218179e+0004, - 5.675111228949473657576693406600265778689e+0004, - 3.597675384251145011342454247417399490174e+0004, - -5.354342756019447546671440667961399442388e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 4.377410140897386263955149197672576223054e-0009, - 7.324111800429115152536250525131924283018e-0002, - 3.344231375161707158666412987337679317358e+0000, - 4.262184407454126175974453269277100206290e+0001, - 1.708080913405656078640701512007621675724e+0002, - 1.667339486966511691019925923456050558293e+0002, -}; -static double qs3[6] = { - 4.875887297245871932865584382810260676713e+0001, - 7.096892210566060535416958362640184894280e+0002, - 3.704148226201113687434290319905207398682e+0003, - 6.460425167525689088321109036469797462086e+0003, - 2.516333689203689683999196167394889715078e+0003, - -1.492474518361563818275130131510339371048e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.504444448869832780257436041633206366087e-0007, - 7.322342659630792930894554535717104926902e-0002, - 1.998191740938159956838594407540292600331e+0000, - 1.449560293478857407645853071687125850962e+0001, - 3.166623175047815297062638132537957315395e+0001, - 1.625270757109292688799540258329430963726e+0001, -}; -static double qs2[6] = { - 3.036558483552191922522729838478169383969e+0001, - 2.693481186080498724211751445725708524507e+0002, - 8.447837575953201460013136756723746023736e+0002, - 8.829358451124885811233995083187666981299e+0002, - 2.126663885117988324180482985363624996652e+0002, - -5.310954938826669402431816125780738924463e+0000, -}; - -static double qzero(double x) -{ - double *p,*q, s,r,z; - if (x >= 8.00) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else if (x >= 2.00) {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (-.125 + r/s)/x; -} diff --git a/lib/libm/noieee_src/n_j1.c b/lib/libm/noieee_src/n_j1.c deleted file mode 100644 index 0bb2181b68b..00000000000 --- a/lib/libm/noieee_src/n_j1.c +++ /dev/null @@ -1,439 +0,0 @@ -/* $OpenBSD: n_j1.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_j1.c,v 1.1 1995/10/10 23:36:53 ragge 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. 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. - */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* double j1(double x), y1(double x) - * Bessel function of the first and second kinds of order zero. - * Method -- j1(x): - * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... - * 2. Reduce x to |x| since j1(x)=-j1(-x), and - * for x in (0,2) - * j1(x) = x/2 + x*z*R0/S0, where z = x*x; - * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) - * for x in (2,inf) - * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * as follows: - * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (sin(x) + cos(x)) - * (To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one.) - * - * 3 Special cases - * j1(nan)= nan - * j1(0) = 0 - * j1(inf) = 0 - * - * Method -- y1(x): - * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN - * 2. For x<2. - * Since - * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) - * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. - * We use the following function to approximate y1, - * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 - * where for x in [0,2] (abs err less than 2**-65.89) - * U(z) = u0 + u1*z + ... + u4*z^4 - * V(z) = 1 + v1*z + ... + v5*z^5 - * Note: For tiny x, 1/x dominate y1 and hence - * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) - * 3. For x>=2. - * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) - * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) - * by method mentioned above. - */ - -#include -#include -#include - -#if defined(__vax__) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double pone(), qone(); - -static double -huge = 1e300, -zero = 0.0, -one = 1.0, -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -tpi = 0.636619772367581343075535053490057448, - - /* R0/S0 on [0,2] */ -r00 = -6.250000000000000020842322918309200910191e-0002, -r01 = 1.407056669551897148204830386691427791200e-0003, -r02 = -1.599556310840356073980727783817809847071e-0005, -r03 = 4.967279996095844750387702652791615403527e-0008, -s01 = 1.915375995383634614394860200531091839635e-0002, -s02 = 1.859467855886309024045655476348872850396e-0004, -s03 = 1.177184640426236767593432585906758230822e-0006, -s04 = 5.046362570762170559046714468225101016915e-0009, -s05 = 1.235422744261379203512624973117299248281e-0011; - -#define two_129 6.80564733841876926e+038 /* 2^129 */ -#define two_m54 5.55111512312578270e-017 /* 2^-54 */ - -double -j1(double x) -{ - double z, s,c,ss,cc,r,u,v,y; - y = fabs(x); - if (!isfinite(x)) /* Inf or NaN */ - if (isnan(x)) - return(x); - else - return (copysign(x, zero)); - y = fabs(x); - if (y >= 2) /* |x| >= 2.0 */ - { - s = sin(y); - c = cos(y); - ss = -s-c; - cc = s-c; - if (y < .5*DBL_MAX) { /* make sure y+y not overflow */ - z = cos(y+y); - if ((s*c) two_129) /* x > 2^129 */ - z = (invsqrtpi*cc)/sqrt(y); - else -#endif /* !defined(__vax__) */ - { - u = pone(y); v = qone(y); - z = invsqrtpi*(u*cc-v*ss)/sqrt(y); - } - if (x < 0) return -z; - else return z; - } - if (y < 7.450580596923828125e-009) { /* |x|<2**-27 */ - if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ - } - z = x*x; - r = z*(r00+z*(r01+z*(r02+z*r03))); - s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); - r *= x; - return (x*0.5+r/s); -} - -static double u0[5] = { - -1.960570906462389484206891092512047539632e-0001, - 5.044387166398112572026169863174882070274e-0002, - -1.912568958757635383926261729464141209569e-0003, - 2.352526005616105109577368905595045204577e-0005, - -9.190991580398788465315411784276789663849e-0008, -}; -static double v0[5] = { - 1.991673182366499064031901734535479833387e-0002, - 2.025525810251351806268483867032781294682e-0004, - 1.356088010975162198085369545564475416398e-0006, - 6.227414523646214811803898435084697863445e-0009, - 1.665592462079920695971450872592458916421e-0011, -}; - -double -y1(double x) -{ - double z, s, c, ss, cc, u, v; - /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if (!isfinite(x)) - if (x < 0) - return(zero/zero); - else if (x > 0) - return (0); - else - return(x); - if (x <= 0) { - if (_IEEE && x == 0) return -one/zero; - else if(x == 0) return(infnan(-ERANGE)); - else if(_IEEE) return (zero/zero); - else return(infnan(EDOM)); - } - if (x >= 2) /* |x| >= 2.0 */ - { - s = sin(x); - c = cos(x); - ss = -s-c; - cc = s-c; - if (x < .5 * DBL_MAX) /* make sure x+x not overflow */ - { - z = cos(x+x); - if ((s*c)>zero) cc = z/ss; - else ss = z/cc; - } - /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) - * where x0 = x-3pi/4 - * Better formula: - * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) - * = 1/sqrt(2) * (sin(x) - cos(x)) - * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) - * = -1/sqrt(2) * (cos(x) + sin(x)) - * To avoid cancellation, use - * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - * to compute the worse one. - */ - if (_IEEE && x>two_129) - z = (invsqrtpi*ss)/sqrt(x); - else { - u = pone(x); v = qone(x); - z = invsqrtpi*(u*ss+v*cc)/sqrt(x); - } - return z; - } - if (x <= two_m54) { /* x < 2**-54 */ - return (-tpi/x); - } - z = x*x; - u = u0[0]+z*(u0[1]+z*(u0[2]+z*(u0[3]+z*u0[4]))); - v = one+z*(v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])))); - return (x*(u/v) + tpi*(j1(x)*log(x)-one/x)); -} - -/* For x >= 8, the asymptotic expansions of pone is - * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. - * We approximate pone by - * pone(x) = 1 + (R/S) - * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 - * S = 1 + ps0*s^2 + ... + ps4*s^10 - * and - * | pone(x)-1-R/S | <= 2 ** ( -60.06) - */ - -static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - 1.171874999999886486643746274751925399540e-0001, - 1.323948065930735690925827997575471527252e+0001, - 4.120518543073785433325860184116512799375e+0002, - 3.874745389139605254931106878336700275601e+0003, - 7.914479540318917214253998253147871806507e+0003, -}; -static double ps8[5] = { - 1.142073703756784104235066368252692471887e+0002, - 3.650930834208534511135396060708677099382e+0003, - 3.695620602690334708579444954937638371808e+0004, - 9.760279359349508334916300080109196824151e+0004, - 3.080427206278887984185421142572315054499e+0004, -}; - -static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - 1.319905195562435287967533851581013807103e-0011, - 1.171874931906140985709584817065144884218e-0001, - 6.802751278684328781830052995333841452280e+0000, - 1.083081829901891089952869437126160568246e+0002, - 5.176361395331997166796512844100442096318e+0002, - 5.287152013633375676874794230748055786553e+0002, -}; -static double ps5[5] = { - 5.928059872211313557747989128353699746120e+0001, - 9.914014187336144114070148769222018425781e+0002, - 5.353266952914879348427003712029704477451e+0003, - 7.844690317495512717451367787640014588422e+0003, - 1.504046888103610723953792002716816255382e+0003, -}; - -static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - 3.025039161373736032825049903408701962756e-0009, - 1.171868655672535980750284752227495879921e-0001, - 3.932977500333156527232725812363183251138e+0000, - 3.511940355916369600741054592597098912682e+0001, - 9.105501107507812029367749771053045219094e+0001, - 4.855906851973649494139275085628195457113e+0001, -}; -static double ps3[5] = { - 3.479130950012515114598605916318694946754e+0001, - 3.367624587478257581844639171605788622549e+0002, - 1.046871399757751279180649307467612538415e+0003, - 8.908113463982564638443204408234739237639e+0002, - 1.037879324396392739952487012284401031859e+0002, -}; - -static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - 1.077108301068737449490056513753865482831e-0007, - 1.171762194626833490512746348050035171545e-0001, - 2.368514966676087902251125130227221462134e+0000, - 1.224261091482612280835153832574115951447e+0001, - 1.769397112716877301904532320376586509782e+0001, - 5.073523125888185399030700509321145995160e+0000, -}; -static double ps2[5] = { - 2.143648593638214170243114358933327983793e+0001, - 1.252902271684027493309211410842525120355e+0002, - 2.322764690571628159027850677565128301361e+0002, - 1.176793732871470939654351793502076106651e+0002, - 8.364638933716182492500902115164881195742e+0000, -}; - -static double pone(double x) -{ - double *p,*q,z,r,s; - if (x >= 8.0) {p = pr8; q= ps8;} - else if (x >= 4.54545211791992188) {p = pr5; q= ps5;} - else if (x >= 2.85714149475097656) {p = pr3; q= ps3;} - else /* if (x >= 2.0) */ {p = pr2; q= ps2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); - return (one + r/s); -} - - -/* For x >= 8, the asymptotic expansions of qone is - * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. - * We approximate pone by - * qone(x) = s*(0.375 + (R/S)) - * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 - * S = 1 + qs1*s^2 + ... + qs6*s^12 - * and - * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) - */ - -static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ - 0.0, - -1.025390624999927207385863635575804210817e-0001, - -1.627175345445899724355852152103771510209e+0001, - -7.596017225139501519843072766973047217159e+0002, - -1.184980667024295901645301570813228628541e+0004, - -4.843851242857503225866761992518949647041e+0004, -}; -static double qs8[6] = { - 1.613953697007229231029079421446916397904e+0002, - 7.825385999233484705298782500926834217525e+0003, - 1.338753362872495800748094112937868089032e+0005, - 7.196577236832409151461363171617204036929e+0005, - 6.666012326177764020898162762642290294625e+0005, - -2.944902643038346618211973470809456636830e+0005, -}; - -static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ - -2.089799311417640889742251585097264715678e-0011, - -1.025390502413754195402736294609692303708e-0001, - -8.056448281239359746193011295417408828404e+0000, - -1.836696074748883785606784430098756513222e+0002, - -1.373193760655081612991329358017247355921e+0003, - -2.612444404532156676659706427295870995743e+0003, -}; -static double qs5[6] = { - 8.127655013843357670881559763225310973118e+0001, - 1.991798734604859732508048816860471197220e+0003, - 1.746848519249089131627491835267411777366e+0004, - 4.985142709103522808438758919150738000353e+0004, - 2.794807516389181249227113445299675335543e+0004, - -4.719183547951285076111596613593553911065e+0003, -}; - -static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ - -5.078312264617665927595954813341838734288e-0009, - -1.025378298208370901410560259001035577681e-0001, - -4.610115811394734131557983832055607679242e+0000, - -5.784722165627836421815348508816936196402e+0001, - -2.282445407376317023842545937526967035712e+0002, - -2.192101284789093123936441805496580237676e+0002, -}; -static double qs3[6] = { - 4.766515503237295155392317984171640809318e+0001, - 6.738651126766996691330687210949984203167e+0002, - 3.380152866795263466426219644231687474174e+0003, - 5.547729097207227642358288160210745890345e+0003, - 1.903119193388108072238947732674639066045e+0003, - -1.352011914443073322978097159157678748982e+0002, -}; - -static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ - -1.783817275109588656126772316921194887979e-0007, - -1.025170426079855506812435356168903694433e-0001, - -2.752205682781874520495702498875020485552e+0000, - -1.966361626437037351076756351268110418862e+0001, - -4.232531333728305108194363846333841480336e+0001, - -2.137192117037040574661406572497288723430e+0001, -}; -static double qs2[6] = { - 2.953336290605238495019307530224241335502e+0001, - 2.529815499821905343698811319455305266409e+0002, - 7.575028348686454070022561120722815892346e+0002, - 7.393932053204672479746835719678434981599e+0002, - 1.559490033366661142496448853793707126179e+0002, - -4.959498988226281813825263003231704397158e+0000, -}; - -static double qone(double x) -{ - double *p,*q, s,r,z; - if (x >= 8.0) {p = qr8; q= qs8;} - else if (x >= 4.54545211791992188) {p = qr5; q= qs5;} - else if (x >= 2.85714149475097656) {p = qr3; q= qs3;} - else /* if (x >= 2.0) */ {p = qr2; q= qs2;} - z = one/(x*x); - r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); - s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); - return (.375 + r/s)/x; -} diff --git a/lib/libm/noieee_src/n_jn.c b/lib/libm/noieee_src/n_jn.c deleted file mode 100644 index d8c6111f54f..00000000000 --- a/lib/libm/noieee_src/n_jn.c +++ /dev/null @@ -1,306 +0,0 @@ -/* $OpenBSD: n_jn.c,v 1.8 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_jn.c,v 1.1 1995/10/10 23:36:54 ragge 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. 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. - */ - -/* - * 16 December 1992 - * Minor modifications by Peter McIlroy to adapt non-IEEE architecture. - */ - -/* - * ==================================================== - * Copyright (C) 1992 by Sun Microsystems, Inc. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - * ******************* WARNING ******************** - * This is an alpha version of SunPro's FDLIBM (Freely - * Distributable Math Library) for IEEE double precision - * arithmetic. FDLIBM is a basic math library written - * in C that runs on machines that conform to IEEE - * Standard 754/854. This alpha version is distributed - * for testing purpose. Those who use this software - * should report any bugs to - * - * fdlibm-comments@sunpro.eng.sun.com - * - * -- K.C. Ng, Oct 12, 1992 - * ************************************************ - */ - -/* - * jn(int n, double x), yn(int n, double x) - * floating point Bessel's function of the 1st and 2nd kind - * of order n - * - * Special cases: - * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; - * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. - * Note 2. About jn(n,x), yn(n,x) - * For n=0, j0(x) is called, - * for n=1, j1(x) is called, - * for nx, a continued fraction approximation to - * j(n,x)/j(n-1,x) is evaluated and then backward - * recursion is used starting from a supposed value - * for j(n,x). The resulting value of j(0,x) is - * compared with the actual value to correct the - * supposed value of j(n,x). - * - * yn(n,x) is similar in all respects, except - * that forward recursion is used for all - * values of n>1. - * - */ - -#include -#include -#include - -#if defined(__vax__) -#define _IEEE 0 -#else -#define _IEEE 1 -#define infnan(x) (0.0) -#endif - -static double -invsqrtpi= 5.641895835477562869480794515607725858441e-0001, -two = 2.0, -zero = 0.0, -one = 1.0; - -double -jn(int n, double x) -{ - int i, sgn; - double a, b, temp; - double z, w; - - /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x) - * Thus, J(-n,x) = J(n,-x) - */ - /* if J(n,NaN) is NaN */ - if (_IEEE && isnan(x)) return x+x; - if (n<0){ - n = -n; - x = -x; - } - if (n==0) return(j0(x)); - if (n==1) return(j1(x)); - sgn = (n&1)&(x < zero); /* even n -- 0, odd n -- sign(x) */ - x = fabs(x); - if (x == 0 || !isfinite (x)) /* if x is 0 or inf */ - b = zero; - else if ((double) n <= x) { - /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ - if (_IEEE && x >= 8.148143905337944345e+090) { - /* x >= 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch(n&3) { - case 0: temp = cos(x)+sin(x); break; - case 1: temp = -cos(x)+sin(x); break; - case 2: temp = -cos(x)-sin(x); break; - case 3: temp = cos(x)-sin(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = j0(x); - b = j1(x); - for(i=1;i 33) /* underflow */ - b = zero; - else { - temp = x*0.5; b = temp; - for (a=one,i=2;i<=n;i++) { - a *= (double)i; /* a = n! */ - b *= temp; /* b = (x/2)^n */ - } - b = b/a; - } - } else { - /* use backward recurrence */ - /* x x^2 x^2 - * J(n,x)/J(n-1,x) = ---- ------ ------ ..... - * 2n - 2(n+1) - 2(n+2) - * - * 1 1 1 - * (for large x) = ---- ------ ------ ..... - * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - - * x x x - * - * Let w = 2n/x and h=2/x, then the above quotient - * is equal to the continued fraction: - * 1 - * = ----------------------- - * 1 - * w - ----------------- - * 1 - * w+h - --------- - * w+2h - ... - * - * To determine how many terms needed, let - * Q(0) = w, Q(1) = w(w+h) - 1, - * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple - */ - /* determine k */ - double t,v; - double q0,q1,h,tmp; int k,m; - w = (n+n)/(double)x; h = 2.0/(double)x; - q0 = w; z = w+h; q1 = w*z - 1.0; k=1; - while (q1<1.0e9) { - k += 1; z += h; - tmp = z*q1 - q0; - q0 = q1; - q1 = tmp; - } - m = n+n; - for(t=zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t); - a = t; - b = one; - /* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n) - * Hence, if n*(log(2n/x)) > ... - * single 8.8722839355e+01 - * double 7.09782712893383973096e+02 - * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result will - * likely underflow to zero - */ - tmp = n; - v = two/x; - tmp = tmp*log(fabs(v*tmp)); - for (i=n-1;i>0;i--){ - temp = b; - b = ((i+i)/x)*b - a; - a = temp; - /* scale b to avoid spurious overflow */ -# if defined(__vax__) -# define BMAX 1e13 -# else -# define BMAX 1e100 -# endif /* defined(__vax__) */ - if (b > BMAX) { - a /= b; - t /= b; - b = one; - } - } - b = (t*j0(x)/b); - } - } - return ((sgn == 1) ? -b : b); -} - -double -yn(int n, double x) -{ - int i, sign; - double a, b, temp; - - /* Y(n,NaN), Y(n, x < 0) is NaN */ - if (x <= 0 || isnan(x)) - if (_IEEE && x < 0) return zero/zero; - else if (x < 0) return (infnan(EDOM)); - else if (_IEEE) return -one/zero; - else return(infnan(-ERANGE)); - else if (!isfinite(x)) return(0); - sign = 1; - if (n<0){ - n = -n; - sign = 1 - ((n&1)<<2); - } - if (n == 0) return(y0(x)); - if (n == 1) return(sign*y1(x)); - if(_IEEE && x >= 8.148143905337944345e+090) { /* x > 2**302 */ - /* (x >> n**2) - * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), - * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then - * - * n sin(xn)*sqt2 cos(xn)*sqt2 - * ---------------------------------- - * 0 s-c c+s - * 1 -s-c -c+s - * 2 -s+c -c-s - * 3 s+c c-s - */ - switch (n&3) { - case 0: temp = sin(x)-cos(x); break; - case 1: temp = -sin(x)-cos(x); break; - case 2: temp = -sin(x)+cos(x); break; - case 3: temp = sin(x)+cos(x); break; - } - b = invsqrtpi*temp/sqrt(x); - } else { - a = y0(x); - b = y1(x); - /* quit if b is -inf */ - for (i = 1; i < n && !isfinite(b); i++){ - temp = b; - b = ((double)(i+i)/x)*b - a; - a = temp; - } - } - if (!_IEEE && !isfinite(b)) - return (infnan(-sign * ERANGE)); - return ((sign > 0) ? b : -b); -} diff --git a/lib/libm/noieee_src/n_lgamma.c b/lib/libm/noieee_src/n_lgamma.c deleted file mode 100644 index 44a1922e000..00000000000 --- a/lib/libm/noieee_src/n_lgamma.c +++ /dev/null @@ -1,327 +0,0 @@ -/* $OpenBSD: n_lgamma.c,v 1.10 2016/09/12 04:39:47 guenther 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. 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. - */ - -/* - * Coded by Peter McIlroy, Nov 1992; - * - * The financial support of UUNET Communications Services is greatfully - * acknowledged. - */ - -#include -#include - -#include "mathimpl.h" - -/* Log gamma function. - * Error: x > 0 error < 1.3ulp. - * x > 4, error < 1ulp. - * x > 9, error < .6ulp. - * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0) - * Method: - * x > 6: - * Use the asymptotic expansion (Stirling's Formula) - * 0 < x < 6: - * Use tgamma(x+1) = x*tgamma(x) for argument reduction. - * Use rational approximation in - * the range 1.2, 2.5 - * Two approximations are used, one centered at the - * minimum to ensure monotonicity; one centered at 2 - * to maintain small relative error. - * x < 0: - * Use the reflection formula, - * G(1-x)G(x) = PI/sin(PI*x) - * Special values: - * non-positive integer returns +Inf. - * NaN returns NaN -*/ -#if defined(__vax__) -#define _IEEE 0 -/* double and float have same size exponent field */ -#define TRUNC(x) x = (double) (float) (x) -#else -static int endian; -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -static double small_lgam(double); -static double large_lgam(double); -static double neg_lgam(double); -static const double one = 1.0; -extern int signgam; - -#define UNDERFL (1e-1020 * 1e-1020) - -#define LEFT (1.0 - (x0 + .25)) -#define RIGHT (x0 - .218) -/* - * Constants for approximation in [1.244,1.712] - */ -#define x0 0.461632144968362356785 -#define x0_lo -.000000000000000015522348162858676890521 -#define a0_hi -0.12148629128932952880859 -#define a0_lo .0000000007534799204229502 -#define r0 -2.771227512955130520e-002 -#define r1 -2.980729795228150847e-001 -#define r2 -3.257411333183093394e-001 -#define r3 -1.126814387531706041e-001 -#define r4 -1.129130057170225562e-002 -#define r5 -2.259650588213369095e-005 -#define s0 1.714457160001714442e+000 -#define s1 2.786469504618194648e+000 -#define s2 1.564546365519179805e+000 -#define s3 3.485846389981109850e-001 -#define s4 2.467759345363656348e-002 -/* - * Constants for approximation in [1.71, 2.5] -*/ -#define a1_hi 4.227843350984671344505727574870e-01 -#define a1_lo 4.670126436531227189e-18 -#define p0 3.224670334241133695662995251041e-01 -#define p1 3.569659696950364669021382724168e-01 -#define p2 1.342918716072560025853732668111e-01 -#define p3 1.950702176409779831089963408886e-02 -#define p4 8.546740251667538090796227834289e-04 -#define q0 1.000000000000000444089209850062e+00 -#define q1 1.315850076960161985084596381057e+00 -#define q2 6.274644311862156431658377186977e-01 -#define q3 1.304706631926259297049597307705e-01 -#define q4 1.102815279606722369265536798366e-02 -#define q5 2.512690594856678929537585620579e-04 -#define q6 -1.003597548112371003358107325598e-06 -/* - * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf]. -*/ -#define lns2pi .418938533204672741780329736405 -#define pb0 8.33333333333333148296162562474e-02 -#define pb1 -2.77777777774548123579378966497e-03 -#define pb2 7.93650778754435631476282786423e-04 -#define pb3 -5.95235082566672847950717262222e-04 -#define pb4 8.41428560346653702135821806252e-04 -#define pb5 -1.89773526463879200348872089421e-03 -#define pb6 5.69394463439411649408050664078e-03 -#define pb7 -1.44705562421428915453880392761e-02 - -__pure double -lgamma(double x) -{ - double r; - - int signgam = 1; -#if _IEEE - endian = ((*(int *) &one)) ? 1 : 0; -#endif - - if (!isfinite(x)) - if (_IEEE) - return (x+x); - else return (infnan(EDOM)); - - if (x > 6 + RIGHT) { - r = large_lgam(x); - return (r); - } else if (x > 1e-16) - return (small_lgam(x)); - else if (x > -1e-16) { - if (x < 0) { - signgam = -1; - x = -x; - } - return (-log(x)); - } else - return (neg_lgam(x)); -} - -float -lgammaf(float x) -{ - return lgamma(x); -} - -/* - * The gamma() function performs identically to lgamma(), including - * the use of signgam. - */ - -double -gamma(double x) -{ - return lgamma(x); -} - -float -gammaf(float x) -{ - return lgammaf(x); -} - -static double -large_lgam(double x) -{ - double z, p, x1; - struct Double t, u, v; - u = __log__D(x); - u.a -= 1.0; - if (x > 1e15) { - v.a = x - 0.5; - TRUNC(v.a); - v.b = (x - v.a) - 0.5; - t.a = u.a*v.a; - t.b = x*u.b + v.b*u.a; - if (_IEEE == 0 && !isfinite(t.a)) - return(infnan(ERANGE)); - return(t.a + t.b); - } - x1 = 1./x; - z = x1*x1; - p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7)))))); - /* error in approximation = 2.8e-19 */ - - p = p*x1; /* error < 2.3e-18 absolute */ - /* 0 < p < 1/64 (at x = 5.5) */ - v.a = x = x - 0.5; - TRUNC(v.a); /* truncate v.a to 26 bits. */ - v.b = x - v.a; - t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ - t.b = v.b*u.a + x*u.b; - t.b += p; t.b += lns2pi; /* return t + lns2pi + p */ - return (t.a + t.b); -} - -static double -small_lgam(double x) -{ - int x_int; - double y, z, t, r = 0, p, q, hi, lo; - struct Double rr; - x_int = (x + .5); - y = x - x_int; - if (x_int <= 2 && y > RIGHT) { - t = y - x0; - y--; x_int++; - goto CONTINUE; - } else if (y < -LEFT) { - t = y +(1.0-x0); -CONTINUE: - z = t - x0_lo; - p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); - q = s0+z*(s1+z*(s2+z*(s3+z*s4))); - r = t*(z*(p/q) - x0_lo); - t = .5*t*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); /* FALLTHROUGH */ - case 5: z *= (y + 4); /* FALLTHROUGH */ - case 4: z *= (y + 3); /* FALLTHROUGH */ - case 3: z *= (y + 2); - rr = __log__D(z); - rr.b += a0_lo; rr.a += a0_hi; - return(((r+rr.b)+t+rr.a)); - case 2: return(((r+a0_lo)+t)+a0_hi); - case 0: r -= log1p(x); /* FALLTHROUGH */ - default: rr = __log__D(x); - rr.a -= a0_hi; rr.b -= a0_lo; - return(((r - rr.b) + t) - rr.a); - } - } else { - p = p0+y*(p1+y*(p2+y*(p3+y*p4))); - q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6))))); - p = p*(y/q); - t = (double)(float) y; - z = y-t; - hi = (double)(float) (p+a1_hi); - lo = a1_hi - hi; lo += p; lo += a1_lo; - r = lo*y + z*hi; /* q + r = y*(a0+p/q) */ - q = hi*t; - z = 1.0; - switch (x_int) { - case 6: z = (y + 5); /* FALLTHROUGH */ - case 5: z *= (y + 4); /* FALLTHROUGH */ - case 4: z *= (y + 3); /* FALLTHROUGH */ - case 3: z *= (y + 2); /* FALLTHROUGH */ - rr = __log__D(z); - r += rr.b; r += q; - return(rr.a + r); - case 2: return (q+ r); - case 0: rr = __log__D(x); - r -= rr.b; r -= log1p(x); - r += q; r-= rr.a; - return(r); - default: rr = __log__D(x); - r -= rr.b; - q -= rr.a; - return (r+q); - } - } -} - -static double -neg_lgam(double x) -{ - int xi; - double y, z, zero = 0.0; - - /* avoid destructive cancellation as much as possible */ - if (x > -170) { - xi = x; - if (xi == x) - if (_IEEE) - return(one/zero); - else - return(infnan(ERANGE)); - y = tgamma(x); - if (y < 0) { - y = -y; - signgam = -1; - } - return (log(y)); - } - z = floor(x + .5); - if (z == x) { /* convention: G(-(integer)) -> +Inf */ - if (_IEEE) - return (one/zero); - else - return (infnan(ERANGE)); - } - y = .5*ceil(x); - if (y == ceil(y)) - signgam = -1; - x = -x; - z = fabs(x + z); /* 0 < z <= .5 */ - if (z < .25) - z = sin(M_PI*z); - else - z = cos(M_PI*(0.5-z)); - z = log(M_PI/(z*x)); - y = large_lgam(x); - return (z - y); -} diff --git a/lib/libm/noieee_src/n_log.c b/lib/libm/noieee_src/n_log.c deleted file mode 100644 index ce5a9f40801..00000000000 --- a/lib/libm/noieee_src/n_log.c +++ /dev/null @@ -1,482 +0,0 @@ -/* $OpenBSD: n_log.c,v 1.10 2016/09/12 19:47:02 guenther 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. 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. - */ - -#include -#include - -#include "mathimpl.h" - -/* Table-driven natural logarithm. - * - * This code was derived, with minor modifications, from: - * Peter Tang, "Table-Driven Implementation of the - * Logarithm in IEEE Floating-Point arithmetic." ACM Trans. - * Math Software, vol 16. no 4, pp 378-400, Dec 1990). - * - * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256, - * where F = j/128 for j an integer in [0, 128]. - * - * log(2^m) = log2_hi*m + log2_tail*m - * since m is an integer, the dominant term is exact. - * m has at most 10 digits (for subnormal numbers), - * and log2_hi has 11 trailing zero bits. - * - * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h - * logF_hi[] + 512 is exact. - * - * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ... - * the leading term is calculated to extra precision in two - * parts, the larger of which adds exactly to the dominant - * m and F terms. - * There are two cases: - * 1. when m, j are non-zero (m | j), use absolute - * precision for the leading term. - * 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1). - * In this case, use a relative precision of 24 bits. - * (This is done differently in the original paper) - * - * Special cases: - * 0 return signalling -Inf - * neg return signalling NaN - * +Inf return +Inf -*/ - -#if defined(__vax__) -#define _IEEE 0 -#define TRUNC(x) x = (double) (float) (x) -#else -#define _IEEE 1 -#define endian (((*(int *) &one)) ? 1 : 0) -#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -#define N 128 - -/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128. - * Used for generation of extend precision logarithms. - * The constant 35184372088832 is 2^45, so the divide is exact. - * It ensures correct reading of logF_head, even for inaccurate - * decimal-to-binary conversion routines. (Everybody gets the - * right answer for integers less than 2^53.) - * Values for log(F) were generated using error < 10^-57 absolute - * with the bc -l package. -*/ -static const double A1 = .08333333333333178827; -static const double A2 = .01250000000377174923; -static const double A3 = .002232139987919447809; -static const double A4 = .0004348877777076145742; - -static const double logF_head[N+1] = { - 0., - .007782140442060381246, - .015504186535963526694, - .023167059281547608406, - .030771658666765233647, - .038318864302141264488, - .045809536031242714670, - .053244514518837604555, - .060624621816486978786, - .067950661908525944454, - .075223421237524235039, - .082443669210988446138, - .089612158689760690322, - .096729626458454731618, - .103796793681567578460, - .110814366340264314203, - .117783035656430001836, - .124703478501032805070, - .131576357788617315236, - .138402322859292326029, - .145182009844575077295, - .151916042025732167530, - .158605030176659056451, - .165249572895390883786, - .171850256926518341060, - .178407657472689606947, - .184922338493834104156, - .191394852999565046047, - .197825743329758552135, - .204215541428766300668, - .210564769107350002741, - .216873938300523150246, - .223143551314024080056, - .229374101064877322642, - .235566071312860003672, - .241719936886966024758, - .247836163904594286577, - .253915209980732470285, - .259957524436686071567, - .265963548496984003577, - .271933715484010463114, - .277868451003087102435, - .283768173130738432519, - .289633292582948342896, - .295464212893421063199, - .301261330578199704177, - .307025035294827830512, - .312755710004239517729, - .318453731118097493890, - .324119468654316733591, - .329753286372579168528, - .335355541920762334484, - .340926586970454081892, - .346466767346100823488, - .351976423156884266063, - .357455888922231679316, - .362905493689140712376, - .368325561158599157352, - .373716409793814818840, - .379078352934811846353, - .384411698910298582632, - .389716751140440464951, - .394993808240542421117, - .400243164127459749579, - .405465108107819105498, - .410659924985338875558, - .415827895143593195825, - .420969294644237379543, - .426084395310681429691, - .431173464818130014464, - .436236766774527495726, - .441274560805140936281, - .446287102628048160113, - .451274644139630254358, - .456237433481874177232, - .461175715122408291790, - .466089729924533457960, - .470979715219073113985, - .475845904869856894947, - .480688529345570714212, - .485507815781602403149, - .490303988045525329653, - .495077266798034543171, - .499827869556611403822, - .504556010751912253908, - .509261901790523552335, - .513945751101346104405, - .518607764208354637958, - .523248143765158602036, - .527867089620485785417, - .532464798869114019908, - .537041465897345915436, - .541597282432121573947, - .546132437597407260909, - .550647117952394182793, - .555141507540611200965, - .559615787935399566777, - .564070138285387656651, - .568504735352689749561, - .572919753562018740922, - .577315365035246941260, - .581691739635061821900, - .586049045003164792433, - .590387446602107957005, - .594707107746216934174, - .599008189645246602594, - .603290851438941899687, - .607555250224322662688, - .611801541106615331955, - .616029877215623855590, - .620240409751204424537, - .624433288012369303032, - .628608659422752680256, - .632766669570628437213, - .636907462236194987781, - .641031179420679109171, - .645137961373620782978, - .649227946625615004450, - .653301272011958644725, - .657358072709030238911, - .661398482245203922502, - .665422632544505177065, - .669430653942981734871, - .673422675212350441142, - .677398823590920073911, - .681359224807238206267, - .685304003098281100392, - .689233281238557538017, - .693147180560117703862 -}; - -static const double logF_tail[N+1] = { - 0., - -.00000000000000543229938420049, - .00000000000000172745674997061, - -.00000000000001323017818229233, - -.00000000000001154527628289872, - -.00000000000000466529469958300, - .00000000000005148849572685810, - -.00000000000002532168943117445, - -.00000000000005213620639136504, - -.00000000000001819506003016881, - .00000000000006329065958724544, - .00000000000008614512936087814, - -.00000000000007355770219435028, - .00000000000009638067658552277, - .00000000000007598636597194141, - .00000000000002579999128306990, - -.00000000000004654729747598444, - -.00000000000007556920687451336, - .00000000000010195735223708472, - -.00000000000017319034406422306, - -.00000000000007718001336828098, - .00000000000010980754099855238, - -.00000000000002047235780046195, - -.00000000000008372091099235912, - .00000000000014088127937111135, - .00000000000012869017157588257, - .00000000000017788850778198106, - .00000000000006440856150696891, - .00000000000016132822667240822, - -.00000000000007540916511956188, - -.00000000000000036507188831790, - .00000000000009120937249914984, - .00000000000018567570959796010, - -.00000000000003149265065191483, - -.00000000000009309459495196889, - .00000000000017914338601329117, - -.00000000000001302979717330866, - .00000000000023097385217586939, - .00000000000023999540484211737, - .00000000000015393776174455408, - -.00000000000036870428315837678, - .00000000000036920375082080089, - -.00000000000009383417223663699, - .00000000000009433398189512690, - .00000000000041481318704258568, - -.00000000000003792316480209314, - .00000000000008403156304792424, - -.00000000000034262934348285429, - .00000000000043712191957429145, - -.00000000000010475750058776541, - -.00000000000011118671389559323, - .00000000000037549577257259853, - .00000000000013912841212197565, - .00000000000010775743037572640, - .00000000000029391859187648000, - -.00000000000042790509060060774, - .00000000000022774076114039555, - .00000000000010849569622967912, - -.00000000000023073801945705758, - .00000000000015761203773969435, - .00000000000003345710269544082, - -.00000000000041525158063436123, - .00000000000032655698896907146, - -.00000000000044704265010452446, - .00000000000034527647952039772, - -.00000000000007048962392109746, - .00000000000011776978751369214, - -.00000000000010774341461609578, - .00000000000021863343293215910, - .00000000000024132639491333131, - .00000000000039057462209830700, - -.00000000000026570679203560751, - .00000000000037135141919592021, - -.00000000000017166921336082431, - -.00000000000028658285157914353, - -.00000000000023812542263446809, - .00000000000006576659768580062, - -.00000000000028210143846181267, - .00000000000010701931762114254, - .00000000000018119346366441110, - .00000000000009840465278232627, - -.00000000000033149150282752542, - -.00000000000018302857356041668, - -.00000000000016207400156744949, - .00000000000048303314949553201, - -.00000000000071560553172382115, - .00000000000088821239518571855, - -.00000000000030900580513238244, - -.00000000000061076551972851496, - .00000000000035659969663347830, - .00000000000035782396591276383, - -.00000000000046226087001544578, - .00000000000062279762917225156, - .00000000000072838947272065741, - .00000000000026809646615211673, - -.00000000000010960825046059278, - .00000000000002311949383800537, - -.00000000000058469058005299247, - -.00000000000002103748251144494, - -.00000000000023323182945587408, - -.00000000000042333694288141916, - -.00000000000043933937969737844, - .00000000000041341647073835565, - .00000000000006841763641591466, - .00000000000047585534004430641, - .00000000000083679678674757695, - -.00000000000085763734646658640, - .00000000000021913281229340092, - -.00000000000062242842536431148, - -.00000000000010983594325438430, - .00000000000065310431377633651, - -.00000000000047580199021710769, - -.00000000000037854251265457040, - .00000000000040939233218678664, - .00000000000087424383914858291, - .00000000000025218188456842882, - -.00000000000003608131360422557, - -.00000000000050518555924280902, - .00000000000078699403323355317, - -.00000000000067020876961949060, - .00000000000016108575753932458, - .00000000000058527188436251509, - -.00000000000035246757297904791, - -.00000000000018372084495629058, - .00000000000088606689813494916, - .00000000000066486268071468700, - .00000000000063831615170646519, - .00000000000025144230728376072, - -.00000000000017239444525614834 -}; - -double -log(double x) -{ - int m, j; - double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0; - volatile double u1; - - /* Catch special cases */ - if (x <= 0) - if (_IEEE && x == zero) /* log(0) = -Inf */ - return (-one/zero); - else if (_IEEE) /* log(neg) = NaN */ - return (zero/zero); - else if (x == zero) /* NOT REACHED IF _IEEE */ - return (infnan(-ERANGE)); - else - return (infnan(EDOM)); - else if (!isfinite(x)) - if (_IEEE) /* x = NaN, Inf */ - return (x+x); - else - return (infnan(ERANGE)); - - /* Argument reduction: 1 <= g < 2; x/2^m = g; */ - /* y = F*(1 + f/F) for |f| <= 2^-8 */ - - m = logb(x); - g = ldexp(x, -m); - if (_IEEE && m == -1022) { - j = logb(g); - m += j; - g = ldexp(g, -j); - } - j = N*(g-1) + .5; - F = (1.0/N) * j + 1; /* F*128 is an integer in [128, 512] */ - f = g - F; - - /* Approximate expansion for log(1+f/F) ~= u + q */ - g = 1/(2*F+f); - u = 2*f*g; - v = u*u; - q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); - - /* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8, - * u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits. - * It also adds exactly to |m*log2_hi + log_F_head[j] | < 750 - */ - if (m | j) { - u1 = u + 513; - u1 -= 513; - } - - /* case 2: |1-x| < 1/256. The m- and j- dependent terms are zero; - * u1 = u to 24 bits. - */ - else { - u1 = u; - TRUNC(u1); - } - u2 = (2.0*(f - F*u1) - u1*f) * g; - /* u1 + u2 = 2f/(2F+f) to extra precision. */ - - /* log(x) = log(2^m*F*(1+f/F)) = */ - /* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q); */ - /* (exact) + (tiny) */ - - u1 += m*logF_head[N] + logF_head[j]; /* exact */ - u2 = (u2 + logF_tail[j]) + q; /* tiny */ - u2 += logF_tail[N]*m; - return (u1 + u2); -} -DEF_STD(log); - -/* - * Extra precision variant, returning struct {double a, b;}; - * log(x) = a+b to 63 bits, with a rounded to 26 bits. - */ -struct Double -__log__D(double x) -{ - int m, j; - double F, f, g, q, u, v, u2, one = 1.0; - volatile double u1; - struct Double r; - - /* Argument reduction: 1 <= g < 2; x/2^m = g; */ - /* y = F*(1 + f/F) for |f| <= 2^-8 */ - - m = logb(x); - g = ldexp(x, -m); - if (_IEEE && m == -1022) { - j = logb(g); - m += j; - g = ldexp(g, -j); - } - j = N*(g-1) + .5; - F = (1.0/N) * j + 1; - f = g - F; - - g = 1/(2*F+f); - u = 2*f*g; - v = u*u; - q = u*v*(A1 + v*(A2 + v*(A3 + v*A4))); - if (m | j) { - u1 = u + 513; - u1 -= 513; - } - else { - u1 = u; - TRUNC(u1); - } - u2 = (2.0*(f - F*u1) - u1*f) * g; - - u1 += m*logF_head[N] + logF_head[j]; - - u2 += logF_tail[j]; u2 += q; - u2 += logF_tail[N]*m; - r.a = u1 + u2; /* Only difference is here */ - TRUNC(r.a); - r.b = (u1 - r.a) + u2; - return (r); -} diff --git a/lib/libm/noieee_src/n_log10.c b/lib/libm/noieee_src/n_log10.c deleted file mode 100644 index 5c0126a0691..00000000000 --- a/lib/libm/noieee_src/n_log10.c +++ /dev/null @@ -1,83 +0,0 @@ -/* $OpenBSD: n_log10.c,v 1.8 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_log10.c,v 1.1 1995/10/10 23:36:58 ragge 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. 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. - */ - -/* LOG10(X) - * RETURN THE BASE 10 LOGARITHM OF x - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/20/85; - * REVISED BY K.C. NG on 1/23/85, 3/7/85, 4/16/85. - * - * Required kernel function: - * log(x) - * - * Method : - * log(x) - * log10(x) = --------- or [1/log(10)]*log(x) - * log(10) - * - * Note: - * [log(10)] rounded to 56 bits has error .0895 ulps, - * [1/log(10)] rounded to 53 bits has error .198 ulps; - * therefore, for better accuracy, in VAX D format, we divide - * log(x) by log(10), but in IEEE Double format, we multiply - * log(x) by [1/log(10)]. - * - * Special cases: - * log10(x) is NaN with signal if x < 0; - * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; - * log10(NaN) is that NaN with no signal. - * - * Accuracy: - * log10(X) returns the exact log10(x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was 1.74 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double ln10hi = 2.3025850929940456790E0; - -double -log10(double x) -{ -#if defined(__vax__) - return(log(x)/ln10hi); -#else /* defined(__vax__) */ - return(ivln10*log(x)); -#endif /* defined(__vax__) */ -} diff --git a/lib/libm/noieee_src/n_log1p.c b/lib/libm/noieee_src/n_log1p.c deleted file mode 100644 index 154cb1c4901..00000000000 --- a/lib/libm/noieee_src/n_log1p.c +++ /dev/null @@ -1,154 +0,0 @@ -/* $OpenBSD: n_log1p.c,v 1.12 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_log1p.c,v 1.1 1995/10/10 23:37:00 ragge 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. 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. - */ - -/* LOG1P(x) - * RETURN THE LOGARITHM OF 1+x - * DOUBLE PRECISION (VAX D FORMAT 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. NG on 2/6/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions: - * scalbn(x,n) - * copysign(x,y) - * logb(x) - * isfinite(x) - * - * Required kernel function: - * log__L(z) - * - * Method : - * 1. Argument Reduction: find k and f such that - * 1+x = 2^k * (1+f), - * where sqrt(2)/2 < 1+f < sqrt(2) . - * - * 2. Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * log(1+f) is computed by - * - * log(1+f) = 2s + s*log__L(s*s) - * where - * log__L(z) = z*(L1 + z*(L2 + z*(... (L6 + z*L7)...))) - * - * See log__L() for the values of the coefficients. - * - * 3. Finally, log(1+x) = k*ln2 + log(1+f). - * - * Remarks 1. In step 3 n*ln2 will be stored in two floating point numbers - * n*ln2hi + n*ln2lo, where ln2hi is chosen such that the last - * 20 bits (for VAX D format), or the last 21 bits ( for IEEE - * double) is 0. This ensures n*ln2hi is exactly representable. - * 2. In step 1, f may not be representable. A correction term c - * for f is computed. It follows that the correction term for - * f - t (the leading term of log(1+f) in step 2) is c-c*x. We - * add this correction term to n*ln2lo to attenuate the error. - * - * - * Special cases: - * log1p(x) is NaN with signal if x < -1; log1p(NaN) is NaN with no signal; - * log1p(INF) is +INF; log1p(-1) is -INF with signal; - * only log1p(0)=0 is exact for finite argument. - * - * Accuracy: - * log1p(x) returns the exact log(1+x) nearly rounded. In a test run - * with 1,536,000 random arguments on a VAX, the maximum observed - * error was .846 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include -#include "math.h" -#include "mathimpl.h" - -static const double ln2hi = 6.9314718055829871446E-1; -static const double ln2lo = 1.6465949582897081279E-12; -static const double sqrt2 = 1.4142135623730950622E0; - -double -log1p(double x) -{ - static const double zero=0.0, negone= -1.0, one=1.0, - half=1.0/2.0, small=1.0E-20; /* 1+small == 1 */ - double z,s,t,c; - int k; - - if (isnan(x)) - return (x); - - if(isfinite(x)) { - if( x > negone ) { - - /* argument reduction */ - if(copysign(x,one)= sqrt2 ) - { k += 1 ; z *= half; t *= half; } - t += negone; x = z + t; - c = (t-x)+z ; /* correction term for x */ - - /* compute log(1+x) */ - s = x/(2+x); t = x*x*half; - c += (k*ln2lo-c*x); - z = c+s*(t+__log__L(s*s)); - x += (z - t) ; - - return(k*ln2hi+x); - } - /* end of if (x > negone) */ - - else { -#if defined(__vax__) - if ( x == negone ) - return (infnan(-ERANGE)); /* -INF */ - else - return (infnan(EDOM)); /* NaN */ -#else /* defined(__vax__) */ - /* x = -1, return -INF with signal */ - if ( x == negone ) return( negone/zero ); - - /* negative argument for log, return NaN with signal */ - else return ( zero / zero ); -#endif /* defined(__vax__) */ - } - } - /* end of if (isfinite(x)) */ - - /* log(-INF) is NaN */ - else if(x<0) - return(zero/zero); - - /* log(+INF) is INF */ - else return(x); -} diff --git a/lib/libm/noieee_src/n_log__L.c b/lib/libm/noieee_src/n_log__L.c deleted file mode 100644 index 84fa2be7c8c..00000000000 --- a/lib/libm/noieee_src/n_log__L.c +++ /dev/null @@ -1,86 +0,0 @@ -/* $OpenBSD: n_log__L.c,v 1.9 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_log__L.c,v 1.1 1995/10/10 23:37:01 ragge 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. 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__L(Z) - * LOG(1+X) - 2S X - * RETURN --------------- WHERE Z = S*S, S = ------- , 0 <= Z <= .0294... - * S 2 + X - * - * DOUBLE PRECISION (VAX D FORMAT 56 bits or IEEE DOUBLE 53 BITS) - * KERNEL FUNCTION FOR LOG; TO BE USED IN LOG1P, LOG, AND POW FUNCTIONS - * CODED IN C BY K.C. NG, 1/19/85; - * REVISED BY K.C. Ng, 2/3/85, 4/16/85. - * - * Method : - * 1. Polynomial approximation: let s = x/(2+x). - * Based on log(1+x) = log(1+s) - log(1-s) - * = 2s + 2/3 s**3 + 2/5 s**5 + ....., - * - * (log(1+x) - 2s)/s is computed by - * - * z*(L1 + z*(L2 + z*(... (L7 + z*L8)...))) - * - * where z=s*s. (See the listing below for Lk's values.) The - * coefficients are obtained by a special Remes algorithm. - * - * Accuracy: - * Assuming no rounding error, the maximum magnitude of the approximation - * error (absolute) is 2**(-58.49) for IEEE double, and 2**(-63.63) - * for VAX D format. - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double L1 = 6.6666666666666703212E-1; -static const double L2 = 3.9999999999970461961E-1; -static const double L3 = 2.8571428579395698188E-1; -static const double L4 = 2.2222221233634724402E-1; -static const double L5 = 1.8181879517064680057E-1; -static const double L6 = 1.5382888777946145467E-1; -static const double L7 = 1.3338356561139403517E-1; -static const double L8 = 1.2500000000000000000E-1; - -double -__log__L(double z) -{ -#if defined(__vax__) - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*(L7+z*L8)))))))); -#else /* defined(__vax__) */ - return(z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7))))))); -#endif /* defined(__vax__) */ -} diff --git a/lib/libm/noieee_src/n_lrint.c b/lib/libm/noieee_src/n_lrint.c deleted file mode 100644 index a0feab5313f..00000000000 --- a/lib/libm/noieee_src/n_lrint.c +++ /dev/null @@ -1,57 +0,0 @@ -/* $OpenBSD: n_lrint.c,v 1.2 2016/09/12 19:47:02 guenther Exp $ */ -/* - * Copyright (c) 2013 Marc Espie. - * - * 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. - * - * THIS SOFTWARE IS PROVIDED BY THE OPENBSD PROJECT 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 OPENBSD - * PROJECT 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. - */ - -#include - -long -lrint(double x) -{ - return (long)rint(x); -} -DEF_STD(lrint); -LDBL_UNUSED_CLONE(lrint); - -long -lrintf(float x) -{ - return (long)rintf(x); -} -DEF_STD(lrintf); - -long long -llrint(double x) -{ - return (long long)rint(x); -} -DEF_STD(llrint); -LDBL_UNUSED_CLONE(llrint); - -long long -llrintf(float x) -{ - return (long long)rintf(x); -} -DEF_STD(llrintf); diff --git a/lib/libm/noieee_src/n_nan.c b/lib/libm/noieee_src/n_nan.c deleted file mode 100644 index 3b0cb57a4ba..00000000000 --- a/lib/libm/noieee_src/n_nan.c +++ /dev/null @@ -1,29 +0,0 @@ -/* $OpenBSD: n_nan.c,v 1.1 2009/03/28 16:16:30 martynas Exp $ */ - -/* - * Written by Martynas Venckus. Public domain - */ - -#include - -/* ARGSUSED */ -double -nan(const char *tagp) -{ - return (0); -} - -/* ARGSUSED */ -float -nanf(const char *tagp) -{ - return (0); -} - -/* ARGSUSED */ -long double -nanl(const char *tagp) -{ - return (0); -} - diff --git a/lib/libm/noieee_src/n_pow.c b/lib/libm/noieee_src/n_pow.c deleted file mode 100644 index eeb06e6c845..00000000000 --- a/lib/libm/noieee_src/n_pow.c +++ /dev/null @@ -1,205 +0,0 @@ -/* $OpenBSD: n_pow.c,v 1.13 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_pow.c,v 1.1 1995/10/10 23:37:02 ragge 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. 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. - */ - -/* POW(X,Y) - * RETURN X**Y - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 7/10/85. - * KERNEL pow_P() REPLACED BY P. McILROY 7/22/92. - * Required system supported functions: - * scalbn(x,n) - * logb(x) - * copysign(x,y) - * isfinite(x) - * remainder(x,y) - * - * Required kernel functions: - * exp__D(a,c) exp(a + c) for |a| << |c| - * struct d_double dlog(x) r.a + r.b, |r.b| < |r.a| - * - * Method - * 1. Compute and return log(x) in three pieces: - * log(x) = n*ln2 + hi + lo, - * where n is an integer. - * 2. Perform y*log(x) by simulating muti-precision arithmetic and - * return the answer in three pieces: - * y*log(x) = m*ln2 + hi + lo, - * where m is an integer. - * 3. Return x**y = exp(y*log(x)) - * = 2^m * ( exp(hi+lo) ). - * - * Special cases: - * (anything) ** 0 is 1 ; - * (anything) ** 1 is itself; - * (anything) ** NaN is NaN; - * NaN ** (anything except 0) is NaN; - * +(anything > 1) ** +INF is +INF; - * -(anything > 1) ** +INF is NaN; - * +-(anything > 1) ** -INF is +0; - * +-(anything < 1) ** +INF is +0; - * +(anything < 1) ** -INF is +INF; - * -(anything < 1) ** -INF is NaN; - * +-1 ** +-INF is NaN and signal INVALID; - * +0 ** +(anything except 0, NaN) is +0; - * -0 ** +(anything except 0, NaN, odd integer) is +0; - * +0 ** -(anything except 0, NaN) is +INF and signal DIV-BY-ZERO; - * -0 ** -(anything except 0, NaN, odd integer) is +INF with signal; - * -0 ** (odd integer) = -( +0 ** (odd integer) ); - * +INF ** +(anything except 0,NaN) is +INF; - * +INF ** -(anything except 0,NaN) is +0; - * -INF ** (odd integer) = -( +INF ** (odd integer) ); - * -INF ** (even integer) = ( +INF ** (even integer) ); - * -INF ** -(anything except integer,NaN) is NaN with signal; - * -(x=anything) ** (k=integer) is (-1)**k * (x ** k); - * -(anything except 0) ** (non-integer) is NaN with signal; - * - * Accuracy: - * pow(x,y) returns x**y nearly rounded. In particular, on a SUN, a VAX, - * and a Zilog Z8000, - * pow(integer,integer) - * always returns the correct integer provided it is representable. - * In a test run with 100,000 random arguments with 0 < x, y < 20.0 - * on a VAX, the maximum observed error was 1.79 ulps (units in the - * last place). - * - * Constants : - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include -#include - -#include "mathimpl.h" - -#if defined(__vax__) -#define TRUNC(x) x = (double) (float) x -#define _IEEE 0 -#else -#define _IEEE 1 -#define endian (((*(int *) &one)) ? 1 : 0) -#define TRUNC(x) *(((int *) &x)+endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif /* defined(__vax__) */ - -static const double zero=0.0, one=1.0, two=2.0, negone= -1.0; - -static double pow_P(double, double); - -double -pow(double x, double y) -{ - double t; - if (y==zero) - return (one); - else if (y==one || isnan(x)) - return (x); /* if x is NaN or y=1 */ - else if (isnan(y)) /* if y is NaN */ - return (y); - else if (!isfinite(y)) /* if y is INF */ - if ((t=fabs(x))==one) /* +-1 ** +-INF is NaN */ - return (y - y); - else if (t>one) - return ((y<0)? zero : ((x0)? zero : ((x<0)? y-y : -y)); - else if (y==two) - return (x*x); - else if (y==negone) - return (one/x); - /* x > 0, x == +0 */ - else if (copysign(one, x) == one) - return (pow_P(x, y)); - - /* sign(x)= -1 */ - /* if y is an even integer */ - else if ( (t=remainder(y,two)) == zero) - return (pow_P(-x, y)); - - /* if y is an odd integer */ - else if (copysign(t,one) == one) - return (-pow_P(-x, y)); - - /* Henceforth y is not an integer */ - else if (x==zero) /* x is -0 */ - return ((y>zero)? -x : one/(-x)); - else if (_IEEE) - return (zero/zero); - else - return (infnan(EDOM)); -} -/* kernel function for x >= 0 */ -static double -pow_P(double x, double y) -{ - struct Double s, t; - double huge = 1e300, tiny = 1e-300; - - if (x == zero) - if (y > zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (x == one) - return (one); - if (!isfinite(x)) - if (y < zero) - return (zero); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - if (y >= 7e18) /* infinity */ - if (x < 1) - return(tiny*tiny); - else if (_IEEE) - return (huge*huge); - else - return (infnan(ERANGE)); - - /* Return exp(y*log(x)), using simulated extended */ - /* precision for the log and the multiply. */ - - s = __log__D(x); - t.a = y; - TRUNC(t.a); - t.b = y - t.a; - t.b = s.b*y + t.b*s.a; - t.a *= s.a; - s.a = t.a + t.b; - s.b = (t.a - s.a) + t.b; - return (__exp__D(s.a, s.b)); -} diff --git a/lib/libm/noieee_src/n_sincos.c b/lib/libm/noieee_src/n_sincos.c deleted file mode 100644 index a3da96c6b9e..00000000000 --- a/lib/libm/noieee_src/n_sincos.c +++ /dev/null @@ -1,110 +0,0 @@ -/* $OpenBSD: n_sincos.c,v 1.16 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_sincos.c,v 1.1 1995/10/10 23:37:04 ragge Exp $ */ -/* - * Copyright (c) 1987, 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. 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. - */ - -#include - -#include "mathimpl.h" - -float -sinf(float x) -{ - return (float)sin((double) x); -} -DEF_STD(sinf); - -double -sin(double x) -{ - double a,c,z; - - if(!isfinite(x)) /* sin(NaN) and sin(INF) must be NaN */ - return x-x; - x=remainder(x,PI2); /* reduce x into [-PI,PI] */ - a=copysign(x,one); - if (a >= PIo4) { - if(a >= PI3o4) /* ... in [3PI/4,PI] */ - x = copysign((a = PI-a),x); - else { /* ... in [PI/4,3PI/4] */ - a = PIo2-a; /* rtn. sign(x)*C(PI/2-|x|) */ - z = a*a; - c = cos__C(z); - z *= half; - a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); - return copysign(a,x); - } - } - - if (a < small) { /* rtn. S(x) */ - big+a; - return x; - } - return x+x*sin__S(x*x); -} -DEF_STD(sin); -LDBL_CLONE(sin); - -float -cosf(float x) -{ - return (float)cos((double) x); -} - -double -cos(double x) -{ - double a,c,z,s = 1.0; - - if(!isfinite(x)) /* cos(NaN) and cos(INF) must be NaN */ - return x-x; - x=remainder(x,PI2); /* reduce x into [-PI,PI] */ - a=copysign(x,one); - if (a >= PIo4) { - if (a >= PI3o4) { /* ... in [3PI/4,PI] */ - a = PI-a; - s = negone; - } - else { /* ... in [PI/4,3PI/4] */ - a = PIo2-a; - return a+a*sin__S(a*a); /* rtn. S(PI/2-|x|) */ - } - } - if (a < small) { - big+a; - return s; /* rtn. s*C(a) */ - } - z = a*a; - c = cos__C(z); - z *= half; - a = (z >= thresh ? half-((z-half)-c) : one-(z-c)); - return copysign(a,s); -} -END_STD(cos); -LDBL_CLONE(cos); diff --git a/lib/libm/noieee_src/n_sinh.c b/lib/libm/noieee_src/n_sinh.c deleted file mode 100644 index 1fb55addabe..00000000000 --- a/lib/libm/noieee_src/n_sinh.c +++ /dev/null @@ -1,106 +0,0 @@ -/* $OpenBSD: n_sinh.c,v 1.10 2009/10/27 23:59:29 deraadt Exp $ */ -/* $NetBSD: n_sinh.c,v 1.1 1995/10/10 23:37:05 ragge 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. 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. - */ - -/* SINH(X) - * RETURN THE HYPERBOLIC SINE OF X - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 3/7/85, 3/24/85, 4/16/85. - * - * Required system supported functions : - * copysign(x,y) - * scalbn(x,N) - * - * Required kernel functions: - * expm1(x) ...return exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by sinh(-x) = - sinh(x). - * 2. - * - * expm1(x) + expm1(x)/(expm1(x)+1) - * 0 <= x <= lnovfl : sinh(x) := -------------------------------- - * 2 - * lnovfl <= x <= lnovfl+ln2 : sinh(x) := expm1(x)/2 (avoid overflow) - * lnovfl+ln2 < x < INF : overflow to INF - * - * - * Special cases: - * sinh(x) is x if x is +INF, -INF, or NaN. - * only sinh(0)=0 is exact for finite argument. - * - * Accuracy: - * sinh(x) returns the exact hyperbolic sine of x nearly rounded. In - * a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 1.93 ulps (units in the last place). - * - * Constants: - * The hexadecimal values are the intended ones for the following constants. - * The decimal values may be used, provided that the compiler will convert - * from decimal to binary accurately enough to produce the hexadecimal values - * shown. - */ - -#include "math.h" -#include "mathimpl.h" - -static const double mln2hi = 8.8029691931113054792E1; -static const double mln2lo = -4.9650192275318476525E-16; -static const double lnovfl = 8.8029691931113053016E1; - -#if defined(__vax__) -static max = 126 ; -#else /* defined(__vax__) */ -static max = 1023 ; -#endif /* defined(__vax__) */ - -double -sinh(double x) -{ - static const double one=1.0, half=1.0/2.0 ; - double t, sign; - - if (isnan(x)) - return (x); - - sign=copysign(one,x); - x=copysign(x,one); - if(x - -#include "mathimpl.h" - -#if defined(__vax__) /* VAX D format */ -#include - static const unsigned short msign=0x7fff, mexp =0x7f80 ; - static const short prep1=57, gap=7, bias=129 ; - static const double novf=1.7E38, nunf=3.0E-39; -#else /* defined(__vax__) */ - static const unsigned short msign=0x7fff, mexp =0x7ff0 ; - static const short prep1=54, gap=4, bias=1023 ; - static const double novf=1.7E308, nunf=3.0E-308; -#endif /* defined(__vax__) */ - -static const double zero = 0.0; - -double -scalbn(double x, int N) -{ - int k; - unsigned short *px=(unsigned short *) &x; - - if( x == zero ) return(x); - -#if defined(__vax__) - if( (k= *px & mexp ) != ~msign ) { - if (N < -260) - return(nunf*nunf); - else if (N > 260) { - return(copysign(infnan(ERANGE),x)); - } -#else /* defined(__vax__) */ - if( (k= *px & mexp ) != mexp ) { - if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); - if( k == 0 ) { - x *= scalbn(1.0,(int)prep1); N -= prep1; return(scalbn(x,N));} -#endif /* defined(__vax__) */ - - if((k = (k>>gap)+ N) > 0 ) - if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k< -prep1 ) - /* gradual underflow */ - {*px=(*px&~mexp)|(short)(1<>gap)-bias); -#else /* defined(__vax__) */ - if( (k= *px & mexp ) != mexp ) - if ( k != 0 ) - return ( (k>>gap) - bias ); - else if( x != zero) - return ( -1022.0 ); - else - return(-(1.0/zero)); - else if(isnan(x)) - return(x); - else - {*px &= msign; return(x);} -#endif /* defined(__vax__) */ -} -DEF_STD(logb); -LDBL_UNUSED_CLONE(logb); - -double -remainder(double x, double p) -{ - short sign; - double hp, dp, tmp; - unsigned short k; - unsigned short - *px=(unsigned short *) &x , - *pp=(unsigned short *) &p , - *pd=(unsigned short *) &dp , - *pt=(unsigned short *) &tmp; - - *pp &= msign ; - -#if defined(__vax__) - if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ -#else /* defined(__vax__) */ - if( ( *px & mexp ) == mexp ) -#endif /* defined(__vax__) */ - return (x-p)-(x-p); /* create nan if x is inf */ - if (p == zero) { -#if defined(__vax__) - return(infnan(EDOM)); -#else /* defined(__vax__) */ - return zero/zero; -#endif /* defined(__vax__) */ - } - -#if defined(__vax__) - if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ -#else /* defined(__vax__) */ - if( ( *pp & mexp ) == mexp ) -#endif /* defined(__vax__) */ - { if (p != p) return p; else return x;} - - else if ( ((*pp & mexp)>>gap) <= 1 ) - /* subnormal p, or almost subnormal p */ - { double b; b=scalbn(1.0,(int)prep1); - p *= b; x = remainder(x,p); x *= b; return(remainder(x,p)/b);} - else if ( p >= novf/2) - { p /= 2 ; x /= 2; return(remainder(x,p)*2);} - else - { - dp=p+p; hp=p/2; - sign= *px & ~msign ; - *px &= msign ; - while ( x > dp ) - { - k=(*px & mexp) - (*pd & mexp) ; - tmp = dp ; - *pt += k ; - -#if defined(__vax__) - if( x < tmp ) *pt -= 128 ; -#else /* defined(__vax__) */ - if( x < tmp ) *pt -= 16 ; -#endif /* defined(__vax__) */ - - x -= tmp ; - } - if ( x > hp ) - { x -= p ; if ( x >= hp ) x -= p ; } - -#if defined(__vax__) - if (x) -#endif /* defined(__vax__) */ - *px ^= sign; - return( x); - - } -} -DEF_STD(remainder); - -/* The drem() function is a deprecated alias for remainder(). */ - -double -drem(double x, double p) -{ - return remainder(x, p); -} - -float -sqrtf(float x) -{ - return (float)sqrt((double) x); -} -DEF_STD(sqrtf); - -double -sqrt(double x) -{ - double q, s, b, r; - double t; - int m, n, i; -#if defined(__vax__) - int k=54; -#else /* defined(__vax__) */ - int k=51; -#endif /* defined(__vax__) */ - - /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ - if(isnan(x) || x == zero) return(x); - - /* sqrt(negative) is invalid */ - if(x= 100) { /* trigger inexact flag */ - if(s1.0) t=1; /* b>1 : Round-to-(+INF) */ - if(t>=0) q+=r; } /* else: Round-to-nearest */ - else { - s *= 2; x *= 4; - t = (x-s)-1; - b=1.0+3*r/4; if(b==1.0) goto end; - b=1.0+r/4; if(b>1.0) t=1; - if(t>=0) q+=r; } - } - -end: return(scalbn(q,n)); -} -DEF_STD(sqrt); -LDBL_CLONE(sqrt); - -#if 0 -/* REMAINDER(X,Y) - * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) - * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) - * INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions (not on a VAX): - * swapINX(i): save inexact flag and reset it to "i" - * swapENI(e): save inexact enable and reset it to "e" - */ - -double -remainder(double x, double y) -{ - static const n0=0, n1=1, n2=2, n3=3; - static const unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; - double hy, y1, t, t1; - short k; - long n; - int i, e; - unsigned short xexp,yexp, *px =(unsigned short *) &x , - nx,nf, *py =(unsigned short *) &y , - sign, *pt =(unsigned short *) &t , - *pt1 =(unsigned short *) &t1 ; - - xexp = px[n0] & mexp ; /* exponent of x */ - yexp = py[n0] & mexp ; /* exponent of y */ - sign = px[n0] &0x8000; /* sign of x */ - -/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ - if(isnan(x)) return(x); if(isnan(y)) return(y); /* x or y is NaN */ - if( xexp == mexp ) return(zero/zero); /* x is INF */ - if(y==zero) return(y/y); - -/* save the inexact flag and inexact enable in i and e respectively - * and reset them to zero - */ - i=swapINX(0); e=swapENI(0); - -/* subnormal number */ - nx=0; - if (yexp == 0) { - t = 1.0; - pt[n0] += m57; - y *= t; - nx = m57; - } - -/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ - if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} - - nf=nx; - py[n0] &= 0x7fff; - px[n0] &= 0x7fff; - -/* mask off the least significant 27 bits of y */ - t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; - -/* LOOP: argument reduction on x whenever x > y */ -loop: - while ( x > y ) - { - t=y; - t1=y1; - xexp=px[n0]&mexp; /* exponent of x */ - k=xexp-yexp-m25; - if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ - {pt[n0]+=k;pt1[n0]+=k;} - n=x/t; x=(x-n*t1)-n*(t-t1); - } - /* end while (x > y) */ - - if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} - -/* final adjustment */ - - hy=y/2.0; - if(x>hy||((x==hy)&&n%2==1)) x-=y; - px[n0] ^= sign; - if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} - -/* restore inexact flag and inexact enable */ - swapINX(i); swapENI(e); - - return(x); -} -DEF_STD(remainder); -#endif - -#if 0 -/* SQRT - * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT - * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE - * CODED IN C BY K.C. NG, 3/22/85. - * - * Warning: this code should not get compiled in unless ALL of - * the following machine-dependent routines are supplied. - * - * Required machine dependent functions: - * swapINX(i) ...return the status of INEXACT flag and reset it to "i" - * swapRM(r) ...return the current Rounding Mode and reset it to "r" - * swapENI(e) ...return the status of inexact enable and reset it to "e" - * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer - * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer - */ - -static const unsigned long table[] = { -0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, -58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, -21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; - -double -newsqrt(double x) -{ - double y, z, t, addc(), subc(); - double const b54=134217728.*134217728.; /* b54=2**54 */ - long mx, scalx; - long const mexp=0x7ff00000; - int i, j, r, e, swapINX(), swapRM(), swapENI(); - unsigned long *py=(unsigned long *) &y , - *pt=(unsigned long *) &t , - *px=(unsigned long *) &x ; - const int n0=0, n1=1; -/* Rounding Mode: RN ...round-to-nearest - * RZ ...round-towards 0 - * RP ...round-towards +INF - * RM ...round-towards -INF - */ - const int RN=0, RZ=1, RP=2, RM=3; - /* machine dependent: work on a Zilog Z8070 - * and a National 32081 & 16081 - */ - -/* exceptions */ - if(isnan(x) || x == 0.0) return(x); /* sqrt(NaN) is NaN, - sqrt(+-0) = +-0 */ - if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ - if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ - -/* save, reset, initialize */ - e=swapENI(0); /* ...save and reset the inexact enable */ - i=swapINX(0); /* ...save INEXACT flag */ - r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ - scalx=0; - -/* subnormal number, scale up x to x*2**54 */ - if(mx==0) {x *= b54 ; scalx-=0x01b00000;} - -/* scale x to avoid intermediate over/underflow: - * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ - if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} - if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} - -/* magic initial approximation to almost 8 sig. bits */ - py[n0]=(px[n0]>>1)+0x1ff80000; - py[n0]=py[n0]-table[(py[n0]>>15)&31]; - -/* Heron's rule once with correction to improve y to almost 18 sig. bits */ - t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; - -/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ - t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; - t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; - -/* twiddle last bit to force y correctly rounded */ - swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ - swapINX(0); /* ...clear INEXACT flag */ - swapENI(e); /* ...restore inexact enable status */ - t=x/y; /* ...chopped quotient, possibly inexact */ - j=swapINX(i); /* ...read and restore inexact flag */ - if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ - x=b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ - if(r==RN) t=addc(t); /* ...t=t+ulp */ - else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ - y=y+t; /* ...chopped sum */ - py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ -end: py[n0]=py[n0]+scalx; /* ...scale back y */ - swapRM(r); /* ...restore Rounding Mode */ - return(y); -} -#endif diff --git a/lib/libm/noieee_src/n_tan.c b/lib/libm/noieee_src/n_tan.c deleted file mode 100644 index 971ebe7a89f..00000000000 --- a/lib/libm/noieee_src/n_tan.c +++ /dev/null @@ -1,74 +0,0 @@ -/* $OpenBSD: n_tan.c,v 1.16 2016/09/12 19:47:02 guenther Exp $ */ -/* $NetBSD: n_tan.c,v 1.1 1995/10/10 23:37:07 ragge Exp $ */ -/* - * Copyright (c) 1987, 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. 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. - */ - -#include - -#include "mathimpl.h" - -float -tanf(float x) -{ - return (float)tan((double) x); -} - -double -tan(double x) -{ - double a,z,ss,cc,c; - int k; - - if(!isfinite(x)) /* tan(NaN) and tan(INF) must be NaN */ - return x-x; - x = remainder(x,PI); /* reduce x into [-PI/2, PI/2] */ - a = copysign(x,one); /* ... = abs(x) */ - if (a >= PIo4) { - k = 1; - x = copysign(PIo2-a,x); - } - else { - k = 0; - if (a < small) { - big+a; - return x; - } - } - z = x*x; - cc = cos__C(z); - ss = sin__S(z); - z *= half; /* Next get c = cos(x) accurately */ - c = (z >= thresh ? half-((z-half)-cc) : one-(z-cc)); - if (k == 0) - return x+(x*(z-(cc-ss)))/c; /* ... sin/cos */ - else - return c/(x+x*ss); /* ... cos/sin */ -} -DEF_STD(tan); -LDBL_UNUSED_CLONE(tan); diff --git a/lib/libm/noieee_src/n_tanh.c b/lib/libm/noieee_src/n_tanh.c deleted file mode 100644 index 3db4a516e63..00000000000 --- a/lib/libm/noieee_src/n_tanh.c +++ /dev/null @@ -1,95 +0,0 @@ -/* $OpenBSD: n_tanh.c,v 1.10 2016/09/12 04:39:47 guenther Exp $ */ -/* $NetBSD: n_tanh.c,v 1.1 1995/10/10 23:37:08 ragge 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. 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. - */ - -#include "math.h" - -/* TANH(X) - * RETURN THE HYPERBOLIC TANGENT OF X - * DOUBLE PRECISION (VAX D FORMAT 56 BITS, IEEE DOUBLE 53 BITS) - * CODED IN C BY K.C. NG, 1/8/85; - * REVISED BY K.C. NG on 2/8/85, 2/11/85, 3/7/85, 3/24/85. - * - * Required system supported functions : - * copysign(x,y) - * isfinite(x) - * - * Required kernel function: - * expm1(x) ...exp(x)-1 - * - * Method : - * 1. reduce x to non-negative by tanh(-x) = - tanh(x). - * 2. - * 0 < x <= 1.e-10 : tanh(x) := x - * -expm1(-2x) - * 1.e-10 < x <= 1 : tanh(x) := -------------- - * expm1(-2x) + 2 - * 2 - * 1 <= x <= 22.0 : tanh(x) := 1 - --------------- - * expm1(2x) + 2 - * 22.0 < x <= INF : tanh(x) := 1. - * - * Note: 22 was chosen so that fl(1.0+2/(expm1(2*22)+2)) == 1. - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - * - * Accuracy: - * tanh(x) returns the exact hyperbolic tangent of x nealy rounded. - * In a test run with 1,024,000 random arguments on a VAX, the maximum - * observed error was 2.22 ulps (units in the last place). - */ - -double -tanh(double x) -{ - static double one=1.0, two=2.0, small = 1.0e-10, big = 1.0e10; - double t, sign; - - if (isnan(x)) - return (x); - - sign=copysign(one,x); - x=copysign(x,one); - if(x < 22.0) - if( x > one ) - return(copysign(one-two/(expm1(x+x)+two),sign)); - else if ( x > small ) - {t= -expm1(-(x+x)); return(copysign(t/(two-t),sign));} - else { /* raise the INEXACT flag for non-zero x */ - t = big + x; - return(copysign(x,sign) - (t-(big+x))); - } - else if(isfinite(x)) - return (sign+1.0E-37); /* raise the INEXACT flag */ - else - return(sign); /* x is +- INF */ -} diff --git a/lib/libm/noieee_src/n_tgamma.c b/lib/libm/noieee_src/n_tgamma.c deleted file mode 100644 index 1bd2f13bbee..00000000000 --- a/lib/libm/noieee_src/n_tgamma.c +++ /dev/null @@ -1,363 +0,0 @@ -/* $OpenBSD: n_tgamma.c,v 1.5 2016/09/12 04:39:47 guenther 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. 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. - */ - -/* - * This code by P. McIlroy, Oct 1992; - * - * The financial support of UUNET Communications Services is greatfully - * acknowledged. - */ - -#include -#include "mathimpl.h" -#include - -/* METHOD: - * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) - * At negative integers, return +Inf and raise overflow. - * - * x < 6.5: - * Use argument reduction G(x+1) = xG(x) to reach the - * range [1.066124,2.066124]. Use a rational - * approximation centered at the minimum (x0+1) to - * ensure monotonicity. - * - * x >= 6.5: Use the asymptotic approximation (Stirling's formula) - * adjusted for equal-ripples: - * - * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x)) - * - * Keep extra precision in multiplying (x-.5)(log(x)-1), to - * avoid premature round-off. - * - * Special values: - * -Inf: return 0 and raise invalid; - * negative integer: return +Inf and raise overflow; - * other x ~< -177.79: return 0 and raise underflow; - * +-0: return +Inf and raise overflow; - * finite x ~> 171.63: return +Inf and raise div-by-0; - * +Inf: return +Inf and raise div-by-0; - * NaN: return 0 and raise invalid. - * - * Accuracy: tgamma(x) is accurate to within - * x > 0: error provably < 0.9ulp. - * Maximum observed in 1,000,000 trials was .87ulp. - * x < 0: - * Maximum observed error < 4ulp in 1,000,000 trials. - */ - -static double neg_gam(double); -static double small_gam(double); -static double smaller_gam(double); -static struct Double large_gam(double); -static struct Double ratfun_gam(double, double); - -/* - * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval - * [1.066.., 2.066..] accurate to 4.25e-19. - */ -#define LEFT -.3955078125 /* left boundary for rat. approx */ -#define x0 .461632144968362356785 /* xmin - 1 */ - -#define a0_hi 0.88560319441088874992 -#define a0_lo -.00000000000000004996427036469019695 -#define P0 6.21389571821820863029017800727e-01 -#define P1 2.65757198651533466104979197553e-01 -#define P2 5.53859446429917461063308081748e-03 -#define P3 1.38456698304096573887145282811e-03 -#define P4 2.40659950032711365819348969808e-03 -#define Q0 1.45019531250000000000000000000e+00 -#define Q1 1.06258521948016171343454061571e+00 -#define Q2 -2.07474561943859936441469926649e-01 -#define Q3 -1.46734131782005422506287573015e-01 -#define Q4 3.07878176156175520361557573779e-02 -#define Q5 5.12449347980666221336054633184e-03 -#define Q6 -1.76012741431666995019222898833e-03 -#define Q7 9.35021023573788935372153030556e-05 -#define Q8 6.13275507472443958924745652239e-06 -/* - * Constants for large x approximation (x in [6, Inf]) - * (Accurate to 2.8*10^-19 absolute) - */ -#define lns2pi_hi 0.418945312500000 -#define lns2pi_lo -.000006779295327258219670263595 -#define Pa0 8.33333333333333148296162562474e-02 -#define Pa1 -2.77777777774548123579378966497e-03 -#define Pa2 7.93650778754435631476282786423e-04 -#define Pa3 -5.95235082566672847950717262222e-04 -#define Pa4 8.41428560346653702135821806252e-04 -#define Pa5 -1.89773526463879200348872089421e-03 -#define Pa6 5.69394463439411649408050664078e-03 -#define Pa7 -1.44705562421428915453880392761e-02 - -static const double zero = 0., one = 1.0, tiny = 1e-300; - -/* - * TRUNC sets trailing bits in a floating-point number to zero. - * is a temporary variable. - */ -#if defined(__vax__) -#define _IEEE 0 -#define TRUNC(x) x = (double) (float) (x) -#else -static int endian; -#define _IEEE 1 -#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000 -#define infnan(x) 0.0 -#endif - -double -tgamma(double x) -{ - struct Double u; -#if _IEEE - endian = (*(int *) &one) ? 1 : 0; -#endif - - if (x >= 6) { - if(x > 171.63) - if (_IEEE) - return (x/zero); - else - return (infnan(ERANGE)); - u = large_gam(x); - return(__exp__D(u.a, u.b)); - } else if (x >= 1.0 + LEFT + x0) - return (small_gam(x)); - else if (x > 1.e-17) - return (smaller_gam(x)); - else if (x > -1.e-17) { - if (x == 0.0) { - if (!_IEEE) - return (infnan(ERANGE)); - } else { - u.a = one - tiny; /* raise inexact */ - } - return (one/x); - } else if (!isfinite(x)) { - if (_IEEE) /* x = NaN, -Inf */ - return (x - x); - else - return (infnan(EDOM)); - } else - return (neg_gam(x)); -} - -/* - * We simply call tgamma() rather than bloating the math library - * with a float-optimized version of it. The reason is that tgammaf() - * is essentially useless, since the function is superexponential - * and floats have very limited range. -- das@freebsd.org - */ - -float -tgammaf(float x) -{ - return tgamma(x); -} - -/* - * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error. - */ - -static struct Double -large_gam(double x) -{ - double z, p; - struct Double t, u, v; - - z = one/(x*x); - p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7)))))); - p = p/x; - - u = __log__D(x); - u.a -= one; - v.a = (x -= .5); - TRUNC(v.a); - v.b = x - v.a; - t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ - t.b = v.b*u.a + x*u.b; - /* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */ - t.b += lns2pi_lo; t.b += p; - u.a = lns2pi_hi + t.b; u.a += t.a; - u.b = t.a - u.a; - u.b += lns2pi_hi; u.b += t.b; - return (u); -} - -/* - * Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.) - * It also has correct monotonicity. - */ - -static double -small_gam(double x) -{ - double y, ym1, t; - struct Double yy, r; - y = x - one; - ym1 = y - one; - if (y <= 1.0 + (LEFT + x0)) { - yy = ratfun_gam(y - x0, 0); - return (yy.a + yy.b); - } - r.a = y; - TRUNC(r.a); - yy.a = r.a - one; - y = ym1; - yy.b = r.b = y - yy.a; - /* Argument reduction: G(x+1) = x*G(x) */ - for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) { - t = r.a*yy.a; - r.b = r.a*yy.b + y*r.b; - r.a = t; - TRUNC(r.a); - r.b += (t - r.a); - } - /* Return r*tgamma(y). */ - yy = ratfun_gam(y - x0, 0); - y = r.b*(yy.a + yy.b) + r.a*yy.b; - y += yy.a*r.a; - return (y); -} - -/* - * Good on (0, 1+x0+LEFT]. Accurate to 1ulp. - */ - -static double -smaller_gam(double x) -{ - double t, d; - struct Double r, xx; - if (x < x0 + LEFT) { - t = x; - TRUNC(t); - d = (t+x)*(x-t); - t *= t; - xx.a = (t + x); - TRUNC(xx.a); - xx.b = x - xx.a; xx.b += t; xx.b += d; - t = (one-x0); t += x; - d = (one-x0); d -= t; d += x; - x = xx.a + xx.b; - } else { - xx.a = x; - TRUNC(xx.a); - xx.b = x - xx.a; - t = x - x0; - d = (-x0 -t); d += x; - } - r = ratfun_gam(t, d); - d = r.a/x; - TRUNC(d); - r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b; - return (d + r.a/x); -} - -/* - * returns (z+c)^2 * P(z)/Q(z) + a0 - */ - -static struct Double -ratfun_gam(double z, double c) -{ - double p, q; - struct Double r, t; - - q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8))))))); - p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4))); - - /* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */ - p = p/q; - t.a = z; - TRUNC(t.a); /* t ~= z + c */ - t.b = (z - t.a) + c; - t.b *= (t.a + z); - q = (t.a *= t.a); /* t = (z+c)^2 */ - TRUNC(t.a); - t.b += (q - t.a); - r.a = p; - TRUNC(r.a); /* r = P/Q */ - r.b = p - r.a; - t.b = t.b*p + t.a*r.b + a0_lo; - t.a *= r.a; /* t = (z+c)^2*(P/Q) */ - r.a = t.a + a0_hi; - TRUNC(r.a); - r.b = ((a0_hi-r.a) + t.a) + t.b; - return (r); /* r = a0 + t */ -} - -static double -neg_gam(double x) -{ - int sgn = 1; - struct Double lg, lsine; - double y, z; - - y = ceil(x); - if (y == x) /* Negative integer. */ - if (_IEEE) - return ((x - x) / zero); - else - return (infnan(ERANGE)); - z = y - x; - if (z > 0.5) - z = one - z; - y = 0.5 * y; - if (y == ceil(y)) - sgn = -1; - if (z < .25) - z = sin(M_PI*z); - else - z = cos(M_PI*(0.5-z)); - /* Special case: G(1-x) = Inf; G(x) may be nonzero. */ - if (x < -170) { - if (x < -190) - return ((double)sgn*tiny*tiny); - y = one - x; /* exact: 128 < |x| < 255 */ - lg = large_gam(y); - lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */ - lg.a -= lsine.a; /* exact (opposite signs) */ - lg.b -= lsine.b; - y = -(lg.a + lg.b); - z = (y + lg.a) + lg.b; - y = __exp__D(y, z); - if (sgn < 0) y = -y; - return (y); - } - y = one-x; - if (one-y == x) - y = tgamma(y); - else /* 1-x is inexact */ - y = -x*tgamma(-x); - if (sgn < 0) y = -y; - return (M_PI / (y*z)); -} -- 2.20.1