Commit graph

540 commits

Author SHA1 Message Date
Steve Kargl
41e016289f Improve accuracy of asinf(3) and acosf(3)
This uses a better rational approximation to improve the accuracy of
both functions. For exhaustive testing of asinf(3) in the interval, the
current libm gives:

    % ./tlibm asin -fPED -x 0x1p-12f -X 1
    Interval tested for asinf: [0.000244141,1]
           ulp <= 0.5:  97.916% 98564994 |  97.916% 98564994
    0.5 <  ulp <  0.6:  2.038% 2051023 |  99.953% 100616017
    0.6 <  ulp <  0.7:  0.047%   47254 | 100.000% 100663271
    0.7 <  ulp <  0.8:  0.000%      25 | 100.000% 100663296
    Max ulp: 0.729891 at 5.00732839e-01

which isn't too bad given that much of the computation is actually done
in double floating point.

With the new rational approximation, exhaustive testing yields:

    % ./tlibm asin -fPED -x 0x1p-12f -X 1
    Interval tested for asinf: [0.000244141,1]
           ulp <= 0.5:  99.711% 100372643 |  99.711% 100372643
    0.5 <  ulp <  0.6:  0.288%  290357 | 100.000% 100663000
    0.6 <  ulp <  0.7:  0.000%     296 | 100.000% 100663296
    Max ulp: 0.636344 at 5.09706438e-01

Similarly, for exhaustive testing of asinf(3) in the interval, the
current libm gives:

    % ./tlibm acos -fPED -x -1 -X -0x1p-12f
    Interval tested for acosf: [-1,-0.000244141]
           ulp <= 0.5:  97.008% 97651921 |  97.008% 97651921
    0.5 <  ulp <  0.6:   2.441%  2457242 |  99.450% 100109163
    0.6 <  ulp <  0.7:   0.472%   475503 |  99.922% 100584666
    0.7 <  ulp <  0.8:   0.071%    71309 |  99.993% 100655975
    0.8 <  ulp <  0.9:   0.007%     7319 | 100.000% 100663294
    0.9 <  ulp <  1.0:   0.000%        2 | 100.000% 100663296
    Max ulp: 0.914007 at -5.01484931e-01

    % ./tlibm acos -fPED -x 0x1p-12f -X 1
    Interval tested for acosf: [0.000244141,1]
           ulp <= 0.5:  97.317% 97962530 |  97.317% 97962530
    0.5 <  ulp <  0.6:   2.340%  2355182 |  99.657% 100317712
    0.6 <  ulp <  0.7:   0.314%   316134 |  99.971% 100633846
    0.7 <  ulp <  0.8:   0.029%    29450 | 100.000% 100663296
    Max ulp: 0.796035 at 4.99814630e-01

With the new rational approximation, exhaustive testing yields:

    % ./tlibm acos -fPED -x -1 -X -0x1p-12f
    Interval tested for acosf: [-1,-0.000244141]
           ulp <= 0.5:  97.010% 97653245 |  97.010% 97653245
    0.5 <  ulp <  0.6:   2.442%  2458373 |  99.452% 100111618
    0.6 <  ulp <  0.7:   0.473%   476012 |  99.925% 100587630
    0.7 <  ulp <  0.8:   0.068%    68603 |  99.993% 100656233
    0.8 <  ulp <  0.9:   0.007%     7063 | 100.000% 100663296
    Max ulp: 0.896189 at -5.04511118e-01

    % ./tlibm acos -fPED -x 0x1p-12f -X 1
    Interval tested for acosf: [0.000244141,1]
           ulp <= 0.5:  97.650% 98298175 |  97.650% 98298175
    0.5 <  ulp <  0.6:   2.028%  2041709 |  99.679% 100339884
    0.6 <  ulp <  0.7:   0.292%   293555 |  99.970% 100633439
    0.7 <  ulp <  0.8:   0.030%    29857 | 100.000% 100663296
    Max ulp: 0.775875 at 4.91849005e-01

PR:		281001
MFC after:	1 week
2024-08-29 21:44:48 +02:00
Steve Kargl
213eb102ae
msun: Fix typo in comment
PR:		280965
MFC after:	3 days
2024-08-21 14:59:07 +08:00
Steve Kargl
34f746cc7f libm: fma: correct zero sign with small inputs
This is a fixed version of 888796ade2.

