Synchronize most of libm with head as of r323004. This excludes a few
arch-specific updates for powerpcspe, mips and riscv, for which support
has not been merged yet.

Bump __FreeBSD_version for the addition of cacoshl, cacosl, casinhl,
casinl, catanl, catanhl, sincos, sincosf, and sincosl.

MFC r305382 (by bde):

Add asm versions of fmod(), fmodf() and fmodl() on amd64. Add asm
versions of fmodf() amd fmodl() on i387.

fmod is similar to remainder, and the C versions are 3 to 9 times
slower than the asm versions on x86 for both, but we had the strange
mixture of all 6 variants of remainder in asm and only 1 of 6
variants of fmod in asm.

MFC r305384 (by bde):

Disconnect the "optimized" asm variants of cos(), sin() and tan() from
the build on i386. Leave them in the source tree for regression tests.

The asm functions were always much less accurate (by a factor of more
than 10**18 in the worst case). They were faster on old CPUs. But
with each new generation of CPUs they get relatively slower. The
double precision C version's average advantage is about a factor of 2
on Haswell.

The asm functions were already intentionally avoided in float and long
double precision on i386 and in all precisions on amd64. Float
precision and amd64 give larger advantages to the C version. The long
double precision C code and compilers' understanding of long double
precision are not so good, so the i387 is still slightly faster for
long double precision, except for the unimportant subcase of huge args
where the sub-optimal C code now somehow beats the i387 by about a
factor of 2.

MFC r305385 (by bde):

Oops, the previous i386 version of e_fmodf.S and e_fmodl.S was
actually the amd64 version.

MFC r306409 (by emaste):

libm: fix some unused variable (rcsid) and dangling else warnings

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

MFC r306410 (by emaste):

libm: simplify i387 subdir logic with make's :S substitution

MFC r306527 (by emaste):

libm: remove unused variables for LDBL_MANT_DIG != 113

Sponsored by: The FreeBSD Foundation

MFC r306709 (by emaste):

libm: remove unused variables

Sponsored by: The FreeBSD Foundation

MFC r307066 (by br):

Don't use fmaxl/fminl on platforms with no long double support,
use fmax/fmin instead.

This fixes fmaxmin test failure on MIPS64.

Reviewed by: emaste
Sponsored by: DARPA, AFRL
Sponsored by: HEIF5
Differential Revision: https://reviews.freebsd.org/D8216

MFC r308172 (by emaste):

libm: add braces around initialization of subobjects

This cleans up a warning when building libm at higher WARNS levels and
makes the intent more clear. By the C standard the values are assigned
to subobject members in order so this change introduces no functional
change. (6.7.9 20)

Reviewed by: kib
Sponsored by: The FreeBSD Foundation
Differential Revision: https://reviews.freebsd.org/D8333

MFC r313761 (by mmokhi):

Add casinl() cacosl() catanl() casinhl() cacoshl() catanhl() APIs to msun
to improve C11 conformance.

PR: 216850 216851 216852 216856 216857 216858
Submitted by: mmokhi
Reported by: sgk@troutmask.apl.washington.edu
Reviewed by: bde, mat, theraven
Approved by: bde (src committer), mat (mentor)
Differential Revision: https://reviews.freebsd.org/D9491

MFC r313863 (by mmokhi):

Fix building of r313761 on platforms that
`long double` is alias of `double` (MIPS, etc)

PR: 216850 216851 216852 216856 216857 216858
Reported by: emsate
Reviewed by: bde emaste hselasky
Approved by: bde emaste hselasky
Differential Revision: https://reviews.freebsd.org/D9491

MFC r313864 (by mmokhi):

Add documentations related to new APIs of r313761

PR: 216850 216851 216852 216856 216857 216858
Submitted by: sgk@troutmask.apl.washington.edu
Reported by: sgk@troutmask.apl.washington.edu
Reviewed by: bde emaste hselasky
Approved by: bde emaste hselasky
Differential Revision: https://reviews.freebsd.org/D9491

MFC r314950 (by ngie):

Don't expect :test_large_inputs to fail with i386 anymore

Recent changes (maybe a side-effect of the ATF-ification in r314649)
invalidate the failure expectation.

PR: 205446
Sponsored by: Dell EMC Isilon

MFC r317349 (by pfg):

msun: Remove trailing space in Sunsoft copyright statement.

Submittedby: kargl

MFC r319047 (by mmel):

Implement sincos, sincosf, and sincosl.
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>
Differential Revision: https://reviews.freebsd.org/D10765

MFC r321457 (by ngie):

Mark :reduction as an expected failure

It fails with clang 5.0+.

PR: 220989
Reported by: Jenkins

MFC r322418 (by rlibby):

lib/msun: avoid referring to broken LDBL_MAX

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

MFC r322435 (by rlibby):

Revert r322418, LDBL_MAX_EXP unsuitable for macro pasting on some arches

