This removes the __ieee754_ prefix from a number of the math functions.
msun/src/math_private.h contains the statement that
/*
* ieee style elementary functions
*
* We rename functions here to improve other sources' diffability
* against fdlibm.
*/
#define __ieee754_sqrt sqrt
...
Here, fdlibm refers to https://netlib.org/fdlibm. It is seen from
https://netlib.org/fdlibm/readme that this prefix was used to
differentiate between different standards:
Wrapper functions will twist the result of the ieee754
function to comply to the standard specified by the value
of _LIB_VERSION
if _LIB_VERSION = _IEEE_, return the ieee754 result;
if _LIB_VERSION = _SVID_, return SVID result;
if _LIB_VERSION = _XOPEN_, return XOPEN result;
if _LIB_VERSION = _POSIX_, return POSIX/ANSI result.
(These are macros, see fdlibm.h for their definition.)
AFAICT, FreeBSD has never supported these wrappers. In addition, as C99,
principally the long double, functions were added to libm, this
convention was not maintained. Given that only 148 of 324 files under
lib/msun contain a "Copyright (C) 1993 by Sun Microsystems" statement,
the removal of the __ieee754_ prefix provides consistency across all
source files.
The last time someone compared lib/msun to fdlibm appears to be
commit 3f70824172
Author: David Schultz <das@FreeBSD.org>
Date: Fri Feb 4 18:26:06 2005 +0000
Reduce diffs against vendor source (Sun fdlibm 5.3).
The most recent fdlibm RCS string that appears in a Sun Microsystem
copyrighted file is date "95/01/18". With Oracle Corporation's
acquisition of Sun Microsystems in 2009, it is unlikely that fdlibm will
ever be updated. A search for fdlibm at https://opensource.oracle.com/
yields no hits.
Finally, OpenBSD removed the use of this prefix over 21 years ago. pSee
revision 1.6 of OpenBSD's math_private.h.
Note: this does not drop the __ieee754_ prefix from the trigonometric
argument reduction functions, e.g., __ieee754_rem_pio2. These functions
are internal to the libm and exported through Symbol.map; and thus,
reserved for the implementation.
PR: 272783
MFC after: 1 week
patch to fix half-cycle trigonometric functions
Paul Zimmermann, a MPFR developer, contacted me about large errors in
the half-cycle trigonometric functions. I've have investigated these
issues and developed the attached patch. The float, double, and ld80
(long double) changes have been tested.
Caveat emptor: The ld128 changes have not been compiled. The ld128
changes have not been tested. I do not have access to a system that
uses ld128 floating point.
Here is an itemized list of changes:
* lib/msun/src/math_private.h:
. Add fast floor macros to compute the integer part of |x| for
0 <= |x| 01xp(N-1), where N is the precision of the type of x.
These macros are used in the half-cycle trigonometric functions
(e.g., sinpi(x)).
. The FFLOOR80 macros is used with the Intel 80-bit extended double
functions. This macors corrects an off-by-one error, which led to
enormous error for |x| > 0x1p32.
* lib/msun/src/s_cospif.c:
* lib/msun/src/s_cospi.c:
* lib/msun/ld80/s_cospil.c:
. Update Copyright years.
. Use FFLOOR*() macro to get integer part of |x|.
. Correct handle the range 0x1p(N-1) <= |x| < 0x1pN. Here, one needs
to determine if the integral value of |x| is even or odd to choose
+1 or -1. If |x| >= 0x1pN, always return +1.
* lib/msun/src/s_sinpif.c:
* lib/msun/src/s_sinpi.c:
* lib/msun/ld80/s_sinpil.c:
. Update Copyright years.
. Use FFLOOR*() macro to get integer part of |x|.
* lib/msun/src/s_tanpif.c:
* lib/msun/src/s_tanpi.c:
* lib/msun/ld80/s_tanpil.c:
. Update Copyright years.
. For +-0.5, return +-inf. Previously, tanpi[fl]() returned an NaN.
. Use FFLOOR*() to get integer part of |x|. Need to determine if the
integer part is even or odd. This is used to set +-0 for |x|
integral
and +-inf for (n+1/2).
. For 0x1p(N-1) <= |x| < 0x1pN need to determine if x is an even or
odd
integer to select +0 or -0. For |x| >= 0x1pN, it is always an even
integer, select 0.
. Note, tanpi[fl](x) is an odd function, so one needs to consider
tanpi[fl](-|x|) = - tanpi[fl](|x|).
* lib/msun/ld128/s_cospil.c:
* lib/msun/ld128/s_sinpil.c:
* lib/msun/ld128/s_tanpil.c:
. Update Copyright years.
. These routines use an FFLOOR128 macros, which likely should be
replaced by a bit twiddling algorithm.
. The same considerations above are applied to 0x1p112 <= |x| <
0x1p113,
and |x| >= 0x1p113 cases.
. Note, even and odd determination used fmodl(x,2.), which is likely
slow.
PR: 272742
MFC after: 1 week
David Das (das@) committed Bruce Evan's (bde's) WIP code for
expl() and logl() in git revision 25a4d6bfda. That code
included instrumentation that allowed bde to generate pari
scripts used in testing/debugging. This patch removes that
instrumentation as it is unlikely that others will ever use it.
* math/libm/msun/src/math_private.h:
. Remove bde's macros for the generation of pari scripts.
* math/libm/msun/ld128/s_expl.c:
* math/libm/msun/ld128/s_logl.c:
* math/libm/msun/ld80/s_expl.c:
* math/libm/msun/ld80/s_logl.c:
. Remove bde's DOPRINT_START macro.
. Change RETURNP to RETURNF.
. Change RETURN2P to RETURNF. Adjust arguments as needed.
. Change RETURNPI to RETURNI.
. Change RETURN2PI to RETURNI. Adjust arguments as needed.
PR: 272765
MFC after: 1 week
In order to compile lib32 libraries and other 32-bit code on arm64,
<machine/foo.h> needs to be redirected to an arm header rather
than arm64 when building with -m32. Ifdef the arm64 headers that
are installed in /usr/include/machine and used by user-level software
(including references from /usr/include/*.h) so that if __arm__ is
defined when including the arm64 version, <arm/foo.h> is included
rather than using the rest of the file's contents. Some arm headers
had no arm64 equivalent; headers were added just to do the redirection.
These files use #error if __arm__ is not defined to guard against
confusion. Also add an include/arm Makefile, and modify Makefiles
as needed to install everything, including the arm files in
/usr/include/arm. fenv.h comes from lib/msun/arm/fenv.h.
The new arm64 headers are:
acle-compat.h
cpuinfo.h
sysreg.h
Reviewed by: jrtc27, imp
Differential Revision: https://reviews.freebsd.org/D40944
The current versions of lib/msun/src/s_cospi.c, s_sinpi.c and s_tanpi.c
all exhibit the same defect. After checking for various numeric ranges,
they check to see whether the input argument is a NaN or an Infinity.
However, the code uses a value of 0x7f80000 instead of the correct value
of 0x7ff00000.
If you review s_cospif.c, s_sinpif.c, and s_tanpif.c, you will see that
the equivalent statements in these functions are accurate and have
appropriate source comments.
The impact of these defects is to flag some valid input values as
invalid and raise a pole error (divide by zero).
Reported by: Paul Green <Paul.Green@stratus.com>
PR: 272539
MFC after: 1 week
The sincos() man page notes the function was added to msun in FreeBSD
9.0 which must have been an oversight in the review as it was commited
to 12.0 and then backported to the 11 branch.
So I have provided a diff to correct this to the first FreeBSD version
it did ship with which was 11.2.
Reviewed by: dim, imp
MFC after: 3 days
Differential Revision: https://reviews.freebsd.org/D40308
The SPDX folks have obsoleted the BSD-2-Clause-FreeBSD identifier. Catch
up to that fact and revert to their recommended match of BSD-2-Clause.
Discussed with: pfg
MFC After: 3 days
Sponsored by: Netflix
Simplify the tests for 32-bit arm soft float support. For the files
included only on arm, drop the test entirely. For others, test
MACHINE_CPUARCH against arm.
No functional change intended. File lists appear the same before / after
the change.
Sponsored by: Netflix
Reviewed by: emaste
Differential Revision: https://reviews.freebsd.org/D38582
Summary:
This knob can be used to make buildsystem prefer generic C implentations of
various functions, instead of machine-specific assembler ones.
Test Plan: `make buildworld` on amd64
Reviewed by: imp, emaste
Differential Revision: https://reviews.freebsd.org/D36076
MFC after: 3 days
- s/modfied/modified/
- s/minimun/minimum/
While here, fix some mandoc warnings:
- whitespace at end of input line
- unusual Xr punctuation
- missing comma before name
Obtained from: NetBSD
MFC after: 5 days
Damian McGuckin <damianm at esi dot com dot au> noted that the accuracy
claims in the code for cbrt(3) and cbrtl(3) were incorrect. Fix the
comments to more accurately describe the accuracies.
PR: 265603
MFC after: 3 days
Since https://github.com/llvm/llvm-project/commit/ca75ac5f04f2, clang 15
has a new warning about _Generic selection expressions, such as used in
math.h:
lib/libc/gdtoa/_ldtoa.c:82:10: error: due to lvalue conversion of the controlling expression, association of type 'volatile float' will never be selected because it is qualified [-Werror,-Wunreachable-code-generic-assoc]
switch (fpclassify(u.e)) {
^
lib/msun/src/math.h:109:2: note: expanded from macro 'fpclassify'
__fp_type_select(x, __fpclassifyf, __fpclassifyd, __fpclassifyl)
^
lib/msun/src/math.h:85:14: note: expanded from macro '__fp_type_select'
volatile float: f(x), \
^
This is because the controlling expression always undergoes lvalue
conversion first, dropping any cv-qualifiers. The 'const', 'volatile',
and 'volatile const' associations will therefore never be used.
MFC after: 1 week
Reviewed by: theraven
Differential Revision: https://reviews.freebsd.org/D35815
Summary:
These functions are missing from the library itself, and exist solely in
the header. This breaks a few ports that expect libm to have the
symbols in the library itself.
Questions on MFC-ability: Can this be MFC'd to 13.2, and how?
Reviewers: imp, emaste, kib
Reviewed By: kib
Differential Revision: https://reviews.freebsd.org/D35204
There are some sections which could be improved
and work to do so is on going. The work will be
covered via 'X-MFC-WITH' commits.
Obtained from: OpenBSD
MFC after: 1 month
Differential Revision: https://reviews.freebsd.org/D34759
All supported compilers (modern versions of GCC and clang) support
this.
Many places didn't have an #else so would just silently do the wrong
thing. Ancient versions of icc (the original motivation for this) are
no longer a compiler FreeBSD supports.
PR: 263102 (exp-run)
Reviewed by: brooks, imp
Differential Revision: https://reviews.freebsd.org/D34797
. Disconnect imprecise.c from the build. This file can be deleted.
. Add b_tgammal.c to the build for ld80 and ld128 targets. The ld128
is a 'git mv' of imprecise.c to ld128/b_tgammal.c.
* lib/msun/ld80/b_expl.c:
. New file. Implement __exp__D for ld80 targets. This is based on
bsdsrc/b_exp.c.
* lib/msun/ld80/b_logl.c:
. New file. Implement __log__D for ld80 targets. This is based on
bsdsrc/b_log.c.
* lib/msun/ld80/b_tgammal.c b/lib/msun/ld80/b_tgammal.c
. New file. Implement tgammal(x) for ld80 targets.
Submitted by: Steve Kargl
Differential Revision: https://reviews.freebsd.org/D33444
Reviewed by: pfg
. Disconnect b_exp.c and b_log.c from the build.
* lib/msun/bsdsrc/b_exp.c:
. Replace scalb() usage with C99's ldexp().
. Replace finite(x) usage with C99's isfinite().
. Whitespace changes towards style(9).
. Remove include of "mathimpl.h". It is no longer needed.
. Remove #if 0 ... #endif code, which has been present since svn r93211
(2002-03-26).
. New minimax polynomial coefficients.
. Add comments to explain origins of some constants.
. Use ansi-C prototype. Remove K&R prototype. Add static to prototype.
* lib/msun/bsdsrc/b_log.c:
. Remove include of "mathimpl.h". It is no longer needed.
. Fix comments to actually describe the code.
. Reduce minimax polynomial from degree 4 to degree 3.
This uses newly computed coefficients.
. Use ansi-C prototype. Remove K&R prototype. Add static to prototype.
. Remove volatile in declaration of u1.
. Alphabetize decalaration list.
. Whitespace changes towards style(9).
. In argument reduction of x to g and m, replace use of logb() and
ldexp() with a single call to frexp(). Add code to get 1 <= g < 2.
. Remove #if 0 ... #endif code, which has been present since svn r93211
(2002-03-26).
. The special case m == -1022, replace logb() with ilogb().
* lib/msun/bsdsrc/b_tgamma.c:
. Update comments. Fix comments where needed.
. Add float.h to get LDBL_MANT_DIG for weak reference of tgammal to tgamma.
. Remove include of "mathimpl.h". It is no longer needed.
. Use "math.h" instead of <math.h>.
. Add '#include math_private.h"
. Add struct Double from mathimpl.h and include b_log.c and b_exp.c.
. Remove forward declarations of neg_gam(), small_gam(), smaller_gam,
large_gam() and ratfun_gam() by re-arranging the code to move these
function above their first reference.
. New minimax coefficients for polynomial in large_gam().
. New splitting of a0 into a0hi nd a0lo, which include additional
bits of precision.
. Use ansi-C prototype. Remove K&R prototype.
. Replace the TRUNC() macro with a simple cast of a double entities
to float before assignment (functional changes).
. Replace sin(M_PI*z) with sinpi(z) and cos(M_PI*(0.5-z)) with cospi(0.5-z).
Submitted by: Steve Kargl
Differential Revision: https://reviews.freebsd.org/D33444
Reviewed by: pfg
This look like a copy and paste leftover.
Reported by: enh@google.com (via freebsd-numerics@)
Reviewed by: Steve Kargl <sgk@troutmask.apl.washington.edu>
MFC after: 3 days
Some of the functions in msun can be implemented using a compiler
builtin function to generate a small number of instructions. Implement
this support in fma, fmax, fmin, and sqrt on arm64.
Care must be taken as the builtin can be implemented as a function
call on some architectures that lack direct support. In these cases
we need to use the original code path.
As we don't set errno on failure build with -fno-math-errno so the
toolchain doesn't convert a builtin into a function call when it
detects a failure, e.g. gcc will add a call to sqrt when the input
is negative leading to an infinite loop.
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D32801
In ac76bc1145, I added a few volatiles to work around ctrig_test
failures with {inf,inf}. This is not necessary anymore now, since in
3b00222f15 we added -fp-exception-behavior=maytrap for clang >= 10 in
libm's Makefile. (The flag tells clang to use stricter floating point
semantics, which libm depends on.)
PR: 244732, 254911
Fixes: ac76bc1145
MFC after: 3 days
The change implements cexpl() for both ld80 and ld128 architectures.
Testing was done on x86_64 and aarch64 systems.
Along the way sincos[fl]() use an optimization that reduces the argument
to being done one rather than twice. This optimization actually pointed
to a bug in the ld128 version of sincosl(), which is now fixed. In
addition, the minmax polynomial coefficients for sincosl() have been
updated.
A concise log of the file-by-file changes follows.
* include/complex.h:
. Add a prototype for cexpl().
* lib/msun/Makefile:
. Add s_cexpl.c to the build.
. Setup a link for cexpl.3 to cexp.3.
* lib/msun/Symbol.map:
. Expose cexpl symbol in libm shared library.
* lib/msun/ld128/s_cexpl.c:
* Implementation of cexpl() for 128-bit long double architectures.
Tested on an aarch64 system.
* lib/msun/ld80/s_cexpl.c:
* Implementation of cexpl() for Intel 80-bit long double.
* lib/msun/man/cexp.3:
. Document cexpl().
* lib/msun/man/complex.3:
. Add a BUGS section about cpow[fl].
* lib/msun/src/s_cexp.c:
. Include float.h for weak references on 53-bit long double targets.
. Use sincos() to reduce argument reduction cost.
* lib/msun/src/s_cexpf.c:
. Use sincosf() to reduce argument reduction cost.
* lib/msun/src/k_sincosl.h:
. Catch up with the new minmax polynomial coefficients for the kernel for
the 128-bit cosl() implementation.
. BUG FIX: *cs was used where *sn should have been. This means that sinl()
was no computed correctly when iy != 0.
* lib/msun/src/s_cosl.c:
. Include fpmath.h to get access to IEEEl2bits.
. Replace M_PI_4 with pio4, a 64-bit or 113-bit approximation for pi / 4.
PR: 216862
MFC after: 1 week
As mention previously, the minmax polynomial approximation
in the kernel for cosl() seem to have a bad set of coefficients.
In testing, cosl() in the interval [0.785, pi/4] for 1 million
values and pi/4 written to 37 decimal digits. The old version
on an aarch64 system gave
% tlibm/tlibm_lmath -l -x 0.78 -X
7.85398163397448309615660845819875721e-1L cos
Interval tested for cosl: [0.78,0.785398]
count: 1000000
xm = 7.80213913234863919029058821396125599e-01L
libm = 7.10763080972549562455058499280609083e-01L
mpfr = 7.10763080972549562455058499280608983e-01L
ULP = 1.04431
The max ULP exceeds 1, which is not good. So, I rinsed off a 10
year code and recomputed coefficients. The new minmax polynomial
now yields
% tlibm/tlibm_lmath -l -x 0.78 -X
7.85398163397448309615660845819875721e-1L cos
Interval tested for cosl: [0.78,0.785398]
count: 1000000
xm = 7.82916198746768272588844890973704219e-01L
libm = 7.08859615479571058183956453286628396e-01L
mpfr = 7.08859615479571058183956453286628469e-01L
ULP = 0.75407
which is very good.
PR: 218514
MFC after: 1 week
Mark Murray graciously provided access to an aarch64 system
to test the ld128 implementations. This patch address
* Misuses of copysignl() in sinpil() and tanpil().
* Redo the splitting of argument 'x' into an integer part and
remainder. The remainder must satify 0 <= r < 1.
* Update the reduction of the integer part to something that can
easily be seen as even or odd, e.g., sin(pi*x) = (-1)^n*sin(pi*r)
with n <= 2^112 and we an reduce n by subtracting integer powers
of 2.
* In s_cospil.c, fix typos where 'x' is used where 'ax', the
remainder, is required.
* In tanpil(), fix the use of an uninitialized variable, ax = fabsl(ax),
ax should be x in fabsl().
One item of note, in the limited tested on aarch64, the max ULP
for sinpil() and cospil() were less than 1.1 ULP, which is higher
that the desired max ULP less than 1. This was traced to the
kernel for cosl() in the fundamental interval [0,pi/4].
The coefficients in the minmax polynomial likely need refinement.
PR: 218514
MFC after: 1 week
The patch fixes the omission of '#include <float.h>', which is needed for
the weak reference on systems with LDBL_MANT_DIG == DBL_MANT_DIG.
PR: 218514
MFC after: 1 week
Both IEEE-754 2008 and ISO/IEC TS 18661-4 define the half-cycle
trignometric functions cospi, sinpi, and tanpi. The attached
patch implements cospi[fl], sinpi[fl], and tanpi[fl]. Limited
testing on the cospi and sinpi reveal a max ULP less than 0.89;
while tanpi is more problematic with a max ULP less than 2.01
in the interval [0,0.5]. The algorithms used in these functions
are documented in {ks}_cospi.c, {ks}_sinpi.c, and s_tanpi.c.
Note. I no longer have access to a system with ld128 and
adequate support to compile and test the ld128 implementations
of these functions. Given the almost complete lack of input from
others on improvements to libm, I doubt that anyone cares. If
someone does care, the ld128 files contain a number of FIXME comments,
and in particular, while the polynomial coefficients are given
I did not update the polynomial algorithms to properly use the
coefficients.
PR: 218514
MFC after: 2 weeks
These files were copied from MUSL. Add the standard copyright notice and
SPDX-License-Identifier: MIT consistent with our new draft license
policy. It reads word for word the same as the MIT license on the SPDX
web site. Add a pointer to the MUSL COPYIRGHT file which contains a list
of all authors of MUSL.
Sponsored by: Netflix
Noticed by: Steve Kargl
Summary:
From Steve Kargl:
Paul Zimmermann has identified a bug in Openlibm's powf(),
which is identical to FreeBSD's libm. Both derived from
fdlibm. https://github.com/JuliaMath/openlibm/issues/212.
Consider
% cat h.c
int
main(void)
{
float x, y, z;
x = 0x1.ffffecp-1F;
y = -0x1.000002p+27F;
z = 0x1.557a86p115F;
printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z);
return 0;
}
% cc -o h -fno-builtin h.c -lm && ./h
9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34
Reviewers: manu
Subscribers: imp, andrew, emaste
Differential Revision: https://reviews.freebsd.org/D31865
.Fa is the suitable macro for functions in comparsion to the
.Ar macro, which should be used for commandline arguments.
While here, fix some mandoc warnings.
Reviewed by: imp (earlier version)
Obtained from: OpenBSD (in partial)
MFC after: 3 days
Differential Revision: https://reviews.freebsd.org/D31090
After df3b437c1e, older gcc's such as
4.2.1 (still used on earlier branches for e.g. mips and powerpc) and
6.3.0 (still used for some cross-builds) started throwing bogus errors
like:
In file included from /workspace/src/lib/msun/src/s_llround.c:11:0:
/workspace/src/lib/msun/src/s_lround.c:54:31: error: initializer element is not constant
static const type dtype_min = type_min - 0.5;
^~~~~~~~
/workspace/src/lib/msun/src/s_lround.c:55:31: error: initializer element is not constant
static const type dtype_max = type_max + 0.5;
^~~~~~~~
Since 'type_min' and 'type_max' are constants declared just above these
lines this error is nonsensical, but older gcc's are not smart enough.
Work around the error by reusing the (type)DTYPE_MIN and (type)DTYPE_MAX
macros, so I can MFC this right away, unbreaking a few stable builds.
MFC after: immediately
It turned out that the (type)DTYPE_MAX conversions at the top of
s_lround.c are now emitted as cvtsi2sd instructions, at least on SSE
capable CPUs. This caused the FE_INEXACT flag to always be set, at least
for the double and float variants. Under clang 11, the whole INRANGE()
comparisons were still optimized away, but this has "improved" in clang
12, due to stricter adherence to the -ffp-exception-behavior=maytrap
compiler flag.
To avoid run-time integer to float conversions, use static constants
instead, so they are computed at compile time, and the INRANGE()
statements are optimized away again, if applicable.
While here, use an integer instead of a floating type to store the test
results in lround_test.c, as this is more appropriate, and we can also
drop the volatile hack.
Reported by: arichardson
MFC after: 3 days
Summary:
Fix FPU exception management for powerpcspe. Bits are in a different place from
the standard FPSCR, so we need to handle the shifting differences. Also,
there's no concept of a "software exception" raise, so we need to do exceptional
math to trigger the exception from software.
Reviewed By: alfredo
Differential Revision: https://reviews.freebsd.org/D22824
For some reason the ld128 log1pl() implementation is less accurate than
logl(), but does at least guarantee precision >= the ld80 implementation.
Mark log1p_accuracy_tests as XFAIL for ld128 and increase the log1p tolerance
to the ld80 equivalent in accuracy_tests to avoid losing test coverage for
the other functions.
PR: 253984
Reviewed By: ngie, dim
Differential Revision: https://reviews.freebsd.org/D29039