PR:		277783
Reported by:	Victor Stinner
Reviewed by:	emaste
MFC after:	1 week
2024-07-28 17:37:45 -04:00
Ed Maste
001606523a libm: add parens to clarify expressions in fma, fmal
Obtained from:	NetBSD
2024-07-15 09:30:02 -04:00
Ryan Libby
07cc7ea738 libmsun: remove duplicates after cdefs.h added inline to __always_inline
Reviewed by:	kib, olce
Sponsored by:	Dell EMC Isilon
Differential Revision:	https://reviews.freebsd.org/D45712
2024-06-25 10:40:14 -07:00
Warner Losh
fa1e8538d1 math.h: Remove support for old gcc versions
Remove support for old gcc versions by unconditionally defining
__MATH_BUILTIN_RELOPS and __MATH_BUILTIN_CONSTANTS. Per kib's request,
don't #undef those so it's easier to understand what the builtins are
doing.

Reviewed by:		brooks
Differential Revision:	https://reviews.freebsd.org/D45655
Sponsored by:		Netflix
2024-06-20 20:41:09 -06:00
Ed Maste
e77ad954bb Revert "libm: fma: correct zero sign with small inputs"
This change introduced a test failure, so revert until that can be
addressed.

This reverts commit 888796ade2.

PR:		277783
Reported by:	rlibby
Sponsored by:	The FreeBSD Foundation
2024-06-11 21:36:12 -04:00
Ed Maste
92927b8bcf msun: update Clang bug reference in fma test
LLVM bugzilla bug 8100 became issue #8472 with the migration to GitHub.

https://github.com/llvm/llvm-project/issues/8472
2024-06-11 20:29:27 -04:00
Ed Maste
888796ade2 libm: fma: correct zero sign with small inputs
PR:		277783
Reported by:	Victor Stinner
Submitted by:	kargl
MFC after:	1 week
Differential Revision: https://reviews.freebsd.org/D44433
2024-06-08 11:55:36 -04:00
Karl Tomlinson
e75a1bbc23 pow,powf(3),__ieee754_rem_pio2(f): Avoid negative integer left shift UB
A compiler clever enough to know that z is positive with a non-zero
biased exponent could, for example, optimize away the scalbnf(z,n) in
pow() because behavior for left shift of negative values is undefined.
`n` is negative when y*log2(|x|) < -0.5.  i.e. |x^y| < sqrt(0.5)

The intended behavior for operator<< in this code is to shift the two's
complement representation of the first operand.

In the pow() functions, the result is added to the IEEE 754 exponent of
z = 2^y'.  n may be negative enough to underflow the biased IEEE 754
exponent below zero, which is manifested in the sign bit of j
(which would correspond to the IEEE 754 sign bit).

The conversion from uint32_t to int32_t for out-of-int32_t-range values
is implementation defined.  The assumed behavior of interpreting the
uint32_t value as a two's complement representation of a signed value
is already assumed in many parts of the code, such as uses of
GET_FLOAT_WORD() with signed integers.

This code passes all the current tests, and makes some out of tree
fuzzing tests pass again rather than hit UB (detailed in the commentary
of the pull request).

Signed-off-by: Karl Tomlinson <karlt+@karlt.net>
Reviewed by: imp, steve kargl, dim
Pull Request: https://github.com/freebsd/freebsd-src/pull/1137
2024-04-23 14:04:07 -06:00
Henri Chataing
2b2cd97844 msun: Fix math error in comment explaining y reduction
x = k + y for some integer k and |y| < 1/2
exp2(x) = exp2(k + y) = exp2(k) * exp2(y)
which can be written as 2**k * exp2(y)

The original had x = 2**k + y, which is has an extra 2** in it which
isn't correct.

Confirmed by forumula 2 in Gal and Bachelis referenced in the comments
for the source of this method
	https://dl.acm.org/doi/pdf/10.1145/103147.103151

The actual code is correct.

Reviewed by: imp (who added s_exp2.c and wrote the commit message)
Pull Request: https://github.com/freebsd/freebsd-src/pull/1127
2024-04-12 16:15:04 -06:00
Collin Funk
22ddb5bb83 math: Add long double constant definitions
These constants are GNU libc extensions that are likely to be adopted
by the next POSIX revision [1]. The definitions can be verified in
a PARI-GP shell session:

