d/freebsd-skq
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add implementations for clog(3), clogf(3), and clog(3).

Browse Source 
PR:	216863
Submitted by:	bde, Steven G. Kargl <sgk@troutmask.apl.washington.edu>
MFC after:	2 weeks
This commit is contained in:
Konstantin Belousov 2018-05-13 09:54:34 +00:00
parent 2ebc882927
commit 0c0288a218
10 changed files with 601 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
include/complex.h
View File

@ -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

2
lib/libc/softfloat/bits64/softfloat-macros
View File

@ -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;

9
lib/msun/Makefile
View File

@ -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

3
lib/msun/Symbol.map
View File

@ -294,6 +294,9 @@ FBSD_1.5 {
casinl;
catanl;
catanhl;
clog;
clogf;
clogl;
sincos;
sincosf;
sincosl;

103
lib/msun/man/clog.3 Normal file
View File

@ -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 .

10
lib/msun/man/complex.3
View File

@ -24,7 +24,7 @@
.\"
.\" $FreeBSD$
.\"
.Dd October 17, 2011
.Dd May 13, 2018
.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.

6
lib/msun/src/math_private.h
View File

@ -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

155
lib/msun/src/s_clog.c Normal file
View File

@ -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

151
lib/msun/src/s_clogf.c Normal file
View File

@ -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));
}

168
lib/msun/src/s_clogl.c Normal file
View File

@ -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));
}