- Handle cases where exp(x) would overflow, but ccosh(x) ~= exp(x) / 2
shouldn't.
- Use the ccosh(x) ~= exp(x) / 2 approximation to simplify the calculation
when x is large.
Similarly for csinh(). Also fixed the return value of csinh(-Inf +- 0i).
exp(x) scaled down by some factor, and the challenge is doing this
accurately when exp(x) would overflow. This change replaces all of
the tricks we've been using with common __ldexp_exp() and
__ldexp_cexp() routines that handle all the scaling.
bde plans to improve on this further by moving the guts of exp() into
k_exp.c and handling the scaling in a more direct manner. But the
current approach is simple and adequate for now.
library," since complex.h, tgmath.h, and fenv.h are also part of the
math library. Replace the outdated sentence with some references to
the other parts.
- Rename __kernel_log() to k_log1p().
- Move some of the work that was previously done in the kernel log into
the callers. This enables further refactoring to improve accuracy or
speed, although I don't recall the details.
- Use extra precision when adding the final scaling term, which improves
accuracy.
- Describe and work around compiler problems that break some of the
multiprecision calculations.
A fix for a small bug is also included:
- Add a special case for log*(1). This is needed to ensure that log*(1) == +0
instead of -0, even when the rounding mode is FE_DOWNWARD.
Submitted by: bde
no longer "fast" on sparc64. (It really wasn't to begin with, since
the old implementation was using long doubles, and long doubles are
emulated in software on sparc64.)
round-to-nearest mode when the result, rounded to twice machine
precision, was exactly halfway between two machine-precision
values. The essence of the fix is to simulate a "sticky bit" in
the pathological cases, which is how hardware implementations
break the ties.
MFC after: 1 month
fenv.h that are currently inlined.
The definitions are provided in fenv.c via 'extern inline'
declaractions. This assumes the compiler handles 'extern inline' as
specified in C99, which has been true under FreeBSD since 8.0.
The goal is to eventually remove the 'static' keyword from the inline
definitions in fenv.h, so that non-inlined references all wind up
pointing to the same external definition like they're supposed to.
I am deferring the second step to provide a window where
newly-compiled apps will still link against old math libraries.
(This isn't supported, but there's no need to cause undue breakage.)
Reviewed by: stefanf, bde
on i386-class hardware for sinl and cosl. The hand-rolled argument
reduction have been replaced by e_rem_pio2l() implementations. To
preserve history the following commands have been executed:
svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h
mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h
svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h
mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h
The ld80 version has been tested by bde, das, and kargl over the
last few years (bde, das) and few months (kargl). An older ld128
version was tested by das. The committed version has only been
compiled tested via 'make universe'.
Approved by: das (mentor)
Obtained from: bde
with r219571 and re-enable building of cbrtl.
Implement the long double version for the cube root function, cbrtl.
The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.
Reviewed by: bde
Approved by: das
implementing accurate logarithms in different bases. This is based
on an approach bde coded up years ago.
This function should always be inlined; it will be used in only a few
places, and rudimentary tests show a 40% performance improvement in
implementations of log2() and log10() on amd64.
The kernel takes a reduced argument x and returns the same polynomial
approximation as e_log.c, but omitting the low-order term. The low-order
term is much larger than the rest of the approximation, so the caller of
the kernel function can scale it to the appropriate base in extra precision
and obtain a much more accurate answer than by using log(x)/log(b).
Explanation by Steve:
jn[f](n,x) for certain ranges of x uses downward recursion to compute
the value of the function. The recursion sequence that is generated is
proportional to the actual desired value, so a normalization step is
taken. This normalization is j0[f](x) divided by the zeroth sequence
member. As Bruce notes, near the zeros of j0[f](x) the computed value
can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only
the leading decimal digit is correct. The solution to the issue is
fairly straight forward. The zeros of j0(x) and j1(x) never coincide,
so as j0(x) approaches a zero, the normalization constant switches to
j1[f](x) divided by the 2nd sequence member. The expectation is that
j1[f](x) is a more accurately computed value.
PR: bin/144306
Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu>
Reviewed by: bde
MFC after: 7 days
and one under lib/msun/amd64. This avoids adding the identifiers to the
.text section, and moves them to the .comment section instead.
Suggested by: bde
Approved by: rpaulo (mentor)