* M_El: exp (1)
* M_LOG2El: log (exp (1)) / log (2)
* M_LOG10El: log (exp (1)) / log (10)
* M_LN2l: log (2)
* M_LN10l: log (10)
* M_PIl: Pi
* M_PI_2l: Pi / 2
* M_PI_4l: Pi / 4
* M_1_PIl: 1 / Pi
* M_2_PIl: 2 / Pi
* M_2_SQRTPIl: 2 / sqrt (Pi)
* M_SQRT2l: sqrt (2)
* M_SQRT1_2l: 1 / sqrt (2)

[1] https://austingroupbugs.net/view.php?id=828

Put these behind __BSD_VISIBLE || __XSI_VISIBLE >= 800 to future-proof
these changes. They shouldn't be defined at lower levels of XSI, but don't
have other XSI 800 stuff in place yet.

Signed-off-by: Collin Funk <collin.funk1@gmail.com>
Reviewed by: imp, allanjude
Pull Request: https://github.com/freebsd/freebsd-src/pull/1121
2024-04-12 12:30:58 -06:00
Martin Oliveira
4309f9f234 include/math.h: fix warning with -Wconversion
The way the __fp_type_select macro uses the _Generic expression causes
gcc to throw a warning on valid code if the -Wconversion flag is used.

For example, consider the following program:

    #include <math.h>
    int main()
    {
    	double x = 1.0;
    	isnan(x);
    	return 0;
    }

which throws a warning:

    $ gcc -Wconversion a.c
    a.c:5:15: warning: conversion from 'double' to 'float' may change value [-Wfloat-conversion]
        5 |         isnan(x);
          |               ^

