Attachment 193340 Details for Bug 216863 – Update with move of s_clogl.c to LDBL_PREC != 53 section of the Makefile.

[patch] Update with move of s_clogl.c to LDBL_PREC != 53 section of the Makefile.

1.patch (text/plain), 24.95 KB, created by Konstantin Belousov on 2018-05-12 20:41:33 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Konstantin Belousov

Created: 2018-05-12 20:41:33 UTC

Size: 24.95 KB

patch

obsolete

>diff --git a/include/complex.h b/include/complex.h
>index b502613c6ac..0f7228fad86 100644
>--- a/include/complex.h
>+++ b/include/complex.h
>@@ -101,6 +101,10 @@ float complex	cexpf(float complex);
> double		cimag(double complex) __pure2;
> float		cimagf(float complex) __pure2;
> long double	cimagl(long double complex) __pure2;
>+double complex	clog(double complex);
>+float complex	clogf(float complex);
>+long double complex
>+		clogl(long double complex);
> double complex	conj(double complex) __pure2;
> float complex	conjf(float complex) __pure2;
> long double complex
>diff --git a/lib/libc/softfloat/bits64/softfloat-macros b/lib/libc/softfloat/bits64/softfloat-macros
>index f4647def2bd..7fedb0c348d 100644
>--- a/lib/libc/softfloat/bits64/softfloat-macros
>+++ b/lib/libc/softfloat/bits64/softfloat-macros
>@@ -157,7 +157,7 @@ INLINE void
>         z0 = a0>>count;
>     }
>     else {
>-        z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;
>+        z1 = ( count < 128 ) ? ( a0>>( count & 63 ) ) : 0;
>         z0 = 0;
>     }
>     *z1Ptr = z1;
>diff --git a/lib/msun/Makefile b/lib/msun/Makefile
>index ed9f96e61ae..7f4283c01ed 100644
>--- a/lib/msun/Makefile
>+++ b/lib/msun/Makefile
>@@ -57,7 +57,7 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
> 	k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \
> 	k_tan.c k_tanf.c \
> 	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \
>-	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \
>+	s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c s_clog.c s_clogf.c \
> 	s_copysign.c s_copysignf.c s_cos.c s_cosf.c \
> 	s_csqrt.c s_csqrtf.c s_erf.c s_erff.c \
> 	s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c \
>@@ -101,7 +101,8 @@ COMMON_SRCS+=	catrigl.c \
> 	e_lgammal.c e_lgammal_r.c \
> 	e_remainderl.c e_sinhl.c e_sqrtl.c \
> 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
>-	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
>+	s_asinhl.c s_atanl.c s_cbrtl.c s_ceill.c \
>+	s_clogl.c s_cosl.c s_cprojl.c \
> 	s_csqrtl.c s_erfl.c s_exp2l.c s_expl.c s_floorl.c s_fmal.c \
> 	s_fmaxl.c s_fminl.c s_frexpl.c s_logbl.c s_logl.c s_nanl.c \
> 	s_nextafterl.c s_nexttoward.c s_remquol.c s_rintl.c s_roundl.c \
>@@ -133,7 +134,8 @@ INCS+=	fenv.h math.h
> 
> MAN=	acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
> 	ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
>-	cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
>+	cimag.3 clog.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 \
>+	exp.3 fabs.3 fdim.3 \
> 	feclearexcept.3 feenableexcept.3 fegetenv.3 \
> 	fegetround.3 fenv.3 floor.3 \
> 	fma.3 fmax.3 fmod.3 hypot.3 ieee.3 ieee_test.3 ilogb.3 j0.3 \
>@@ -166,6 +168,7 @@ MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
> 	cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \
> 	cimag.3 cproj.3 cimag.3 cprojf.3 cimag.3 cprojl.3 \
> 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
>+MLINKS+=clog.3 clogf.3 clog.3 clogl.3
> MLINKS+=copysign.3 copysignf.3 copysign.3 copysignl.3
> MLINKS+=cos.3 cosf.3 cos.3 cosl.3
> MLINKS+=cosh.3 coshf.3 cosh.3 coshl.3
>diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
>index f7d9e03d0f8..a9837f6c4d4 100644
>--- a/lib/msun/Symbol.map
>+++ b/lib/msun/Symbol.map
>@@ -294,6 +294,9 @@ FBSD_1.5 {
> 	casinl;
> 	catanl;
> 	catanhl;
>+	clog;
>+	clogf;
>+	clogl;
> 	sincos;
> 	sincosf;
> 	sincosl;
>diff --git a/lib/msun/man/clog.3 b/lib/msun/man/clog.3
>new file mode 100644
>index 00000000000..0d23a845fb6
>--- /dev/null
>+++ b/lib/msun/man/clog.3
>@@ -0,0 +1,103 @@
>+.\" Copyright (c) 2017 Steven G. Kargl <kargl@FreeBSD.org>
>+.\" 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.
>+.\"
>+.\" THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
>+.\"
>+.\" $FreeBSD$
>+.\"
>+.Dd February 13, 2017
>+.Dt CLOG 3
>+.Os
>+.Sh NAME
>+.Nm clog ,
>+.Nm clogf ,
>+and
>+.Nm clogl
>+.Nd complex natural logrithm functions
>+.Sh LIBRARY
>+.Lb libm
>+.Sh SYNOPSIS
>+.In complex.h
>+.Ft double complex
>+.Fn clog "double complex z"
>+.Ft float complex
>+.Fn clogf "float complex z"
>+.Ft long double complex
>+.Fn clogl "long double complex z"
>+.Sh DESCRIPTION
>+The
>+.Fn clog ,
>+.Fn clogf ,
>+and 
>+.Fn clogl
>+functions compute the complex natural logrithm of
>+.Fa z .
>+with a branch cut along the negative real axis .
>+.Sh RETURN VALUES
>+The
>+.Fn clog
>+function returns the complex natural logarithm value, in the
>+range of a strip mathematically unbounded along the real axis and in
>+the interval [-I* \*(Pi , +I* \*(Pi ] along the imaginary axis.
>+The function satisfies the relationship: 
>+.Fo clog
>+.Fn conj "z" Fc
>+=
>+.Fo conj
>+.Fn clog "z" Fc .
>+.Pp
>+
>+.\" Table is formatted for an 80-column xterm.
>+.Bl -column ".Sy +\*(If + I*\*(Na" ".Sy Return value" ".Sy Divide-by-zero exception"
>+.It Sy Argument          Ta Sy Return value Ta Sy Comment
>+.It -0 + I*0             Ta -\*(If + I*\*(Pi    Ta Divide-by-zero exception
>+.It                      Ta                     Ta raised
>+.It +0 + I*0             Ta -\*(If + I*0        Ta Divide by zero exception
>+.It                      Ta                     Ta raised
>+.It x + I*\*(If          Ta +\*(If + I*\*(Pi/2  Ta For finite x
>+.It x + I*\*(Na          Ta  \*(Na + I*\*(Na    Ta Optionally raises invalid
>+.It                      Ta                     Ta floating-point exception
>+.It                      Ta                     Ta for finite x
>+.It -\*(If + I*y         Ta +\*(If + I*\*(Pi    Ta For finite positive-signed y
>+.It +\*(If + I*y         Ta +\*(If + I*0        Ta For finite positive-signed y
>+.It -\*(If + I*\*(If     Ta +\*(If + I*3\*(Pi/4
>+.It +\*(If + I*\*(If     Ta +\*(If + I*\*(Pi/4
>+.It \*(Pm\*(If + I*\*(Na Ta +\*(If + I*\*(Na
>+.It \*(Na + I*y          Ta \*(Na + I*\*(Na    Ta Optionally raises invalid
>+.It                      Ta                    Ta floating-point exception
>+.It                      Ta                    Ta for finite y
>+.It \*(Na + I*\*(If      Ta +\*(If + I*\*(Na
>+.It \*(Na + I*\*(Na      Ta \*(Na + I*\*(Na 
>+.El
>+
>+.Sh SEE ALSO
>+.Xr complex 3 ,
>+.Xr log 3 ,
>+.Xr math 3
>+.Sh STANDARDS
>+The
>+.Fn clog ,
>+.Fn cexpf ,
>+and
>+.Fn clogl
>+functions conform to
>+.St -isoC-99 .
>diff --git a/lib/msun/man/complex.3 b/lib/msun/man/complex.3
>index d2864944cf0..42f6e8d3d9b 100644
>--- a/lib/msun/man/complex.3
>+++ b/lib/msun/man/complex.3
>@@ -24,7 +24,7 @@
> .\"
> .\" $FreeBSD$
> .\"
>-.Dd October 17, 2011
>+.Dd February 13, 2017
> .Dt COMPLEX 3
> .Os
> .Sh NAME
>@@ -77,6 +77,10 @@ csqrt	complex square root
> .Cl
> cexp	exponential base e
> .El
>+.Ss Natural logrithm Function
>+.Cl
>+clog	natural logrithm
>+.El
> .\" Section 7.3.9 of ISO C99 standard
> .Ss Manipulation Functions
> .Cl
>@@ -117,8 +121,6 @@ The
> functions described here conform to
> .St -isoC-99 .
> .Sh BUGS
>-The logarithmic functions
>-.Fn clog
>-and the power functions
>+The power functions
> .Fn cpow
> are not implemented.
>diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h
>index 74879027dba..7f81381a242 100644
>--- a/lib/msun/src/math_private.h
>+++ b/lib/msun/src/math_private.h
>@@ -294,8 +294,9 @@ do {								\
> 
> /* Support switching the mode to FP_PE if necessary. */
> #if defined(__i386__) && !defined(NO_FPSETPREC)
>-#define	ENTERI()				\
>-	long double __retval;			\
>+#define	ENTERI() ENTERIT(long double)
>+#define	ENTERIT(returntype)			\
>+	returntype __retval;			\
> 	fp_prec_t __oprec;			\
> 						\
> 	if ((__oprec = fpgetprec()) != FP_PE)	\
>@@ -318,6 +319,7 @@ do {								\
> } while (0)
> #else
> #define	ENTERI()
>+#define	ENTERIT(x)
> #define	RETURNI(x)	RETURNF(x)
> #define	ENTERV()
> #define	RETURNV()	return
>diff --git a/lib/msun/src/s_clog.c b/lib/msun/src/s_clog.c
>new file mode 100644
>index 00000000000..8150fa7f8bc
>--- /dev/null
>+++ b/lib/msun/src/s_clog.c
>@@ -0,0 +1,155 @@
>+/*-
>+ * Copyright (c) 2013 Bruce D. Evans
>+ * 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 unmodified, 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 AUTHOR ``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 AUTHOR 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 <sys/cdefs.h>
>+__FBSDID("$FreeBSD$");
>+
>+#include <complex.h>
>+#include <float.h>
>+
>+#include "fpmath.h"
>+#include "math.h"
>+#include "math_private.h"
>+
>+#define	MANT_DIG	DBL_MANT_DIG
>+#define	MAX_EXP		DBL_MAX_EXP
>+#define	MIN_EXP		DBL_MIN_EXP
>+
>+static const double
>+ln2_hi = 6.9314718055829871e-1,		/*  0x162e42fefa0000.0p-53 */
>+ln2_lo = 1.6465949582897082e-12;	/*  0x1cf79abc9e3b3a.0p-92 */
>+
>+double complex
>+clog(double complex z)
>+{
>+	double_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
>+	double x, y, v;
>+	uint32_t hax, hay;
>+	int kx, ky;
>+
>+	x = creal(z);
>+	y = cimag(z);
>+	v = atan2(y, x);
>+
>+	ax = fabs(x);
>+	ay = fabs(y);
>+	if (ax < ay) {
>+		t = ax;
>+		ax = ay;
>+		ay = t;
>+	}
>+
>+	GET_HIGH_WORD(hax, ax);
>+	kx = (hax >> 20) - 1023;
>+	GET_HIGH_WORD(hay, ay);
>+	ky = (hay >> 20) - 1023;
>+
>+	/* Handle NaNs and Infs using the general formula. */
>+	if (kx == MAX_EXP || ky == MAX_EXP)
>+		return (CMPLX(log(hypot(x, y)), v));
>+
>+	/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
>+	if (ax == 1) {
>+		if (ky < (MIN_EXP - 1) / 2)
>+			return (CMPLX((ay / 2) * ay, v));
>+		return (CMPLX(log1p(ay * ay) / 2, v));
>+	}
>+
>+	/* Avoid underflow when ax is not small.  Also handle zero args. */
>+	if (kx - ky > MANT_DIG || ay == 0)
>+		return (CMPLX(log(ax), v));
>+
>+	/* Avoid overflow. */
>+	if (kx >= MAX_EXP - 1)
>+		return (CMPLX(log(hypot(x * 0x1p-1022, y * 0x1p-1022)) +
>+		    (MAX_EXP - 2) * ln2_lo + (MAX_EXP - 2) * ln2_hi, v));
>+	if (kx >= (MAX_EXP - 1) / 2)
>+		return (CMPLX(log(hypot(x, y)), v));
>+
>+	/* Reduce inaccuracies and avoid underflow when ax is denormal. */
>+	if (kx <= MIN_EXP - 2)
>+		return (CMPLX(log(hypot(x * 0x1p1023, y * 0x1p1023)) +
>+		    (MIN_EXP - 2) * ln2_lo + (MIN_EXP - 2) * ln2_hi, v));
>+
>+	/* Avoid remaining underflows (when ax is small but not denormal). */
>+	if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
>+		return (CMPLX(log(hypot(x, y)), v));
>+
>+	/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
>+	t = (double)(ax * (0x1p27 + 1));
>+	axh = (double)(ax - t) + t;
>+	axl = ax - axh;
>+	ax2h = ax * ax;
>+	ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
>+	t = (double)(ay * (0x1p27 + 1));
>+	ayh = (double)(ay - t) + t;
>+	ayl = ay - ayh;
>+	ay2h = ay * ay;
>+	ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
>+
>+	/*
>+	 * When log(|z|) is far from 1, accuracy in calculating the sum
>+	 * of the squares is not very important since log() reduces
>+	 * inaccuracies.  We depended on this to use the general
>+	 * formula when log(|z|) is very far from 1.  When log(|z|) is
>+	 * moderately far from 1, we go through the extra-precision
>+	 * calculations to reduce branches and gain a little accuracy.
>+	 *
>+	 * When |z| is near 1, we subtract 1 and use log1p() and don't
>+	 * leave it to log() to subtract 1, since we gain at least 1 bit
>+	 * of accuracy in this way.
>+	 *
>+	 * When |z| is very near 1, subtracting 1 can cancel almost
>+	 * 3*MANT_DIG bits.  We arrange that subtracting 1 is exact in
>+	 * doubled precision, and then do the rest of the calculation
>+	 * in sloppy doubled precision.  Although large cancellations
>+	 * often lose lots of accuracy, here the final result is exact
>+	 * in doubled precision if the large calculation occurs (because
>+	 * then it is exact in tripled precision and the cancellation
>+	 * removes enough bits to fit in doubled precision).  Thus the
>+	 * result is accurate in sloppy doubled precision, and the only
>+	 * significant loss of accuracy is when it is summed and passed
>+	 * to log1p().
>+	 */
>+	sh = ax2h;
>+	sl = ay2h;
>+	_2sumF(sh, sl);
>+	if (sh < 0.5 || sh >= 3)
>+		return (CMPLX(log(ay2l + ax2l + sl + sh) / 2, v));
>+	sh -= 1;
>+	_2sum(sh, sl);
>+	_2sum(ax2l, ay2l);
>+	/* Briggs-Kahan algorithm (except we discard the final low term): */
>+	_2sum(sh, ax2l);
>+	_2sum(sl, ay2l);
>+	t = ax2l + sl;
>+	_2sumF(sh, t);
>+	return (CMPLX(log1p(ay2l + t + sh) / 2, v));
>+}
>+
>+#if (LDBL_MANT_DIG == 53)
>+__weak_reference(clog, clogl);
>+#endif
>diff --git a/lib/msun/src/s_clogf.c b/lib/msun/src/s_clogf.c
>new file mode 100644
>index 00000000000..19a94453ea3
>--- /dev/null
>+++ b/lib/msun/src/s_clogf.c
>@@ -0,0 +1,151 @@
>+/*-
>+ * Copyright (c) 2013 Bruce D. Evans
>+ * 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 unmodified, 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 AUTHOR ``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 AUTHOR 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 <sys/cdefs.h>
>+__FBSDID("$FreeBSD$");
>+
>+#include <complex.h>
>+#include <float.h>
>+
>+#include "fpmath.h"
>+#include "math.h"
>+#include "math_private.h"
>+
>+#define	MANT_DIG	FLT_MANT_DIG
>+#define	MAX_EXP		FLT_MAX_EXP
>+#define	MIN_EXP		FLT_MIN_EXP
>+
>+static const float
>+ln2f_hi =  6.9314575195e-1,		/*  0xb17200.0p-24 */
>+ln2f_lo =  1.4286067653e-6;		/*  0xbfbe8e.0p-43 */
>+
>+float complex
>+clogf(float complex z)
>+{
>+	float_t ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl, sh, sl, t;
>+	float x, y, v;
>+	uint32_t hax, hay;
>+	int kx, ky;
>+
>+	x = crealf(z);
>+	y = cimagf(z);
>+	v = atan2f(y, x);
>+
>+	ax = fabsf(x);
>+	ay = fabsf(y);
>+	if (ax < ay) {
>+		t = ax;
>+		ax = ay;
>+		ay = t;
>+	}
>+
>+	GET_FLOAT_WORD(hax, ax);
>+	kx = (hax >> 23) - 127;
>+	GET_FLOAT_WORD(hay, ay);
>+	ky = (hay >> 23) - 127;
>+
>+	/* Handle NaNs and Infs using the general formula. */
>+	if (kx == MAX_EXP || ky == MAX_EXP)
>+		return (CMPLXF(logf(hypotf(x, y)), v));
>+
>+	/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
>+	if (hax == 0x3f800000) {
>+		if (ky < (MIN_EXP - 1) / 2)
>+			return (CMPLXF((ay / 2) * ay, v));
>+		return (CMPLXF(log1pf(ay * ay) / 2, v));
>+	}
>+
>+	/* Avoid underflow when ax is not small.  Also handle zero args. */
>+	if (kx - ky > MANT_DIG || hay == 0)
>+		return (CMPLXF(logf(ax), v));
>+
>+	/* Avoid overflow. */
>+	if (kx >= MAX_EXP - 1)
>+		return (CMPLXF(logf(hypotf(x * 0x1p-126F, y * 0x1p-126F)) +
>+		    (MAX_EXP - 2) * ln2f_lo + (MAX_EXP - 2) * ln2f_hi, v));
>+	if (kx >= (MAX_EXP - 1) / 2)
>+		return (CMPLXF(logf(hypotf(x, y)), v));
>+
>+	/* Reduce inaccuracies and avoid underflow when ax is denormal. */
>+	if (kx <= MIN_EXP - 2)
>+		return (CMPLXF(logf(hypotf(x * 0x1p127F, y * 0x1p127F)) +
>+		    (MIN_EXP - 2) * ln2f_lo + (MIN_EXP - 2) * ln2f_hi, v));
>+
>+	/* Avoid remaining underflows (when ax is small but not denormal). */
>+	if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
>+		return (CMPLXF(logf(hypotf(x, y)), v));
>+
>+	/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
>+	t = (float)(ax * (0x1p12F + 1));
>+	axh = (float)(ax - t) + t;
>+	axl = ax - axh;
>+	ax2h = ax * ax;
>+	ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
>+	t = (float)(ay * (0x1p12F + 1));
>+	ayh = (float)(ay - t) + t;
>+	ayl = ay - ayh;
>+	ay2h = ay * ay;
>+	ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
>+
>+	/*
>+	 * When log(|z|) is far from 1, accuracy in calculating the sum
>+	 * of the squares is not very important since log() reduces
>+	 * inaccuracies.  We depended on this to use the general
>+	 * formula when log(|z|) is very far from 1.  When log(|z|) is
>+	 * moderately far from 1, we go through the extra-precision
>+	 * calculations to reduce branches and gain a little accuracy.
>+	 *
>+	 * When |z| is near 1, we subtract 1 and use log1p() and don't
>+	 * leave it to log() to subtract 1, since we gain at least 1 bit
>+	 * of accuracy in this way.
>+	 *
>+	 * When |z| is very near 1, subtracting 1 can cancel almost
>+	 * 3*MANT_DIG bits.  We arrange that subtracting 1 is exact in
>+	 * doubled precision, and then do the rest of the calculation
>+	 * in sloppy doubled precision.  Although large cancellations
>+	 * often lose lots of accuracy, here the final result is exact
>+	 * in doubled precision if the large calculation occurs (because
>+	 * then it is exact in tripled precision and the cancellation
>+	 * removes enough bits to fit in doubled precision).  Thus the
>+	 * result is accurate in sloppy doubled precision, and the only
>+	 * significant loss of accuracy is when it is summed and passed
>+	 * to log1p().
>+	 */
>+	sh = ax2h;
>+	sl = ay2h;
>+	_2sumF(sh, sl);
>+	if (sh < 0.5F || sh >= 3)
>+		return (CMPLXF(logf(ay2l + ax2l + sl + sh) / 2, v));
>+	sh -= 1;
>+	_2sum(sh, sl);
>+	_2sum(ax2l, ay2l);
>+	/* Briggs-Kahan algorithm (except we discard the final low term): */
>+	_2sum(sh, ax2l);
>+	_2sum(sl, ay2l);
>+	t = ax2l + sl;
>+	_2sumF(sh, t);
>+	return (CMPLXF(log1pf(ay2l + t + sh) / 2, v));
>+}
>diff --git a/lib/msun/src/s_clogl.c b/lib/msun/src/s_clogl.c
>new file mode 100644
>index 00000000000..e59a1377afd
>--- /dev/null
>+++ b/lib/msun/src/s_clogl.c
>@@ -0,0 +1,168 @@
>+/*-
>+ * Copyright (c) 2013 Bruce D. Evans
>+ * 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 unmodified, 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 AUTHOR ``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 AUTHOR 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 <sys/cdefs.h>
>+__FBSDID("$FreeBSD$");
>+
>+#include <complex.h>
>+#include <float.h>
>+#ifdef __i386__
>+#include <ieeefp.h>
>+#endif
>+
>+#include "fpmath.h"
>+#include "math.h"
>+#include "math_private.h"
>+
>+#define	MANT_DIG	LDBL_MANT_DIG
>+#define	MAX_EXP		LDBL_MAX_EXP
>+#define	MIN_EXP		LDBL_MIN_EXP
>+
>+static const double
>+ln2_hi = 6.9314718055829871e-1;		/*  0x162e42fefa0000.0p-53 */
>+
>+#if LDBL_MANT_DIG == 64
>+#define	MULT_REDUX	0x1p32		/* exponent MANT_DIG / 2 rounded up */
>+static const double
>+ln2l_lo = 1.6465949582897082e-12;	/*  0x1cf79abc9e3b3a.0p-92 */
>+#elif LDBL_MANT_DIG == 113
>+#define	MULT_REDUX	0x1p57
>+static const long double
>+ln2l_lo = 1.64659495828970812809844307550013433e-12L;	/*  0x1cf79abc9e3b39803f2f6af40f343.0p-152L */
>+#else
>+#error "Unsupported long double format"
>+#endif
>+
>+long double complex
>+clogl(long double complex z)
>+{
>+	long double ax, ax2h, ax2l, axh, axl, ay, ay2h, ay2l, ayh, ayl;
>+	long double sh, sl, t;
>+	long double x, y, v;
>+	uint16_t hax, hay;
>+	int kx, ky;
>+
>+	ENTERIT(long double complex);
>+
>+	x = creall(z);
>+	y = cimagl(z);
>+	v = atan2l(y, x);
>+
>+	ax = fabsl(x);
>+	ay = fabsl(y);
>+	if (ax < ay) {
>+		t = ax;
>+		ax = ay;
>+		ay = t;
>+	}
>+
>+	GET_LDBL_EXPSIGN(hax, ax);
>+	kx = hax - 16383;
>+	GET_LDBL_EXPSIGN(hay, ay);
>+	ky = hay - 16383;
>+
>+	/* Handle NaNs and Infs using the general formula. */
>+	if (kx == MAX_EXP || ky == MAX_EXP)
>+		RETURNI(CMPLXL(logl(hypotl(x, y)), v));
>+
>+	/* Avoid spurious underflow, and reduce inaccuracies when ax is 1. */
>+	if (ax == 1) {
>+		if (ky < (MIN_EXP - 1) / 2)
>+			RETURNI(CMPLXL((ay / 2) * ay, v));
>+		RETURNI(CMPLXL(log1pl(ay * ay) / 2, v));
>+	}
>+
>+	/* Avoid underflow when ax is not small.  Also handle zero args. */
>+	if (kx - ky > MANT_DIG || ay == 0)
>+		RETURNI(CMPLXL(logl(ax), v));
>+
>+	/* Avoid overflow. */
>+	if (kx >= MAX_EXP - 1)
>+		RETURNI(CMPLXL(logl(hypotl(x * 0x1p-16382L, y * 0x1p-16382L)) +
>+		    (MAX_EXP - 2) * ln2l_lo + (MAX_EXP - 2) * ln2_hi, v));
>+	if (kx >= (MAX_EXP - 1) / 2)
>+		RETURNI(CMPLXL(logl(hypotl(x, y)), v));
>+
>+	/* Reduce inaccuracies and avoid underflow when ax is denormal. */
>+	if (kx <= MIN_EXP - 2)
>+		RETURNI(CMPLXL(logl(hypotl(x * 0x1p16383L, y * 0x1p16383L)) +
>+		    (MIN_EXP - 2) * ln2l_lo + (MIN_EXP - 2) * ln2_hi, v));
>+
>+	/* Avoid remaining underflows (when ax is small but not denormal). */
>+	if (ky < (MIN_EXP - 1) / 2 + MANT_DIG)
>+		RETURNI(CMPLXL(logl(hypotl(x, y)), v));
>+
>+	/* Calculate ax*ax and ay*ay exactly using Dekker's algorithm. */
>+	t = (long double)(ax * (MULT_REDUX + 1));
>+	axh = (long double)(ax - t) + t;
>+	axl = ax - axh;
>+	ax2h = ax * ax;
>+	ax2l = axh * axh - ax2h + 2 * axh * axl + axl * axl;
>+	t = (long double)(ay * (MULT_REDUX + 1));
>+	ayh = (long double)(ay - t) + t;
>+	ayl = ay - ayh;
>+	ay2h = ay * ay;
>+	ay2l = ayh * ayh - ay2h + 2 * ayh * ayl + ayl * ayl;
>+
>+	/*
>+	 * When log(|z|) is far from 1, accuracy in calculating the sum
>+	 * of the squares is not very important since log() reduces
>+	 * inaccuracies.  We depended on this to use the general
>+	 * formula when log(|z|) is very far from 1.  When log(|z|) is
>+	 * moderately far from 1, we go through the extra-precision
>+	 * calculations to reduce branches and gain a little accuracy.
>+	 *
>+	 * When |z| is near 1, we subtract 1 and use log1p() and don't
>+	 * leave it to log() to subtract 1, since we gain at least 1 bit
>+	 * of accuracy in this way.
>+	 *
>+	 * When |z| is very near 1, subtracting 1 can cancel almost
>+	 * 3*MANT_DIG bits.  We arrange that subtracting 1 is exact in
>+	 * doubled precision, and then do the rest of the calculation
>+	 * in sloppy doubled precision.  Although large cancellations
>+	 * often lose lots of accuracy, here the final result is exact
>+	 * in doubled precision if the large calculation occurs (because
>+	 * then it is exact in tripled precision and the cancellation
>+	 * removes enough bits to fit in doubled precision).  Thus the
>+	 * result is accurate in sloppy doubled precision, and the only
>+	 * significant loss of accuracy is when it is summed and passed
>+	 * to log1p().
>+	 */
>+	sh = ax2h;
>+	sl = ay2h;
>+	_2sumF(sh, sl);
>+	if (sh < 0.5 || sh >= 3)
>+		RETURNI(CMPLXL(logl(ay2l + ax2l + sl + sh) / 2, v));
>+	sh -= 1;
>+	_2sum(sh, sl);
>+	_2sum(ax2l, ay2l);
>+	/* Briggs-Kahan algorithm (except we discard the final low term): */
>+	_2sum(sh, ax2l);
>+	_2sum(sl, ay2l);
>+	t = ax2l + sl;
>+	_2sumF(sh, t);
>+	RETURNI(CMPLXL(log1pl(ay2l + t + sh) / 2, v));
>+}

Actions: View | Diff

Attachments on bug 216863: 180028 | 187872 | 193337 | 193340