As a component of atan2(y, x), the case of x == 1.0 is farmed out to
atan(y). The current implementation of this comparison is vulnerable
to signed integer underflow (that is, undefined behavior), and it's
performed in a somewhat more complicated way than it need be. Change
it to not be quite so cute, rather directly comparing the high/low
bits of x to the specific IEEE-754 bit pattern that encodes 1.0.
Note that while there are three different e_atan* files in the
relevant directory, only this one needs fixing. e_atan2f.c already
compares against the full bit pattern encoding 1.0f, while
e_atan2l.cuses bitwise-ands/ors/nots and so doesn't require a change.
Closes#130
Submitted by: Jeff Walden (@jswalden github PR #130)
Reviewed by: bde
MFC After: 1 month
Mainly focus on files that use BSD 2-Clause license, however the tool I
was using mis-identified many licenses so this was mostly a manual - error
prone - task.
The Software Package Data Exchange (SPDX) group provides a specification
to make it easier for automated tools to detect and summarize well known
opensource licenses. We are gradually adopting the specification, noting
that the tags are considered only advisory and do not, in any way,
superceed or replace the license texts.
LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html
Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'". Avoid this by
referring to proxy values instead.
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
LDBL_MAX is broken on i386:
https://lists.freebsd.org/pipermail/freebsd-numerics/2012-September/000288.html
Gcc has produced +Infinity for LDBL_MAX on i386 and amd64 with -m32
for some time, and newer versions of gcc are now warning that the
"floating constant exceeds range of 'long double'". Avoid this by
referring to half the value of LDBL_MAX instead.
Reviewed by: bde
Approved by: markj (mentor)
Sponsored by: Dell EMC Isilon
The primary benefit of these functions is that argument
reduction is done once instead of twice in independent
calls to sin() and cos().
* lib/msun/Makefile:
. Add s_sincos[fl].c to the build.
. Add sincos.3 documentation.
. Add appropriate MLINKS.
* lib/msun/Symbol.map:
. Expose sincos[fl] symbols in dynamic libm.so.
* lib/msun/man/sincos.3:
. Documentation for sincos[fl].
* lib/msun/src/k_sincos.h:
. Kernel for sincos() function. This merges the individual kernels
for sin() and cos(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosf.h:
. Kernel for sincosf() function. This merges the individual kernels
for sinf() and cosf(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/k_sincosl.h:
. Kernel for sincosl() function. This merges the individual kernels
for sinl() and cosl(). The merger offered an opportunity to re-arrange
the individual kernels for better performance.
* lib/msun/src/math.h:
. Add prototytpes for sincos[fl]().
* lib/msun/src/math_private.h:
. Add RETURNV macros. This is needed to reset fpsetprec on I386
hardware for a function with type void.
* lib/msun/src/s_sincos.c:
. Implementation of sincos() where sin() and cos() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosf.c:
. Implementation of sincosf() where sinf() and cosf() were merged into
one routine and possibly re-arranged for better performance.
* lib/msun/src/s_sincosl.c:
. Implementation of sincosl() where sinl() and cosl() were merged into
one routine and possibly re-arranged for better performance.
PR: 215977, 218300
Submitted by: Steven G. Kargl <sgk@troutmask.apl.washington.edu>
MFC after: 1 month
Differential Revision: https://reviews.freebsd.org/D10765
an inexact floating point exception. The variable cannot be eliminated,
unfortunately, otherwise the desired addition triggering the exception
will be emitted neither by clang, nor by gcc.
Reviewed by: Steve Kargl, bde
MFC after: 3 days
s_{fabs,fmax,logb,scalb}{,f,l}.c may be built elsewhere with a higher
WARNS setting.
Reviewed by: ed
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D8061
- Fix a case where NaNs were not mixed correctly and signalling NaNs were
not converted to quiet NaNs.
- Eliminate two negations from ctan(z).
In collaboration with: bde
but must still satisfy csinh(conj(z)) == conj(csinh(z)) and csinh(-z) ==
-csinh(z). This allows eliminating two negations from csin(z).
In collaboration with: bde
and qone[f]() were marked as __inline, but their forward
declarations were not updated. Fix the forward declarations
to match the actual function declarations.
Requested by: bde
should raise a divide-by-zero floating point exception for x = +-0
and an invalid floating point exception for x < 0 including x = -Inf.
Update the code to raise the exception and update the documentation
with hopefully better description of the behavior.
Reviewed by: bde (code only)
values for the different invervals were not converted correctly.
Adjust the threshold values to values, which should agree with the
comments.
Reported by: cognet (j1f only)
Discussed with: pfg, bde
Reviewed by: bde
Drop an unnecessary check in some calculations. The check
would have Coverity falsely conclude that variables could
be left undefined.
Discussed with: kargl, bde
Reviewed by: bde
within [INT_MIN, INT_MAX] where the magnitude of the lower
and upper bounds are sufficiently large to span the range of
scalbn[fl].
While here, remove the GNU style bug in the function declarations.
Reviewed by: bde, pfg
The changes unrelated to the bug in r277948 made
the code very difficult to understand to both
coverity and regular humans so take a step back
to something that is much easier to understand
for both and follows better the original code.
CID: 1267992, 1267993, 1267994
Discussed with: kargl
This fixes evaluation of exceptional values in scalblnl().
While here, simplify the code as suggested by Bruce Evans.
Reported by: clang static analyzer
MFC after: 1 week
The C11 standard introduced a set of macros (CMPLX, CMPLXF, CMPLXL) that
can be used to construct complex numbers from a pair of real and
imaginary numbers. Unfortunately, they require some compiler support,
which is why we only define them for Clang and GCC>=4.7.
The cpack() function in libm performs the same task as CMPLX(), but
cannot be used to generate compile-time constants. This means that all
invocations of cpack() can safely be replaced by C11's CMPLX(). To keep
the code building with GCC 4.2, provide copies of CMPLX() that can at
least be used to generate run-time complex numbers.
This makes it easier to build some of the functions outside of libm.
for |x| small.
While here, remove the explicit cast of 0.25 to float. Replace
a multiplication involving 0.25 by a division using an integer
constant 4. Make a similar change in j0() to minimize the diff.
Suggested by: bde
It is automatically set when -fPIC is passed to the compiler.
Reviewed by: dim, kib
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D1179
lgamma(x) = -log(x) - log(1+x) + x*(1-g) + x**2*P(x) with g = 0.57...
being the Euler constant and P(x) a polynomial. Substitution of small
into the RHS shows that the last 3 terms are negligible in comparison to
the leading term. The choice of 3 may be conservative.
The value large=2**(p+3) is detemined from Stirling's approximation
lgamma(x) = x*(log(x)-1) - log(x)/2 + log(2*pi)/2 + P(1/x)/x
Again, substitution of large into the RHS reveals the last 3 terms
are negligible in comparison to the leading term.
Move the x=+-0 special case into the |x|<small block.
In the ld80 and ld128 implementaion, use fdlibm compatible comparisons
involving ix, lx, and llx. This replaces several floating point
comparisons (some involving fabsl()) and also fixes the special cases
x=1 and x=2.
While here
. Remove unnecessary parentheses.
. Fix/improve comments due to the above changes.
. Fix nearby whitespace.
* src/e_lgamma_r.c:
. Sort declaration.
. Remove unneeded explicit cast for type conversion.
. Replace a double literal constant by an integer literal constant.
* src/e_lgammaf_r.c:
. Sort declaration.
* ld128/e_lgammal_r.c:
. Replace a long double literal constant by a double literal constant.
* ld80/e_lgammal_r.c:
. Remove unused '#include float.h'
. Replace a long double literal constant by a double literal constant.
Requested by: bde
. Hook e_lgammal[_r].c to the build.
. Create man page links for lgammal[-r].3.
* Symbol.map:
. Sort lgammal to its rightful place.
. Add FBSD_1.4 section for the new lgamal_r symbol.
* ld128/e_lgammal_r.c:
. 128-bit implementataion of lgammal_r().
* ld80/e_lgammal_r.c:
. Intel 80-bit format implementation of lgammal_r().
* src/e_lgamma.c:
. Expose lgammal as a weak reference to lgamma for platforms
where long double is mapped to double.
* src/e_lgamma_r.c:
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
. Expose lgammal_r as a weak reference to lgamma_r for platforms
where long double is mapped to double.
* src/e_lgammaf_r.c:
. Fixed the Cygnus Support conversion of e_lgamma_r.c to float.
This includes the generation of new polynomial and rational
approximations with fewer terms. For each approximation, include
a comment on an estimate of the accuracy over the relevant domain.
. Use integer literal constants instead of real literal constants.
Let compiler(s) do the job of conversion to the appropriate type.
This allows the removal of several explicit casts of double values
to float.
* src/e_lgammal.c:
. Wrapper for lgammal() about lgammal_r().
* src/imprecise.c:
. Remove the lgamma.
* src/math.h:
. Add a prototype for lgammal_r().
* man/lgamma.3:
. Document the new functions.
Reviewed by: bde
inf and raises the divided-by-zero exception. Compilers
constant fold one/zero to inf but do not raise the exception.
Introduce a volatile vzero to prevent the constant folding.
Move the declaration of zero into the main declaration block.
While here, fix a nearby disordering of 'lx,ix'
Discussed with: bde
Remove parentheses in a return statement to be consistent with the
rest of the file.
Rename sin_pi() in the float version to sin_pif().
Remove large comment that precedes sin_pif(). The comment
duplicates a comment in e_lgamma_r.c where the algorithm
is documented.
Requested by: bde
sin_pi(x) is only called for x < 0 and |x| < 2**(p-1) where p is
the precision of x. The new argument reduction is an optimization
compared to the old code, and it removes a chunk of dead code.
Accuracy tests in the intervals (-21,-20), (-20,-19), ... (-1,0)
show no differences between the old and new code.
Obtained from: bde
