Attachment 180028 Details for Bug 216863 – clog, clogf, clogl implementations with docs

clog, clogf, clogl implementations with docs

clog.diff (text/plain), 24.45 KB, created by sgk on 2017-02-15 19:11:49 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: sgk

Created: 2017-02-15 19:11:49 UTC

Size: 24.45 KB

patch

obsolete

>Index: include/complex.h
>===================================================================
>--- include/complex.h	(revision 313777)
>+++ include/complex.h	(working copy)
>@@ -87,6 +87,10 @@
> 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
>Index: lib/libc/softfloat/bits64/softfloat-macros
>===================================================================
>--- lib/libc/softfloat/bits64/softfloat-macros	(revision 313777)
>+++ lib/libc/softfloat/bits64/softfloat-macros	(working copy)
>@@ -157,7 +157,7 @@
>         z0 = a0>>count;
>     }
>     else {
>-        z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;
>+        z1 = ( count < 128 ) ? ( a0>>( count & 63 ) ) : 0;
>         z0 = 0;
>     }
>     *z1Ptr = z1;
>Index: lib/msun/Makefile
>===================================================================
>--- lib/msun/Makefile	(revision 313777)
>+++ lib/msun/Makefile	(working copy)
>@@ -57,7 +57,7 @@
> 	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_clogl.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 \
>@@ -130,7 +131,8 @@
> 
> 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 \
>@@ -159,6 +164,7 @@
> 	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
>Index: lib/msun/Symbol.map
>===================================================================
>--- lib/msun/Symbol.map	(revision 313777)
>+++ lib/msun/Symbol.map	(working copy)
>@@ -294,4 +294,7 @@
> 	casinl;
> 	catanl;
> 	catanhl;
>+	clog;
>+	clogf;
>+	clogl;
> };
>Index: lib/msun/man/clog.3
>===================================================================
>--- lib/msun/man/clog.3	(nonexistent)
>+++ lib/msun/man/clog.3	(working copy)
>@@ -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 .
>Index: lib/msun/man/complex.3
>===================================================================
>--- lib/msun/man/complex.3	(revision 313777)
>+++ lib/msun/man/complex.3	(working copy)
>@@ -24,7 +24,7 @@
> .\"
> .\" $FreeBSD$
> .\"
>-.Dd October 17, 2011
>+.Dd February 13, 2017
> .Dt COMPLEX 3
> .Os
> .Sh NAME
>@@ -77,6 +77,10 @@
> .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 @@
> 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.
>Index: lib/msun/src/math_private.h
>===================================================================
>--- lib/msun/src/math_private.h	(revision 313777)
>+++ lib/msun/src/math_private.h	(working copy)
>@@ -294,8 +294,9 @@
> 
> /* 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)	\
>@@ -308,6 +309,7 @@
> } while (0)
> #else
> #define	ENTERI(x)
>+#define	ENTERIT(x)
> #define	RETURNI(x)	RETURNF(x)
> #endif
> 
>Index: lib/msun/src/s_clog.c
>===================================================================
>--- lib/msun/src/s_clog.c	(nonexistent)
>+++ lib/msun/src/s_clog.c	(working copy)
>@@ -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
>Index: lib/msun/src/s_clogf.c
>===================================================================
>--- lib/msun/src/s_clogf.c	(nonexistent)
>+++ lib/msun/src/s_clogf.c	(working copy)
>@@ -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));
>+}
>Index: lib/msun/src/s_clogl.c
>===================================================================
>--- lib/msun/src/s_clogl.c	(nonexistent)
>+++ lib/msun/src/s_clogl.c	(working copy)
>@@ -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

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