This happens because the functions are invoked inside of the _Generic.
Looking at the example of _Generic in the C11 specification, one sees
that the parameters are outside of the _Generic expression (see page 79
here: https://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf).

Reference: https://stackoverflow.com/a/68309379
Signed-off-by: Martin Oliveira <martin.oliveira@eideticom.com>
Reviewed by: imp
Pull Request: https://github.com/freebsd/freebsd-src/pull/841
2024-02-02 21:05:49 -07:00
Steve Kargl
0dd5a5603e lib/msun: Cleanup after $FreeBSD$ removal
Remove no longer needed explicit inclusion of sys/cdefs.h.

PR:	276669
MFC after:	1 week
2024-01-28 17:00:23 +02:00
Warner Losh
dc36d6f9bb lib: Remove ancient SCCS tags.
Remove ancient SCCS tags from the tree, automated scripting, with two
minor fixup to keep things compiling. All the common forms in the tree
were removed with a perl script.

Sponsored by:		Netflix
2023-11-26 22:23:28 -07:00
Warner Losh
1b681154f5 math: Move to const instead of __const
There's no reason to use the __const construct here. This is a left-over
from supporting K&R and ANSI compilers in the original Sun msun. All
other K&R crutches have been removed. Remove these as well. There's no
semantic difference. And there's already several others in math.h.

Sponsored by:		Netflix
2023-11-19 20:55:57 -07:00
John Baldwin
eba230afba Purge more stray embedded $FreeBSD$ strings
These do not use __FBSDID but instead use bare char arrays.

Reviewed by:	imp, emaste
Differential Revision:	https://reviews.freebsd.org/D41957
2023-09-25 07:54:56 -07:00
Warner Losh
1d386b48a5 Remove $FreeBSD$: one-line .c pattern
Remove /^[\s*]*__FBSDID\("\$FreeBSD\$"\);?\s*\n/
2023-08-16 11:54:42 -06:00
Warner Losh
2a63c3be15 Remove $FreeBSD$: one-line .c comment pattern
Remove /^/[*/]\s*\$FreeBSD\$.*\n/
2023-08-16 11:54:29 -06:00
Warner Losh
42b388439b Remove $FreeBSD$: one-line .h pattern
Remove /^\s*\*+\s*\$FreeBSD\$.*$\n/
2023-08-16 11:54:23 -06:00
Warner Losh
b3e7694832 Remove $FreeBSD$: two-line .h pattern
Remove /^\s*\*\n \*\s+\$FreeBSD\$$\n/
2023-08-16 11:54:16 -06:00
Steve Kargl
99843eb899 Clean up libm use of the __ieee754_ prefix
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
2023-08-03 21:51:53 +02:00
Steve Kargl
2d3b0a687b Fixes for bugs in sinpi/cospi/tanpi
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
2023-08-03 07:27:58 +03:00
Steve Kargl
c66a499e03 Cleanup debugging code in libm
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
2023-08-03 07:27:58 +03:00
Steve Kargl
f2c94ddd0a
msun: Fix typo in math_private.h comment
PR:		272738
MFC after:	3 days
2023-07-27 02:21:44 +08:00
Steve Kargl
be4c7f2735 libm: correctly test for for NaN and Infinity in sinpi(), cospi(), and tanpi()
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
2023-07-17 08:23:27 +03:00
Warner Losh
4d846d260e spdx: The BSD-2-Clause-FreeBSD identifier is obsolete, drop -FreeBSD
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
2023-05-12 10:44:03 -06:00
Steve Kargl
620d855fac msun: correct comment
The comment in msun/src/e_jn.c lacks proper grammar, and is incorrect on
the choice of normalization entity.

PR:	266503
MFC after:	3 days
2022-09-19 21:40:07 +03:00
Gordon Bergling
82007616d0 msun: Remove a double word in a source code comment
- s/to to/to/

MFC after:	3 days
2022-09-10 12:59:10 +02:00
Gordon Bergling
a52f4499ae libm: Correct some typos in source code comments
- 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
2022-09-03 19:14:02 +02:00
Steve Kargl
369ea0520a [libm] Correct comments in s_cbrt[l].c
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
2022-08-04 19:33:34 +02:00
Dimitry Andric
e50027e38d Remove unnecessary const and volatile qualifiers from __fp_type_select()
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
2022-07-15 20:09:27 +02:00
Yi Kong
7e06f4708c
msun: Rewrite function definitions with identifier lists
This syntax is removed in C2x proposal N2432.

Reviewed by:	pfg
MFC after:	3 days
Differential Revision:	https://reviews.freebsd.org/D35771
2022-07-12 13:17:47 +08:00
John Baldwin
56f5947a71 Remove checks for __GNUCLIKE_ASM assuming it is always true.
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
2022-04-12 10:05:45 -07:00
Gordon Bergling
29fea59e78 math(3): Remove a double word in a source code comment
- s/is is/is/

MFC after:	3 days
2022-04-09 10:13:37 +02:00
Mark Murray
03a88e3de9 * lib/msun/Makefile b/lib/msun/Makefile:
. 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
2021-12-15 18:36:20 +00:00
Dimitry Andric
20d425842a Remove set-but-unused variable from s_sincosl.c
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
2021-12-14 22:50:30 +01:00
Andrew Turner
b2e843161d Use a builtin where possible in msun
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
2021-11-19 11:40:46 +00:00
Dimitry Andric
e2157cd000 Partially revert ac76bc1145 because it is no longer necessary
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
2021-11-05 22:27:20 +01:00
Steve Kargl
046e2d5db1 Implementations of cexpl()
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
2021-11-05 13:51:42 +02:00
Steve Kargl
3bfc837685 sinpi,cospi,tanpi: float.h needed for week reference
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
2021-10-29 03:15:19 +03:00
Steve Kargl
dce5f3abed [LIBM] implementations of sinpi[fl], cospi[fl], and tanpi[fl]
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
2021-10-26 02:50:20 +03:00
Warner Losh
3550a49f68 msun: Add copyright notices
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
2021-10-22 22:00:54 -06:00
Mark Murray
292815eac6 Fix powf().
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
2021-09-06 18:51:31 +01:00
Dimitry Andric
0bcd49c13a Work around bogus old gcc "initializer element is not constant" error
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
2021-06-25 20:43:20 +02:00
Dimitry Andric
df3b437c1e Fix failures in libm's lround_test after clang 12 import
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
2021-06-22 18:38:45 +02:00
Dimitry Andric
7702d940ec Avoid -pedantic warnings about using _Generic in __fp_type_select
When compiling parts of math.h with clang using a C standard before C11,
and using -pedantic, it will result in warnings similar to:

bug254714.c:5:11: warning: '_Generic' is a C11 extension [-Wc11-extensions]
  return !isfinite(1.0);
          ^
/usr/include/math.h:111:21: note: expanded from macro 'isfinite'
                    ^
/usr/include/math.h:82:39: note: expanded from macro '__fp_type_select'
                                      ^

This is because the block that enables use of _Generic is conditional
not only on C11, but also on whether the compiler advertises support for
C generic selections via __has_extension(c_generic_selections).

To work around the warning without having to pessimize the code, use the
__extension__ keyword, which is supported by both clang and gcc. While
here, remove the check for __clang__, as _Generic has been supported for
a long time by gcc too now.

Reported by:	yuri
PR:		254714
MFC after:	1 week
2021-04-08 18:20:32 +02:00
Alex Richardson
f5542795b9 s_scalbn.c: Add missing float.h include
This caused LDBL_MANT_DIG to not be defined and therefore the scalbnl
alias was not being emitted for double==long double platforms.

Fixes:		760b2ffc ("Update scalbn* functions to the musl versions")
Reported by:	Jenkins
2021-03-01 14:22:47 +00:00
Alex Richardson
760b2ffc55 Update scalbn* functions to the musl versions
The only diff compared to musl is a minor change to scalbnl() to replace
musl's union ldshape with union IEEEl2bits.
This fixes the scalbn tests on non-x86 (since x86 has an assembly version
that is used instead).

Musl commit messages:
commit 8c44a060243f04283ca68dad199aab90336141db
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Mon Apr 3 02:38:13 2017 +0200

    fix scalbn when result is in the subnormal range

    in nearest rounding mode scalbn could introduce double rounding error
    when an intermediate value and the final result were both in the
    subnormal range e.g.

      scalbn(0x1.7ffffffffffffp-1, -1073)

    returned 0x1p-1073 instead of 0x1p-1074, because the intermediate
    computation got rounded to 0x1.8p-1023.

    with the fix an intermediate value can only be in the subnormal range
    if the final result is 0 which is correct even after double rounding.
    (there still can be two roundings so signals may be raised twice, but
    that's only observable with trapping exceptions which is not supported.)

commit 2eaed464e2080d8321d3903b71086a1ecfc4ee4a
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Wed Sep 4 15:52:54 2013 +0000

    math: use float_t and double_t in scalbnf and scalbn

    remove STRICT_ASSIGN (c99 semantics is assumed) and use the conventional
    union to prepare the scaling factor (so libm.h is no longer needed)

commit 1b77b9072f374bd26eb0574b83a0d5f18d75ec60
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Thu Aug 15 10:07:46 2013 +0000

    math: minor scalbn*.c simplification

commit c4359e01303da2755fe7e8033826b132eb3659b1
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Tue Nov 13 10:55:35 2012 +0100

    math: excess precision fix modf, modff, scalbn, scalbnf

    old code was correct only if the result was stored (without the
    excess precision) or musl was compiled with -ffloat-store.
    now we use STRICT_ASSIGN to work around the issue.
    (see note 160 in c11 section 6.8.6.4)

commit 666271c105e4137bdfa195e217799d74143370d4
Author: Szabolcs Nagy <nsz@port70.net>
Date:   Tue Nov 13 10:30:40 2012 +0100

    math: fix scalbn and scalbnf on overflow/underflow

    old code was correct only if the result was stored (without the
    excess precision) or musl was compiled with -ffloat-store.
    (see note 160 in n1570.pdf section 6.8.6.4)

commit 8051e08e10d2b739fcfcbc6bc7466e8d77fa49f1
Author: nsz <nsz@port70.net>
Date:   Mon Mar 19 10:54:07 2012 +0100

    simplify scalbn*.c implementations

    The old scalbn.c was wrong and slow, the new one is just slow.
    (scalbn(0x1p+1023,-2097) should give 0x1p-1074, but the old code gave 0)

Reviewed By:	dim
Differential Revision: https://reviews.freebsd.org/D28872
2021-03-01 12:53:45 +00:00
Alex Richardson
a7b42c4b7f msun: ctanh/ctanhf: Import fix from musl libc
This applies musl commit b02eed9c4841913d690a2d0029737d72615384fe by
Szabolcs Nagy and updates the tests accordingly. This also allows
removing an XFAIL from the test.

musl commit message:

complex: fix ctanh(+-0+i*nan) and ctanh(+-0+-i*inf)

These cases were incorrect in C11 as described by
http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1886.htm

PR: 217528

Reviewed By:	dim
MFC after:	1 week
Differential Revision: https://reviews.freebsd.org/D28578
2021-02-15 22:55:12 +00:00