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 gi

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

show more ...

lib/msun: Cleanup after $FreeBSD$ removal

Remove no longer needed explicit inclusion of sys/cdefs.h.

PR: 276669
MFC after: 1 week

Remove $FreeBSD$: one-line .c pattern

Remove /^[\s*]*__FBSDID\("\$FreeBSD\$"\);?\s*\n/

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 elementar

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 3f70824172feb82ea3dcdb3866b54fe0eb7cd890
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

show more ...

Minor improvements:
- Improve the order of some tests.
- Fix style.

Submitted by: bde

Fix some problems with asinf(), acosf(), atanf(), and atan2f():

- Adjust several constants for float precision. Some thresholds
that were appropriate for double precision were never changed
when

Fix some problems with asinf(), acosf(), atanf(), and atan2f():

- Adjust several constants for float precision. Some thresholds
that were appropriate for double precision were never changed
when these routines were converted to float precision. This
has an impact on performance but not accuracy. (Submitted by bde.)

- Reduce the degrees of the polynomials used. A smaller degree
suffices for float precision.

- In asinf(), use double arithmetic in part of the calculation to
avoid a corner case and some complicated arithmetic involving a
division and some buggy constants. This improves performance and
accuracy.

Max error (ulps):
asinf acosf atanf
before 0.925 0.782 0.852
after 0.743 0.804 0.852

As bde points out, it's cheaper for asin*() and acos*() to use
polynomials instead of rational functions, but that's a task for
another day.

show more ...

As in other parts of libm, mark a few constants as volatile to prevent
spurious optimizations. gcc doesn't support FENV_ACCESS, so when it
folds constants, it assumes that the rounding mode is always

As in other parts of libm, mark a few constants as volatile to prevent
spurious optimizations. gcc doesn't support FENV_ACCESS, so when it
folds constants, it assumes that the rounding mode is always the
default and floating point exceptions never matter.

show more ...

s/rcsid/__FBSDID/

Assume __STDC__, remove non-__STDC__ code.

Submitted by: keramida

$Id$ -> $FreeBSD$

Use __ieee754_sqrt() instead of sqrt() internally. Similarly for the
float versions. Using sqrt() was inefficient.

Obtained from: NetBSD

Revert $FreeBSD$ to $Id$

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so long.

Boy, I'm glad we're not using sup anymore. This update would have been
insane otherwise.

show more ...

Remove trailing whitespace.

J.T. Conklin's latest version of the Sun math library.

-- Begin comments from J.T. Conklin:
The most significant improvement is the addition of "float" versions
of the math functions that take float

J.T. Conklin's latest version of the Sun math library.

-- Begin comments from J.T. Conklin:
The most significant improvement is the addition of "float" versions
of the math functions that take float arguments, return floats, and do
all operations in floating point. This doesn't help (performance)
much on the i386, but they are still nice to have.

The float versions were orginally done by Cygnus' Ian Taylor when
fdlibm was integrated into the libm we support for embedded systems.
I gave Ian a copy of my libm as a starting point since I had already
fixed a lot of bugs & problems in Sun's original code. After he was
done, I cleaned it up a bit and integrated the changes back into my
libm.
-- End comments

Reviewed by: jkh
Submitted by: jtc

show more ...

Minor improvements:
- Improve the order of some tests.
- Fix style.

Submitted by: bde

Fix some problems with asinf(), acosf(), atanf(), and atan2f():

- Adjust several constants for float precision. Some thresholds
that were appropriate for double precision were never changed
when

Fix some problems with asinf(), acosf(), atanf(), and atan2f():

- Adjust several constants for float precision. Some thresholds
that were appropriate for double precision were never changed
when these routines were converted to float precision. This
has an impact on performance but not accuracy. (Submitted by bde.)

- Reduce the degrees of the polynomials used. A smaller degree
suffices for float precision.

- In asinf(), use double arithmetic in part of the calculation to
avoid a corner case and some complicated arithmetic involving a
division and some buggy constants. This improves performance and
accuracy.

Max error (ulps):
asinf acosf atanf
before 0.925 0.782 0.852
after 0.743 0.804 0.852

As bde points out, it's cheaper for asin*() and acos*() to use
polynomials instead of rational functions, but that's a task for
another day.

show more ...

As in other parts of libm, mark a few constants as volatile to prevent
spurious optimizations. gcc doesn't support FENV_ACCESS, so when it
folds constants, it assumes that the rounding mode is always

As in other parts of libm, mark a few constants as volatile to prevent
spurious optimizations. gcc doesn't support FENV_ACCESS, so when it
folds constants, it assumes that the rounding mode is always the
default and floating point exceptions never matter.

show more ...

s/rcsid/__FBSDID/

Assume __STDC__, remove non-__STDC__ code.

Submitted by: keramida

$Id$ -> $FreeBSD$

Use __ieee754_sqrt() instead of sqrt() internally. Similarly for the
float versions. Using sqrt() was inefficient.

Obtained from: NetBSD

Revert $FreeBSD$ to $Id$

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so

Make the long-awaited change from $Id$ to $FreeBSD$

This will make a number of things easier in the future, as well as (finally!)
avoiding the Id-smashing problem which has plagued developers for so long.

Boy, I'm glad we're not using sup anymore. This update would have been
insane otherwise.

show more ...

Remove trailing whitespace.