Either need a different way to spell HALF_LDBL_MAX, or a different way
to spell LDBL_MAX_EXP, or a different approach.

Reported by: ian

MFC r322921 (by ngie):

Revert r321457

It doesn't fail after ^/head@r322855 (the releng_50 clang merge).

PR: 220989

Copy head@r302406 to stable/11 as part of the 11.0-RELEASE cycle.
Prune svn:mergeinfo from the new branch, as nothing has been merged
here.

Additional commits post-branch will follow.

Approved by: re (implicit)
Sponsored by: The FreeBSD Foundation

The value small=2**-(p+3), where p is the precision, can be determine from
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

For targets that have a signed zero, lgamma_r(-0, &signgamp) should
set signgamp = -1.

Submitted by: enh at google dot com (e_lgamma[f]_r.c)

* Makefile:
. 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

Remove an initialized, but otherwise, unused variable.

According to the ISO C standard, lgamma(-integer) returns
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

Fix a tab that somehow became 8 spaces.

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

Compute sin(pi*x) without actually doing the pi*x multiplication.
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

Remove some unnecessary initializations.

Obtained from: DragonFlyBSD

s/rcsid/__FBSDID/

Don't assume that int is signed 32-bits in one place. Keep assuming
that ints have >= 31 value bits elsewhere. s/int/int32_t/ seems to
have been done too globally for all other files in msun/src before
msun/ was imported into FreeBSD.

Minor fixes in comments.

e_lgamma_r.c:
Describe special cases in more detail:
- exception for lgamma(0) and lgamma(neg.integer)
- lgamma(-Inf) = Inf. This is wrong but is required by C99 Annex F. I
hope to change this.

Fixed about 50 million errors of infinity ulps and about 3 million errors
of between 1.0 and 1.8509 ulps for lgammaf(x) with x between -2**-21 and
-2**-70.

As usual, the cutoff for tiny args was not correctly translated to
float precision. It was 2**-70 but 2**-21 works. Not as usual, having
a too-small threshold was worse than a pessimization. It was just a
pessimization for (positive) args between 2**-70 and 2**-21, but for
the first ~50 million (negative) args below -2**-70, the general code
overflowed and gave a result of infinity instead of correct (finite)
results near 70*log(2). For the remaining ~361 million negative args
above -2**21, the general code gave almost-acceptable errors (lgamma[f]()
is not very accurate in general) but the pessimization was larger than
for misclassified tiny positive args.

Now the max error for lgammaf(x) with |x| < 2**-21 is 0.7885 ulps, and
speed and accuracy are almost the same for positive and negative args
in this range. The maximum error overall is still infinity ulps.

A cutoff of 2**-70 is probably wastefully small for the double precision
case. Smaller cutoffs can be used to reduce the max error to nearly
0.5 ulps for tiny args, but this is useless since the general algrorithm
for nearly-tiny args is not nearly that accurate -- it has a max error of
about 1 ulp.

Use only double precision for "kernel" cosf and sinf (except for
returning float). The functions are renamed from __kernel_{cos,sin}f()
to __kernel_{cos,sin}df() so that misuses of them will cause link errors
and not crashes.

This version is an almost-routine translation with no special optimizations
for accuracy or efficiency. The not-quite-routine part is that in
__kernel_cosf(), regenerating the minimax polynomial with double
precision coefficients gives a coefficient for the x**2 term that is
not quite -0.5, so the literal 0.5 in the code and the related `hz'
variable need to be modified; also, the special code for reducing the
error in 1.0-x**2*0.5 is no longer needed, so it is convenient to
adjust all the logic for the x**2 term a little. Note that without
extra precision, it would be very bad to use a coefficient of other
than -0.5 for the x**2 term -- the old version depends on multiplication
by -0.5 being infinitely precise so as not to need even more special
code for reducing the error in 1-x**2*0.5.

This gives an unimportant increase in accuracy, from ~0.8 to ~0.501
ulps. Almost all of the error is from the final rounding step, since
the choice of the minimax polynomials so that their contribution to the
error is a bit less than 0.5 ulps just happens to give contributions that
are significantly less (~.001 ulps).

An Athlons, for uniformly distributed args in [-2pi, 2pi], this gives
overall speed increases in the 10-20% range, despite giving a speed
decrease of typically 19% (from 31 cycles up to 37) for sinf() on args
in [-pi/4, pi/4].

Fix formatting, this is hard to explain, so I'll show one example.

- float ynf(int n, float x) /* wrapper ynf */
+float
+ynf(int n, float x) /* wrapper ynf */

This is because the __STDC__ stuff was indented.

Reviewed by: md5

Assume __STDC__, remove non-__STDC__ code.

Reviewed by: md5

$Id$ -> $FreeBSD$

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 long.

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

Remove trailing whitespace.

This commit was generated by cvs2svn to compensate for changes in r2116,
which included commits to RCS files with non-trunk default branches.

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