Remove svn:mergeinfo carried over from stable/9.

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

Copy stable/9 to releng/9.3 as part of the 9.3-RELEASE cycle.

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

MFC r262613:

Merge the projects/clang-sparc64 branch back to head. This brings in
several updates from the llvm and clang trunks to make the sparc64
backend fully functional.

Apart from one patch to sys/sparc64/include/pcpu.h which is still under
discussion, this makes it possible to let clang fully build world and
kernel for sparc64.

Any assistance with testing this on actual sparc64 hardware is greatly
appreciated, as there will unavoidably be bugs left.

Many thanks go to Roman Divacky for his upstream work on getting the
sparc64 backend into shape.

MFC r262985:

Repair a few minor mismerges from r262261 in the clang-sparc64 project
branch. This is also to minimize differences with upstream.

MFC r242879:

Only define isnan, isnanf, __isnan and __isnanf in libc.so, not in
libc.a and libc_p.a. In addition, define isnan in libm.a and libm_p.a,
but not in libm.so.

This makes it possible to statically link executables using both isnan
and isnanf with libc and libm.

Tested by: kargl

MFC r242894:

Add an explanatory comment to lib/libc/gen/isnan.c about the fix to make
static linking with libc and libm work.

Requested by: jilles

MFC r239192:

Change a few extern inline functions in libm to static inline, since
they need to refer to static constants, which C99 does not allow for
extern inline functions.

While here, change a comment in e_rem_pio2f.c to mention the correct
number of bits.

Reviewed by: bde

MFC r239195:

Add __always_inline to __ieee754_rem_pio2() and __ieee754_rem_pio2f(),
since some older versions of gcc refuse to inline these otherwise.

Requested by: bde

Merge r236148. This allows C++ code to include <cmath> after explicitly including math.h.

MFH r234685: utf8 and drop middle name

Merge quick_exit and changes required for C++11 code to compile against FreeBSD headers.

Merges changes from: r227472 r227475 r227475 r227476 r227476 r227490 r227490 r228322 r228323 r228329 r228330 r228528 r228529 r228529 r228901 r228918 r228918 r232971 r232971

Also bump __FreeBSD_version for this and the xlocale merge.

MFC r235286:

General mdoc(7) and typo fixes.

PR: 167734

MFC r233973:
Fix bugs in remquo{,f,l}.

MFC r226595:
Per IEEE754r, pow(1, y) is 1 even if y is NaN, and pow(-1, +-Inf) is 1.

MFC r226594:
Bugfix: feenableexcept() and fedisableexcept() should just return the
old exception mask, not mask | ~FE_ALL_EXCEPT.

MFC various fma{,f,l} improvements:

r226245 - refactoring
r226371 - fix double-rounding bug
r226373 - new math_private.h macros
r226601 - fix nit in r226371

Reparent mergeinfo introduced in r229461 on lib/msun/man and lib/libc/gen to more appropriate directories.

This is an intentional direct commit to the stable/9 branch.

Approved by: lstewart

MFC r227458, r226436:

- change "is is" to "is" or "it is"
- change "the the" to "the"
- other typo fixes

Approved by: lstewart

Copy head to stable/9 as part of 9.0-RELEASE release cycle.

Approved by: re (implicit)

In the libm access macros for the double type, z can sometimes
be used uninitialized. This can lead to spurious exceptions
and bit clobbering.

Submitted by: bde
Approved by: das (mentor)

Fix clang warnings.

Approved by: philip (mentor)

Clean up the unneeded cpp macro INLINE_REM_PIO2L.

Reviewed by: das
Approved by: das (mentor)

Improve the accuracy from a max ULP of ~2000 to max ULP < 0.79
on i386-class hardware for sinl and cosl. The hand-rolled argument
reduction have been replaced by e_rem_pio2l() implementations. To
preserve history the following commands have been executed:

svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h
mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h

svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h
mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h

The ld80 version has been tested by bde, das, and kargl over the
last few years (bde, das) and few months (kargl). An older ld128
version was tested by das. The committed version has only been
compiled tested via 'make universe'.

Approved by: das (mentor)
Obtained from: bde

Take two. Add the missing file that should have been committed
with r219571 and re-enable building of cbrtl.

Implement the long double version for the cube root function, cbrtl.
The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.

Reviewed by: bde
Approved by: das

Temporary disable the building of cbrtl until I
can determine why svn will not allow one to commit
a new file.

Approved by: das (implicit)

Implement the long double version for the cube root function, cbrtl.
The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.

Reviewed by: bde
Approved by: das

Add cexp() to the complex(3) manpage. Thanks to bde for pointing out
that I missed this.

Remove part of an uncommitted change that snuck into the last commit.

Convert log10f() to use __kernel_log(), which is more accurate and simpler.

Convert log10() to use __kernel_log(), which is more accurate and simpler.

Add cexp() and cexpf().

Reviewed by: bde (earlier version)

Fix typos - remove duplicate "the".

PR: bin/154928
Submitted by: Eitan Adler <lists at eitanadler.com>
MFC after: 3 days

Add complex(3) manual page documenting our partial support for C99
complex arithmetic in libm.

Reviewed by: David Schultz <das@FreeBSD.org>
MFC after: 2 weeks

Fix a bug where the wrong argument was passed to SET_FLOAT_WORD().
This bug results in a type mismatch that happens to be harmless
because of the way SET_FLOAT_WORD() works.

Submitted by: bde

Fix a bug where the wrong argument was passed to INSERT_WORDS().
This bug results in a type mismatch that happens to be harmless
because of the way INSERT_WORDS() works.

Submitted by: bde

For small arguments, these functions use simple approximations,
e.g. cos(small) = 1, sin(small) = small. This commit tightens
the thresholds at which the simple approximations are used.

Reviewed by: bde

Fix a bogus threshold that was copied from the double precision version.
This commit should have no effect on correctness; it merely changes the
threshold at which a simpler approximation can be used.

Reviewed by: bde

Remove duplicate .note.GNU-stack section declaration.

Reported by: arundel

Add section .note.GNU-stack for assembly files used by 386 and amd64.

Another minor nit: Make sure the constant here is a float so the compiler
doesn't promote the entire expression to double.

Fix various nits in style and comments that were pointed out by bde.
Code changes verified with md5.

Add log2() and log2f().

Add a "kernel" log function, based on e_log.c, which is useful for
implementing accurate logarithms in different bases. This is based
on an approach bde coded up years ago.

This function should always be inlined; it will be used in only a few
places, and rudimentary tests show a 40% performance improvement in
implementations of log2() and log10() on amd64.

The kernel takes a reduced argument x and returns the same polynomial
approximation as e_log.c, but omitting the low-order term. The low-order
term is much larger than the rest of the approximation, so the caller of
the kernel function can scale it to the appropriate base in extra precision
and obtain a much more accurate answer than by using log(x)/log(b).

Disable gcc's built-in rint() function when compiling s_nearbyint.c.
It results in incorrect optimizations that break nearbyint().

PR: 143358
Reviewed by: bde

Fix bug in jn(3) and jnf(3) that led to -inf results

Explanation by Steve:
jn[f](n,x) for certain ranges of x uses downward recursion to compute
the value of the function. The recursion sequence that is generated is
proportional to the actual desired value, so a normalization step is
taken. This normalization is j0[f](x) divided by the zeroth sequence
member. As Bruce notes, near the zeros of j0[f](x) the computed value
can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only
the leading decimal digit is correct. The solution to the issue is
fairly straight forward. The zeros of j0(x) and j1(x) never coincide,
so as j0(x) approaches a zero, the normalization constant switches to
j1[f](x) divided by the 2nd sequence member. The expectation is that
j1[f](x) is a more accurately computed value.

PR: bin/144306
Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu>
Reviewed by: bde
MFC after: 7 days

Use __FBSDID() instead of RCSID() in most .S files under lib/msun/i386,
and one under lib/msun/amd64. This avoids adding the identifiers to the
.text section, and moves them to the .comment section instead.

Suggested by: bde
Approved by: rpaulo (mentor)

This is exactly the same as the .else, so remove it.

MFtbemd: Move to using MACHINE_CPUARCH, now that it is safe.

Repair some build breakage introduced in r211725 and garbage collect some
code made obsolete in the same commit.

MFtbemd:

Prefer MACHNE_CPUARCH to MACHINE_ARCH in most contexts where you want
to test of all the CPUs of a given family conform.

Fix typos, spelling, formatting and mdoc mistakes found by Nobuyuki while
translating these manual pages. Minor corrections by me.

Submitted by: Nobuyuki Koganemaru <n-kogane@syd.odn.ne.jp>

powerpc64 floating-point is identical to powerpc, so use the same
code on both architectures.

Introduce __isnanf() as an alias for isnanf(), and make the isnan()
macro expand to __isnanf() instead of isnanf() for float arguments.
This change is needed because isnanf() isn't declared in strict POSIX
or C99 mode.

Compatibility note: Apps using isnan(float) that are compiled after
this change won't link against an older libm.

Reported by: Florian Forster <octo@verplant.org>

mdoc: spell out theta, the Unicode glyph is hard to read for terminal fonts

It is referred to as "theta" later in the document anyway,
so stop being fancy.

mdoc: spell macros correctly, there's no need for the backslash escape

mdoc: Garbage collect unused/unneeded macros

mdoc: move remaining sections into consistent order

This pertains mostly to FILES, HISTORY, EXIT STATUS and AUTHORS sections.

Found by: mdocml lint run
Reviewed by: ru

mdoc: move CAVEATS, BUGS and SECURITY CONSIDERATIONS sections to the
bottom of the manpages and order them consistently.

GNU groff doesn't care about the ordering, and doesn't even mention
CAVEATS and SECURITY CONSIDERATIONS as common sections and where to put
them.

Found by: mdocml lint run
Reviewed by: ru

Fix several typos in macros or macro misusage.

Found by: make manlint
Reviewed by: ru
Approved by: philip (mentor)

Placate new binutils, by using 16-bit %ax instead of 32-bit %eax as an
argument for fnstsw. Explicitely specify sizes for the XMM control and
status word and X87 control and status words.

Reviewed by: das
Tested by: avg
MFC after: 2 weeks

Use the documented machine constraint for SSE registers.

The amd64-specific bits of msun use an undocumented constraint, which is
less likely to be supported by other compilers (such as Clang). Change
the code to use a more common machine constraint.

Obtained from: /projects/clangbsd/

Use ISO C99 style inline semantics in msun.

Because we use ISO C99 nowadays, we can just get rid of enforcing
GNU89-style inlining.

Use, in uncovered part, the END() macro in order to improve debugging.
In this specific case, Valgrind won't get confused when analyzing such
functions.

Sponsored by: Sandvine Incorporated
Tested by: emaste
MFC: 3 days

Namespace: scalb() is withdrawn from POSIX.

Eliminate __real__ and __imag__ gccisms.

C99 TC2 now wants FP_FAST_FMA* to be defined to 1, if the macros are
defined at all. See also: defect report #223.

Use __gnu89_inline so that these files will compile with newer versions
of gcc, where the meaning of 'inline' was changed to match C99.

Noticed by: rdivacky

Fix the types of INFINITY and NAN, which were broken in r131851. They
should both be floats, not doubles.

PR: 127795
Submitted by: Christoph Mallon
MFC after: 2 weeks

Add support for the FPA floating-point format on ARM. The
FPA floating-point format is identical to the VFP format,
but is always stored in big-endian.
Introduce _IEEE_WORD_ORDER to describe the byte-order of
the FP representation.

Obtained from: Juniper Networks, Inc

Clean out some empty mergeinfo records, presumably by people doing local
cp/mv operations. The full repo-relative URL should be specified for the
source in these cases.

Remove some unused variables.

Reported by: Intel C Compiler

In the line
#pragma STDC CX_LIMITED_RANGE ON
the "ON" needs to be in caps. gcc doesn't understand this pragma
anyway and assumes it is always on in any case, but icc supports
it and cares about the case.

Implement cproj{,f,l}().

Use cpack() and the gcc extension __imag__ to implement cimag() and
conj() instead of using expressions like z * I. The latter is bad for
several reasons:

1. It is implemented using arithmetic, which is unnecessary, and can
generate floating point exceptions, contrary to the requirements on
these functions.

2. gcc implements complex multiplication using a formula that breaks
down for infinities, e.g., it gives INFINITY * I == nan + inf I.

Fix some style bogosity from fdlibm.

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

Submitted by: bde

A few minor corrections, including some from bde:
- When y/x is huge, it's faster and more accurate to return pi/2
instead of pi - pi/2.
- There's no need for 3 lines of bit fiddling to compute -z.
- Fix a comment.

On i386, gcc truncates long double constants to double precision
at compile time regardless of the dynamic precision, and there's
no way to disable this misfeature at compile time. Hence, it's
impossible to generate the appropriate tables of constants for the
long double inverse trig functions in a straightforward way on i386;
this change hacks around the problem by encoding the underlying bits
in the table.

Note that these functions won't pass the regression test on i386,
even with the FPU set to extended precision, because the regression
test is similarly damaged by gcc. However, the tests all pass when
compiled with a modified version of gcc.

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

Add implementations of acosl(), asinl(), atanl(), atan2l(),
and cargl().

Reviewed by: bde
sparc64 testing resources from: remko

Set WARNS=1.

I believe I've committed all the bits necessary to make this compile
on all supported architectures. :crosses fingers:

The high part of the mantissa is 64 bits on sparc64.

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.

Sort the .PATH entries to give a more reasonable order of precedence:
1. architecture-specific files
2. long double format-specific files
3. bsdsrc
4. src
5. man
The original order was virtually the opposite of this.

This should not cause any functional changes at this time. The
difference is only significant when one wants to override, say, a
generic foo.c with a more specialized foo.c (as opposed to foo.S).

Fix a typo in the cosl() prototype.

Implement fmodl.
Document fmodl and fix some errors in the fmod manpage.

Symbol.map is handled by cpp, so use C-style comments

Approved by: cognet (mentor)

Add mips support to libm, from mips2-jnpr perforce branch.

Fix some corner cases:
- fma(x, y, z) returns z, not NaN, if z is infinite, x and y are finite,
x*y overflows, and x*y and z have opposite signs.
- fma(x, y, z) doesn't generate an overflow, underflow, or inexact exception
if z is NaN or infinite, as per IEEE 754R.
- If the rounding mode is set to FE_DOWNWARD, fma(1.0, 0.0, -0.0) is -0.0,
not +0.0.

Remove a (bogus) remnant of debugging this on sparc64.

Add assembly versions of remquol() and remainderl().

Hook remquol() and remainderl() up to the build.

Implement remainderl() as a wrapper around remquol(). The extra work
remquol() performs to compute the quotient is negligible.

Implement remquol() based on remquo().

Implement csqrtl().

Hook hypotl() and cabsl() up to the build.

Document hypotl().

Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>

Alias hypotl() and cabsl() for platforms where long double is the same
as double.

Implement cabsl() in terms of hypotl().

Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>

Implement hypotl(). This is bde's conversion of fdlibm hypot(), with minor
fixes for ld128 by me.

Use fabs[f]() instead of bit fiddling for setting absolute values.
This makes little difference in float precision, but in double
precision gives a speedup of about 30% on amd64 (A64 CPU) and i386
(A64). This depends on fabs[f]() being inline and efficient. The
bit fiddling (or any use of SET_HIGH_WORD(), which libm does too
much because it was best on old 32-bit machines) always causes
packing overheads and sometimes causes stalls in the packing, since
it operates on only part of a variable in the double precision case.
It apparently did cause stalls in a critical path here.

Use the expression fabs(x+0.0)-fabs(y+0.0) instead of
fabs(x+0.0)+fabs(y+0.0) when mixing NaNs. This improves
consistency of the result by making it harder for the compiler to reorder
the operands. (FP addition is not necessarily commutative because the
order of operands makes a difference on some machines iff the operands are
both NaNs.)

Fix a missing mask in a hi+lo decomposition. Thus bug made the extra
precision in software useless, so hypotf() had some errors in the 1-2
ulp range unless there is extra precision in hardware (as happens on
i386).

Include math.h for the fmaf() prototype.

Fix some rather obscene code that has ambiguous if...if...else...
constructs in it.

1 << 47 needs to be written 1ULL << 47.

Hook up sqrtl() to the build.

MD implementations of sqrtl().

MI implementation of sqrtl(). This is very slow and should
be overridden when hardware sqrt is available.

Fix and improve some magic numbers for the "medium size" case.

e_rem_pio2.c:
This case goes up to about 2**20pi/2, but the comment about it said that
it goes up to about 2**19pi/2.

It went too far above 2**pi/2, giving a multiplier fn with 21 significant
bits in some cases. This would be harmful except for a numerical
accident. It happens that the terms of the approximation to pi/2,
when rounded to 33 bits so that multiplications by 20-bit fn's are
exact, happen to be rounded to 32 bits so multiplications by 21-bit
fn's are exact too, so the bug only complicates the error analysis (we
might lose a bit of accuracy but have bits to spare).

e_rem_pio2f.c:
The bogus comment in e_rem_pio2.c was copied and the code was changed
to be bug-for-bug compatible with it, except the limit was made 90
ulps smaller than necessary. The approximation to pi/2 was not
modified except for discarding some of it.

The same rough error analysis that justifies the limit of 2**20pi/2
for double precision only justifies a limit of 2**18pi/2 for float
precision. We depended on exhaustive testing to check the magic numbers
for float precision. More exaustive testing shows that we can go up
to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision,
with a the maximum error for cosf() and sinf() unchanged at 0.5009
ulps despite the maximum error in rem_pio2f being ~0.25 ulps. Implement
this.

Inline __ieee754__rem_pio2f(). On amd64 (A64) and i386 (A64), this
gives an average speedup of about 12 cycles or 17% for
9pi/4 < |x| <= 2**19pi/2 and a smaller speedup for larger x, and a
small speeddown for |x| <= 9pi/4 (only 1-2 cycles average, but that
is 4%).

Inlining this is less likely to bust caches than inlining the float
version since it is much smaller (about 220 bytes text and rodata) and
has many fewer branches. However, the float version was already large
due to its manual inlining of the branches and also the polynomial
evaluations.

Use a temporary array instead of the arg array y[] for calling
__kernel_rem_pio2(). This simplifies analysis of aliasing and thus
results in better code for the usual case where __kernel_rem_pio2()
is not called. In particular, when __ieee854_rem_pio2[f]() is inlined,
it normally results in y[] being returned in registers. I couldn't
get this to work using the restrict qualifier.

In float precision, this saves 2-3% in most cases on amd64 and i386
(A64) despite it not being inlined in float precision yet. In double
precision, this has high variance, with an average gain of 2% for
amd64 and 0.7% for i386 (but a much larger gain for usual cases) and
some losses.

Change __ieee754_rem_pio2f() to return double instead of float so that
this function and its callers cosf(), sinf() and tanf() don't waste time
converting values from doubles to floats and back for |x| > 9pi/4.
All these functions were optimized a few years ago to mostly use doubles
internally and across the __kernel*() interfaces but not across the
__ieee754_rem_pio2f() interface.

This saves about 40 cycles in cosf(), sinf() and tanf() for |x| > 9pi/4
on amd64 (A64), and about 20 cycles on i386 (A64) (except for cosf()
and sinf() in the upper range). 40 cycles is about 35% for |x| < 9pi/4
<= 2**19pi/2 and about 5% for |x| > 2**19pi/2. The saving is much
larger on amd64 than on i386 since the conversions are not easy to
optimize except on i386 where some of them are automatic and others
are optimized invalidly. amd64 is still about 10% slower in cosf()
and tanf() in the lower range due to conversion overhead.

This also gives a tiny speedup for |x| <= 9pi/4 on amd64 (by simplifying
the code). It also avoids compiler bugs and/or additional slowness
in the conversions on (not yet supported) machines where double_t !=
double.

Fix some off-by-1 errors.

e_rem_pio2.c:
Float and double precision didn't work because init_jk[] was 1 too small.
It needs to be 2 larger than you might expect, and 1 larger than it was
for these precisions, since its test for recomputing needs a margin of
47 bits (almost 2 24-bit units).

init_jk[] seems to be barely enough for extended and quad precisions.
This hasn't been completely verified. Callers now get about 24 bits
of extra precision for float, and about 19 for double, but only about
8 for extended and quad. 8 is not enough for callers that want to
produce extra-precision results, but current callers have rounding
errors of at least 0.8 ulps, so another 1/2**8 ulps of error from the
reduction won't affect them much.

Add a comment about some of the magic for init_jk[].

e_rem_pio2.c:
Double precision worked in practice because of a compensating off-by-1
error here. Extended precision was asked for, and it executed exactly
the same code as the unbroken double precision.

e_rem_pio2f.c:
Float precision worked in practice because of a compensating off-by-1
error here. Double precision was asked for, and was almost needed,
since the cosf() and sinf() callers want to produce extra-precision
results, at least internally so that their error is only 0.5009 ulps.
However, the extra precision provided by unbroken float precision is
enough, and the double-precision code has extra overheads, so the
off-by-1 error cost about 5% in efficiency on amd64 and i386.

Let PowerPC world optionally build with -msoft-float. For FPU-less PowerPC
variations (e500 currently), this provides a gcc-level FPU emulation and is an
alternative approach to the recently introduced kernel-level emulation
(FPU_EMU).

Approved by: cognet (mentor)
MFp4: e500

Optimize the 9pi/2 < |x| <= 2**19pi/2 case some more by avoiding an
fabs(), a conditional branch, and sign adjustments of 3 variables for
x < 0 when the branch is taken. In double precision, even when the
branch is perfectly predicted, this saves about 10 cycles or 10% on
amd64 (A64) and i386 (A64) for the negative half of the range, but
makes little difference for the positive half of the range. In float
precision, it also saves about 4 cycles for the positive half of the
range on i386, and many more cycles in both halves on amd64 (28 in the
negative half and 11 in the positive half for tanf), but the amd64
times for float precision are anomalously slow so the larger
improvement is only a side effect.

Previous commits arranged for the x < 0 case to be handled simply:
- one part of the rounding method uses the magic number 0x1.8p52
instead of the usual 0x1.0p52. The latter is required for large |x|,
but it doesn't work for negative x and we don't need it for large |x|.
- another part of the rounding method no longer needs to add `half'.
It would have needed to add -half for negative x.
- removing the "quick check no cancellation" in the double precision
case removed the need to take the absolute value of the quadrant
number.

Add my noncopyright in e_rem_pio2.c

Avoid using FP-to-integer conversion for !(amd64 || i386) too. Use the
FP-to-FP method to round to an integer on all arches, and convert this
to an int using FP-to-integer conversion iff irint() is not available.
This is cleaner and works well on at least ia64, where it saves 20-30
cycles or about 10% on average for 9Pi/4 < |x| <= 32pi/2 (should be
similar up to 2**19pi/2, but I only tested the smaller range).

After the previous commit to e_rem_pio2.c removed the "quick check no
cancellation" non-optimization, the result of the FP-to-integer
conversion is not needed so early, so using irint() became a much
smaller optimization than when it was committed.

An earlier commit message said that cos, cosf, sin and sinf were equally
fast on amd64 and i386 except for cos and sin on i386. Actually, cos
and sin on amd64 are equally fast to cosf and sinf on i386 (~88 cycles),
while cosf and sinf on amd64 are not quite equally slow to cos and sin
on i386 (average 115 cycles with more variance).

Remove the "quick check no cancellation" optimization for
9pi/2 < |x| < 32pi/2 since it is only a small or negative optimation
and it gets in the way of further optimizations. It did one more
branch to avoid some integer operations and to use a different
dependency on previous results. The branches are fairly predictable
so they are usually not a problem, so whether this is a good
optimization depends mainly on the timing for the previous results,
which is very machine-dependent. On amd64 (A64), this "optimization"
is a pessimization of about 1 cycle or 1%; on ia64, it is an
optimization of about 2 cycles or 1%; on i386 (A64), it is an
optimization of about 5 cycles or 4%; on i386 (Celeron P2) it is an
optimization of about 4 cycles or 3% for cos but a pessimization of
about 5 cycles for sin and 1 cycle for tan. I think the new i386
(A64) slowness is due to an pipeline stall due to an avoidable
load-store mismatch (so the old timing was better), and the i386
(Celeron) variance is due to its branch predictor not being too good.

Optimize the 9pi/2 < |x| <= 2**19pi/2 case on amd64 and i386 by avoiding
the the double to int conversion operation which is very slow on these
arches. Assume that the current rounding mode is the default of
round-to-nearest and use rounding operations in this mode instead of
faking this mode using the round-towards-zero mode for conversion to
int. Round the double to an integer as a double first and as an int
second since the double result is needed much earler.

Double rounding isn't a problem since we only need a rough approximation.
We didn't support other current rounding modes and produce much larger
errors than before if called in a non-default mode.

This saves an average about 10 cycles on amd64 (A64) and about 25 on
i386 (A64) for x in the above range. In some cases the saving is over
25%. Most cases with |x| < 1000pi now take about 88 cycles for cos
and sin (with certain CFLAGS, etc.), except on i386 where cos and sin
(but not cosf and sinf) are much slower at 111 and 121 cycles respectivly
due to the compiler only optimizing well for float precision. A64
hardware cos and sin are slower at 105 cycles on i386 and 110 cycles
on amd64.

Add an irint() function in inline asm for amd64 and i386. irint() is
the same as lrint() except it returns int instead of long. Though the
extern lrint() is fairly fast on these arches, it still takes about
12 cycles longer than the inline version, and 12 cycles is a lot in
applications where [li]rint() is used to avoid slow conversions that
are only a couple of times slower.

This is only for internal use. The libm versions of *rint*() should
also be inline, but that would take would take more header engineering.
Implementing irint() instead of lrint() also avoids a conflict with
the extern declaration of the latter.

Optimize the conversion to bits a little (by about 11 cycles or 16%
on i386 (A64), 5 cycles on amd64 (A64), and 3 cycles on ia64). gcc
tends to generate very bad code for accessing floating point values
as bits except when the integer accesses have the same width as the
floating point values, and direct accesses to bit-fields (as is common
only for long double precision) always gives such accesses. Use the
expsign access method, which is good for 80-bit long doubles and
hopefully no worse for 128-bit long doubles. Now the generated code
is less bad. There is still unnecessary copying of the arg on amd64
and i386 and mysterious extra slowness on amd64.

Optimize the fixup for +-0 by using better classification for this case
and by using a table lookup to avoid a branch when this case occurs.
On i386, this saves 1-4 cycles out of about 64 for non-large args.

Fix rintl() on signaling NaNs and unsupported formats.

s/rcsid/__FBSDID/

Remove an unused variable.

Eliminate some warnings.

Merge cosmetic changes from e_rem_pio2.c 1.10 (convert to __FBSDID();
fix indentation and return type of __ieee754_rem_pio2()).

Remove unused variables.

Optimize for 3pi/4 <= |x| <= 9pi/4 in much the same way as for
pi/4 <= |x| <= 3pi/4. Use the same branch ladder as for float precision.
Remove the optimization for |x| near pi/2 and don't do it near the
multiples of pi/2 in the newly optimized range, since it requires
fairly large code to handle only relativley few cases. Ifdef out
optimization for |x| <= pi/4 since this case can't occur because it
is done in callers.

On amd64 (A64), for cos() and sin() with uniformly distributed args,
no cache misses, some parallelism in the caller, and good but not great
CC and CFLAGS, etc., this saves about 40 cycles or 38% in the newly
optimized range, or about 27% on average across the range |x| <= 2pi
(~65 cycles for most args, while the A64 hardware fcos and fsin take
~75 cycles for half the args and 125 cycles for the other half). The
speedup for tan() is much smaller, especially relatively. The speedup
on i386 (A64) is slightly smaller, especially relatively. i386 is
still much slower than amd64 here (unlike in the float case where it
is slightly faster).

Rearrange the polynomial evaluation for better parallelism. This
saves an average of about 8 cycles or 5% on A64 (amd64 and i386 --
more in cycles but about the same percentage on i386, and more with
old versions of gcc) with good CFLAGS and some parallelism in the
caller. As usual, it takes a couple more multiplications so it will
be slower on old machines.

Convert to __FBSDID().

Document return values better.

Add tgammaf() as a simple wrapper around tgamma().

2 long double constants were missing L suffixes. This helped break tanl()
on !(amd64 || i386). It gave slightly worse than double precision in some
cases. tanl() now passes tests of 2^24 values on ia64.

Fix a typo which broke k_tanl.c on !(amd64 || i386).

Inline __ieee754__rem_pio2(). With gcc4-2, this gives an average
optimization of about 10% for cos(x), sin(x) and tan(x) on
|x| < 2**19*pi/2. We didn't do this before because __ieee754__rem_pio2()
is too large and complicated for gcc-3.3 to inline very well. We don't
do this for float precision because it interferes with optimization
of the usual (?) case (|x| < 9pi/4) which is manually inlined for float
precision only.

This has some rough edges:
- some static data is duplicated unnecessarily. There isn't much after
the recent move of large tables to k_rem_pio2.c, and some static data
is duplicated to good affect (all the data static const, so that the
compiler can evaluate expressions like 2*pio2 at compile time and
generate even more static data for the constant for this).
- extern inline is used (for the same reason as in previous inlining of
k_cosf.c etc.), but C99 apparently doesn't allow extern inline
functions with static data, and gcc will eventually warn about this.

Convert to __FBSDID().

Indent __ieee754_rem_pio2()'s declaration consistently (its style was
made inconsistent with fdlibm a while ago, so complete this).

Fix __ieee754_rem_pio2()'s return type to match its prototype. Someone
changed too many ints to int32_t's when fixing the assumption that all
ints are int32_t's.

Use volatile hacks to make sure exp() generates an underflow
exception when it's supposed to. Previously, gcc -O2 was optimizing
away the statement that generated it.

Hook up sinl(), cosl(), and tanl() to the build.

Add implementations of sinl(), cosl(), and tanl().

Submitted by: Steve Kargl <sgk@apl.washington.edu>

Documentation for sinl(), cosl(), and tanl().

Add kernel functions for 128-bit long doubles. These could be improved
a bit, but access to a freebsd/sparc64 machine is needed.

Submitted by: bde and Steve Kargl <sgk@apl.washington.edu> (earlier version)

Add kernel functions for 80-bit long doubles. Many thanks to Steve and
Bruce for putting lots of effort into these; getting them right isn't
easy, and they went through many iterations.

Submitted by: Steve Kargl <sgk@apl.washington.edu> with revisions from bde

Add more pi for long doubles. Also, avoid storing multiple copies
of the pi/2 array, as it is unlikely to vary, except in Indiana.

Sigh, the weak reference for ceill(), floorl() and truncl() was in
unreachable code due to a missing include. This kept arm and powerpc
broken.

Reported by: sam, grehan

Oops, the weak reference for ceill(), floorl() and truncl() was in the
wrong file. This broke arm and powerpc.

Reported by: grehan

Use the expression fabs(x+0.0)+fabs(y+0.0) instad of a+b (where a is
|x| or |y| and b is |y| or |x|) when mixing NaN arg(s).

hypot*() had its own foot shooting for mixing NaNs -- it swaps the
args so that |x| in bits is largest, but does this before quieting
signaling NaNs, so on amd64 (where the result of adding NaNs depends
on the order) it gets inconsistent results if setting the quiet bit
makes a difference, just like a similar ia64 and i387 hardware comparison.
The usual fix (see e_powf.c 1.13 for more details) of mixing using
(a+0.0)+-(b+0.0) doesn't work on amd64 if the args are swapped (since
the rder makes a difference with SSE). Fortunately, the original args
are unchanged and don't need to be swapped when we let the hardware
decide the mixing after quieting them, but we need to take their
absolute value.

hypotf() doesn't seem to have any real bugs masked by this non-bug.
On amd64, its maximum error in 2^32 trials on amd64 is now 0.8422 ulps,
and on i386 the maximum error is unchanged and about the same, except
with certain CFLAGS it magically drops to 0.5 (perfect rounding).

Convert to __FBSDID().

Fix the hi+lo decomposition for 2/(3ln2). The decomposition needs to
be into 12+24 bits of precision for extra-precision multiplication,
but was into 13+24 bits. On i386 with -O1 the bug was hidden by
accidental extra precision, but on amd64, in 2^32 trials the bug
caused about 200000 errors of more than 1 ulp, with a maximum error
of about 80 ulps. Now the maximum error in 2^32 trials on amd64
is 0.8573 ulps. It is still 0.8316 ulps on i386 with -O1.

The nearby decomposition of 1/ln2 and the decomposition of 2/(3ln2) in
the double precision version seem to be sub-optimal but not broken.

Use the expression (x+0.0)-(y+0.0) instead of x+y when mixing NaN arg(s).
This uses 2 tricks to improve consistency so that more serious problems
aren't hidden in simple regression tests by noise for the NaNs:

- for a signaling NaN, adding 0.0 generates the invalid exception and
converts to a quiet NaN, and doesn't have too many effects for other
types of args (it converts -0 to +0 in some rounding modes, but that
hopefully doesn't change the result after adding the NaN arg). This
avoids some inconsistencies on i386 and ia64. On these arches, the
result of an operation on 2 NaNs is apparently the largest or the
smallest of the NaNs as bits (consistently largest or smallest for
each arch, but the opposite). I forget which way the comparison
goes and if the sign bit affects it. The quiet bit is is handled
poorly by not always setting it before the comparision or ignoring
it. Thus if one of the args was originally a signaling NaN and the
other was originally a quiet NaN, then the result depends too much
on whether the signaling NaN has been quieted at this point, which
in turn depends on optimizations and promotions. E.g., passing float
signaling NaNs to double functions must quiet them on conversion;
on i387, loading a signaling NaN of type float or double (but not
long double) into a register involves a conversion, so it quiets
signaling NaNs, so if the addition has 2 register operands than it
only sees quiet NaNs, but if the addition has a memory operand then
it sees a signaling NaN iff it is in the memory operand.

- subtraction instead of addition is used to avoid a dubious optimization
in old versions of gcc. For SSE operations, mixing of NaNs apparently
always gives the target operand. This is not as good as the i387
and ia64 behaviour. It doesn't mix NaNs at all, and makes addition
not quite commutative. Old versions of gcc sometimes rewrite x+y
to y+x and thus give different results (in bits) for NaNs. gcc-3.3.3
rewrites x+y to y+x for one of pow() and powf() but not the other,
so starting from float NaN args x and y, powf(x, y) was almost always
different from pow(x, y).

These tricks won't give consistency of 2-arg float and double functions
with long double ones on amd64, since long double ones use the i387
which has different semantics from SSE.

Convert to __FBSDID().

s_ceill.c
s_floorl.c
s_truncl.c

On arches where long double is the same as double, alias ceil(), floor()
and trunc() to the corresponding long double functions. This is not
just an optimization for these arches. The full long double functions
have a wrong value for `huge', and the arches without full long doubles
depended on it being wrong.

Fix the C version of ceill(x) for -1 < x <= -0 in all rounding modes.
The result should be -0, but was +0.

Fix exp2*(x) on signaling NaNs by returning x+x as usual.

This has the side effect of confusing gcc-4.2.1's optimizer into more
often doing the right thing. When it does the wrong thing here, it
seems to be mainly making too many copies of x with dependency chains.
This effect is tiny on amd64, but in some cases on i386 it is enormous.
E.g., on i386 (A64) with -O1, the current version of exp2() should
take about 50 cycles, but took 83 cycles before this change and 66
cycles after this change. exp2f() with -O1 only speeded up from 51
to 47 cycles. (exp2f() should take about 40 cycles, on an Athlon in
either i386 or amd64 mode, and now takes 42 on amd64). exp2l() with
-O1 slowed down from 155 cycles to 123 for some args; this is unimportant
since the i386 exp2l() is a fake; the wrong thing for it seems to
involve branch misprediction.

Rearrange the polynomial evaluation for better parallelism. This is
faster on all machines tested (old Celeron (P2), A64 (amd64 and i386)
and ia64) except on ia64 when compiled with -O1. It takes 2 more
multiplications, so it will be slower on old machines. The speedup
is about 8 cycles = 17% on A64 (amd64 and i386) with best CFLAGS
and some parallelism in the caller.

Move the evaluation of 2**k up a bit so that it doesn't compete too
much with the new polynomial evaluation. Unlike the previous
optimization, this rearrangement cannot change the result, so compilers
and CPU schedulers can do it, but they don't do it quite right yet.
This saves a whole 1 or 2 cycles on A64.

Use hardware remainder on amd64 since it is 5 to 10 times faster than
software remainder and is already used for remquo().

Fix remainder() and remainderf() in round-towards-minus-infinity mode
when the result is +-0. IEEE754 requires (in all rounding modes) that
if the result is +-0 then its sign is the same as that of the first
arg, but in round-towards-minus-infinity mode an uncorrected implementation
detail always reversed the sign. (The detail is that x-x with x's
sign positive gives -0 in this mode only, but the algorithm assumed
that x-x always has positive sign for such x.)

remquo() and remquof() seem to need the same fix, but I cannot test them
yet.

Use long doubles when mixing NaN args. This trick improves consistency
of results on at least amd64, so that more serious problems like the
above aren't hidden in simple regression tests by noise for the NaNs.
On amd64, hardware remainder should be used since it is about 10 times
faster than software remainder and is already used for remquo(), but
it involves using the i387 even for floats and doubles, and the i387
does NaN mixing which is better than but inconsistent with SSE NaN mixing.
Software remainder() would probably have been inconsistent with
software remainderl() for the same reason if the latter existed.

Signaling NaNs cause further inconsistencies on at least ia64 and i386.

Use __FBSDID().

Use double precision for z and thus for the entire calculation of
exp2(i/TBLSIZE) * p(z) instead of only for the final multiplication
and addition. This fixes the code to match the comment that the maximum
error is 0.5010 ulps (except on machines that evaluate float expressions
in extra precision, e.g., i386's, where the evaluation was already
in extra precision).

Fix and expand the comment about use of double precision.

The relative roundoff error from evaluating p(z) in non-extra precision
was about 16 times larger than in exp2() because the interval length
is 16 times smaller. Its maximum was at least P1 * (1.0 ulps) *
max(|z|) ~= log(2) * 1.0 * 1/32 ~= 0.0217 ulps (1.0 ulps from the
addition in (1 + P1*z) with a cancelation error when z ~= -1/32). The
actual final maximum was 0.5313 ulps, of which 0.0303 ulps must have
come from the additional roundoff error in p(z). I can't explain why
the additional roundoff error was almost 3/2 times larger than the rough
estimate.

As usual, use a minimax polynomial that is specialized for float
precision. The new polynomial has degree 4 instead of 10, and a maximum
error of 2**-30.04 ulps instead of 2**-33.15. This doesn't affect the
final error significantly; the maximum error was and is about 0.5015
ulps on i386 -O1, and the number of cases with an error of > 0.5 ulps
is increased from 13851 to 14407.

Note that the error is only this close to 0.5 ulps due to excessive
extra precision caused by compiler bugs on i386. The extra precision
could be obtained intentionally, and is useful for keeping the error
of the hyperbolic float functions below 1 ulp, since these functions
are implemented using expm1f. My recent change for scaling by 2**k
had the unintentional side effect of retaining extra precision for
longer, so callers of expm1f see errors of more like 0.0015 ulps than
0.5015 ulps, and for the hyperbolic functions this reduces the maximum
error from nearly about 2 ulps to about 0.75 ulps.

This is about 10% faster on i386 (A64). expm1* is still very slow,
but now the float version is actually significantly faster. The
algorithm is very sophisticated but not very good except on machines
with fast division.

Fix a comment about coefficients and expand a related one.

Fix truncl() when the result should be -0.0L. When the result is +-0.0L,
it must have the same sign as the arg in all rounding modes, but it was
always +0.0L.

Oops, fix the fix in rev.1.10. logb() and logbf() were broken on
denormals, and logb() remained broken after 1.10 because the fix for
logbf() was incompletely translated.

Convert to __FBSDID().

Use a better method of scaling by 2**k. Instead of adding to the
exponent bits of the reduced result, construct 2**k (hopefully in
parallel with the construction of the reduced result) and multiply by
it. This tends to be much faster if the construction of 2**k is
actually in parallel, and might be faster even with no parallelism
since adjustment of the exponent requires a read-modify-wrtite at an
unfortunate time for pipelines.

In some cases involving exp2* on amd64 (A64), this change saves about
40 cycles or 30%. I think it is inherently only about 12 cycles faster
in these cases and the rest of the speedup is from partly-accidentally
avoiding compiler pessimizations (the construction of 2**k is now
manually scheduled for good results, and -O2 doesn't always mess this
up). In most cases on amd64 (A64) and i386 (A64) the speedup is about
20 cycles. The worst case that I found is expf on ia64 where this
change is a pessimization of about 10 cycles or 5%. The manual
scheduling for plain exp[f] is harder and not as tuned.

Details specific to expm1*:
- the saving is closer to 12 cycles than to 40 for expm1* on i386 (A64).
For some reason it is much larger for negative args.
- also convert to __FBSDID().

Use a better method of scaling by 2**k. Instead of adding to the
exponent bits of the reduced result, construct 2**k (hopefully in
parallel with the construction of the reduced result) and multiply by
it. This tends to be much faster if the construction of 2**k is
actually in parallel, and might be faster even with no parallelism
since adjustment of the exponent requires a read-modify-wrtite at an
unfortunate time for pipelines.

In some cases involving exp2* on amd64 (A64), this change saves about
40 cycles or 30%. I think it is inherently only about 12 cycles faster
in these cases and the rest of the speedup is from partly-accidentally
avoiding compiler pessimizations (the construction of 2**k is now
manually scheduled for good results, and -O2 doesn't always mess this
up). In most cases on amd64 (A64) and i386 (A64) the speedup is about
20 cycles. The worst case that I found is expf on ia64 where this
change is a pessimization of about 10 cycles or 5%. The manual
scheduling for plain exp[f] is harder and not as tuned.

This change ld128/s_exp2l.c has not been tested.

As for the float trig functions and logf, use a minimax polynomial
that is specialized for float precision. The new polynomial has degree
5 instead of 11, and a maximum error of 2**-27.74 ulps instead
of 2**-30.64. This doesn't affect the final error significantly; the
maximum error was and is about 0.9101 ulps on amd64 -01 and the number
of cases with an error of > 0.5 ulps is actually reduced by epsilon
despite the larger error in the polynomial.

This is about 15% faster on amd64 (A64), i386 (A64) and ia64. The asm
version is still used instead of this on i386 since it is faster and
more accurate.

Adjust the exponent before converting the result from double to
float precision. This fixes some double rounding problems for
subnormals and simplifies things a bit.

Fix a harmless type error in 1.9.

Fix cutoffs. This is just a cleanup and an optimization for unusual
cases which are used mainly by regression tests.

As usual, the cutoff for tiny args was not correctly translated to
float precision. It was 2**-54 but 2**-24 works. It must be about
2**-precision, since the error from approximating log(1+x) by x is
about the same as |x|. Exhaustive testing shows that 2**-24 gives
perfect rounding in round-to-nearest mode.

Similarly for the cutoff for being small, except this is not used by
so many other functions. It was 2**-29 but 2**-15 works. It must be
a bit smaller than sqrt(2**-precision), since the error from
approximating log(1+x) by x-x*x/2 is about the same as x*x. Exhaustive
testing shows that 2**-15 gives a maximum error of 0.5052 ulps in
round-to-nearest-mode. The algorithm for the general case is only good
for 0.8388 ulps, so this is sufficient (but it loses slightly on i386 --
then extra precision gives 0.5032 ulps for the general case).

While investigating this, I noticed that optimizing the usual case by
falling into a middle case involving a simple polynomial evaluation
(return x-x*x/2 instead of x here) is not such a good idea since it
gives an enormous pessimization of tinier args on machines for which
denormals are slow. Float x*x/2 is denormal when |x| ~< 2**-64 and
x*x/2 is evaluated in float precision, so it can easily be denormal
for normal x. This is even more interesting for general polynomial
evaluations. Multiplying out large powers of x is normally a good
optimization since it reduces dependencies, but it creates denormals
starting with quite large x.

Oops, when merging from the float version to the double versions, don't
forget to translate "float" to "double".

ucbtest didn't detect the bug, but exhaustive testing of the float
case relative to the double case eventually did. The bug only affects
args x with |x| ~> 2**19*(pi/2) on non-i386 (i386 is broken in a
different way for large args).

Remove the float version of the kernel of arg reduction for pi/2, since
it should never have existed and it has not been used for many years
(floats are reduced faster using doubles). All relevant changes (just
the workaround for broken assignment) have been merged to the double
version.

Do an ordinary assignment in STRICT_ASSIGN() except for floats until
there is a problem with non-floats (when i386 defaults to extra
precision). This essentially restores yesterday's behaviour for doubles
on i386 (since generic rint() isn't used and everywhere else assumed
working assignment), but for arches that use the generic rint() it
finishes restoring some of 1995's behaviour (don't waste time doing
unnecessary store/load).

Use STRICT_ASSIGN() for exp2f() and exp2() instead of a volatile
variable hack for exp2f() only.

The volatile variable had a surprisingly large cost for exp2f() -- 19
cycles or 15% on i386 in the worst case observed. This is only partly
explained by there being several references to the variable, only one
of which benefited from it being volatile. Arches that have working
assignment are likely to benefit even more from not having any volatile
variable.

exp2() now has a chance of working with extra precision on i386.

exp2() has even more references to the variable, so it would have been
pessimized more by simply declaring the variable as volatile. Even
the temporary volatile variable for STRICT_ASSIGN costs 5-10% on i386,
(A64) so I will change STRICT_ASSIGN() to do an ordinary assignment
until i386 defaults to extra precision.

Use STRICT_ASSIGN() for _kernel_rem_pio2f() and _kernel_rem_pio2f()
instead of a volatile cast hack for the float version only. The cast
hack broke with gcc-4, but this was harmless since the float version
hasn't been used for a few years. Merge from the float version so
that the double version has a chance of working on i386 with extra
precision.

See k_rem_pio2f.c rev.1.8 for the original hack.

Convert to _FBSDID().

Use STRICT_ASSIGN() for log1pf() and log1p() instead of a volatile cast
hack for log1pf() only. The cast hack broke with gcc-4, resulting in
~1 million errors of more than 1 ulp, with a maximum error of ~1.5 ulps.
Now the maximum error for log1pf() on i386 is 0.5034 ulps again (this
depends on extra precision), and log1p() has a chance of working with
extra precision.

See s_log1pf.c 1.8 for the original hack. (It claims only 62343 large
errors).

Convert to _FBSDID(). Another thing broken with gcc-4 is the static
const hack used for rcsids.

Use STRICT_ASSIGN() instead of assorted direct volatile hacks to work
around assignments not working for gcc on i386. Now volatile hacks
for rint() and rintf() don't needlessly pessimize so many arches
and the remaining pessimizations (for arm and powerpc) can be avoided
centrally.

This cleans up after s_rint.c 1.3 and 1.13 and s_rintf.c 1.3 and 1.9:
- s_rint.c 1.13 broke 1.3 by only using a volatile cast hack in 1 place
when it was needed in 2 places, and the volatile cast hack stopped
working with gcc-4. These bugs only affected correctness tests on
i386 since i386 normally uses asm rint() and doesn't support the
extra precision mode that would break assignments of doubles.
- s_rintf.c 1.9 improved(?) on 1.3 by using a volatile variable hack
instead of an extra-precision variable hack, but it declared 2
variables as volatile when only 1 variable needed to be volatile.
This only affected speed tests on i386 since i386 uses asm rintf().

Use volatile hacks to make sure these functions generate an underflow
exception when they're supposed to. Previously, gcc -O2 was optimizing
away the statement that generated it.

Hook up exp2l() and related docs to the build.

Introduce a new log(3) manpage and move the relevant functions there.
Document exp2l() in exp(3), and remove the quaint discussion of topics
such as what these functions were called on the HP-71B's variant of
BASIC.

Implement exp2l(). There is one version for machines with 80-bit
long doubles (i386, amd64, ia64) and one for machines with 128-bit
long doubles (sparc64). Other platforms use the double version.
I've only done runtime testing on i386.

Thanks to bde@ for helpful discussions and bugfixes.

Add a macro STRICT_ASSIGN() to help avoid the compiler bug that
assignments and casts don't clip extra precision, if any. The
implementation is to assign to a temporary volatile variable and read
the result back to assign to the original lvalue.

lib/msun currently 2 different hard-coded hacks to avoid the problem
in just a few places and needs it in a few more places. One variant
uses volatile for the original lvalue. This works but is slower than
necessary. Another temporarily casts the lvalue to volatile. This
broke with gcc-4.2.1 or earlier (gcc now stores to the lvalue but
doesn't load from it).

Optimize this a bit better.

Submitted by: bde (although these aren't all of his changes)

Implement rintl(), nearbyintl(), lrintl(), and llrintl().
Thanks to bde@ for feedback and testing of rintl().

- Correct the range check in the double version to catch negative values
that would overflow.
- Style fixes and improved handling of NaNs suggested by bde.

Grumble. DO declare logbl(), DON'T declare logl() just yet.
bde is going to commit logl() Real Soon Now.
I'm just trying to slow him down with merge conflicts.

Noticed by: bde

Remove the declaration of logl(). The relevant bits haven't been
committed yet, but the declaration leaked in when I added nan() and
friends.

Reported by: pav

Since nan() is supposed to work the same as strtod("nan(...)", NULL),
my original implementation made both use the same code. Unfortunately,
this meant libm depended on a vendor header at compile time and previously-
unexposed vendor bits in libc at runtime.

Hence, I just wrote my own version of the relevant vendor routine. As it
turns out, mine has a factor of 8 fewer of lines of code, and is a bit more
readable anyway. The strtod() and *scanf() routines still use vendor code.

Reviewed by: bde

Remove z_abs(). The z_*() functions were in libf77, and for some reason
someone thought it would be a good idea to copy z_abs() to libm in 1994.
However, it's never been declared or documented anywhere, and I'm
reasonably confident that nobody uses it.

Discussed with: bde, deischen, kan

Oops, the previous commit was not needed -- the file was committed but
not checked out due to my checkout error.

Translate from the i386 so that this compiles and runs.

I hope that this and the i386 version of it will not be needed, but
this is currently about 16 cycles or 36% faster than the C version,
and the i386 version is about 8 cycles or 19% faster than the C
version, due to poor optimization of the C version.

Don't try to build s_nanl.c before it is committed.

Add logbl(3) to libm.

Document the fact that we have nan(3) now, and make some minor clarifications
in other places.

Implement and document nan(), nanf(), and nanl(). This commit
adds two new directories in msun: ld80 and ld128. These are for
long double functions specific to the 80-bit long double format
used on x86-derived architectures, and the 128-bit format used on
sparc64, respectively.

1. Add csqrt{,f}(3).
2. Put carg{,f}(3) under the FBSD_1.1 namespace where it belongs
(requested by kan@)

Implement and document csqrt(3) and csqrtf(3).

Update the standards section, and make a minor clarification about the
return value of sqrt.

Typo in previous commit

Symbol.map additions for carg and cargf. (They're in C99, so I didn't
add a new version for them.)

s/C90/C99/

Add a "STANDARDS" section.

Implement carg(3) and cargf(3).

Rotting in an old src tree since: March 2005

Oops, back out previous commit since it was backwards to a wrong branch.

MFC: 1.11: fix the threshold for (not) using the simple Taylor approximation.

Fix an aliasing bug which was finally detected by gcc-4.2. fdlibm has
hundreds of similar aliasing bugs, but all except this one seem to have
been fixed by Cygnus and/or NetBSD before the modified version of fdlibm
was imported into FreeBSD in 1994.

PR: standards/113147
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>

Merge the relevant part of rev.1.14 of s_cbrt.c (a micro-optimization
involving moving the check for x == 0). The savings in cycles are
smaller for cbrtf() than for cbrt(), and positive in all measured cases
with gcc-3.4.4, but still very machine/compiler-dependent.

Bump library versions in preparation for 7.0.

Ok'd by: kan

Enable symbol versioning by default. Use WITHOUT_SYMVER to disable it.
Warning, after symbol versioning is enabled, going back is not easy
(use WITHOUT_SYMVER at your own risk).

Change the default thread library to libthr.

There most likely still needs to be a version bump for at least the
thread libraries. If necessary, this will happen later.

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.

Fix tgamma() on some special args:
(1) tgamma(-Inf) returned +Inf and failed to raise any exception, but
should always have raised an exception, and should behave like
tgamma(negative integer).
(2) tgamma(negative integer) returned +Inf and raised divide-by-zero,
but should return NaN and raise "invalid" on any IEEEish system.
(3) About half of the 2**52 negative intgers between -2**53 and -2**52
were misclassified as non-integers by using floor(x + 0.5) to round
to nearest, so tgamma(x) was wrong (+-0 instead of +Inf and now NaN)
on these args. The floor() expression is hard to use since rounding
of (x + 0.5) may give x or x + 1, depending on |x| and the current
rounding mode. The fixed version uses ceil(x) to classify x before
operating on x and ends up being more efficient since ceil(x) is
needed anyway.
(4) On at least the problematic args in (3), tgamma() raised a spurious
inexact.
(5) tgamma(large positive) raised divide-by-zero but should raise overflow.
(6) tgamma(+Inf) raised divide-by-zero but should not raise any exception.
(7) Raise inexact for tiny |x| in a way that has some chance of not being
optimized away.

The fix for (5) and (6), and probably for (2), also prevents -O optimizing
away the exception.

PR: 112180 (2)
Standards: Annex F in C99 (IEC 60559 binding) requires (1), (2) and (6).

Document (in a comment) the current (slightly broken) handling of special
values in more detail, and change the style of this comment to be closer
to fdlibm and C99:
- tgamma(-Inf) was undocumented and is wrong (+Inf, should be NaN)
- tgamma(negative integer) is as intended (+Inf) but not best for IEEE-754
(NaN)
- tgamma(-0) was documented as being wrong (+Inf) but was correct (-Inf)
- documentation of setting of exceptions (overflow, etc.) was more
complete here than in most of libm, but was further from matching
the actual setting than in most of libm, due to various bugs here
(primarily, always evaluating +Inf one/zero and getting unwanted
divide-by-zero exceptions from this). Now the actual behaviour with
gcc -O0 is documented. Optimization still breaks setting of exceptions
all over libm, so nothing can depend on this working.
- tgamma(NaN)'s exception was documented as being wrong (invalid) but was
correct (no exception with IEEEish NaNs).

Finish (?) rev.1.5. gamma was not renamed to tgamma in one place.

Finish (?) rev.1.6. errno.h was not completely removed.

Use C comments since we now preprocess these files with CPP.

Remove California Regent's clause 3, per letter

Implement modfl().

Fix a problem relating to fesetenv() clobbering i387 register stack.

Details: As a side-effect of restoring a saved FP environment,
fesetenv() overwrites the tag word, which indicates which i387
registers are in use. Normally this isn't a problem because
the calling convention requires the register stack to be empty
on function entry and exit. However, fesetenv() is inlined, so we
need to tell gcc explicitly that the i387 registers get clobbered.

PR: 85101

Fix a cut-and-paste-o.

Correctly handle NaN.

Correctly handle inf/nan. This routine is currently unused because we
seem to have assembly versions for all architectures, but it can't
hurt to fix it.

Remove modf from libm's symbol map. It's actually in libc for
hysterical raisins.

Remove an unneeded fnstcw instruction.

Noticed by: bde

Remove a note pertaining to the Alpha.

Moved __BEGIN_DECLS up a little so that it covers __test_sse() and C++
isn't broken,

PR: 104425

Remove alpha left-overs.

Fixed the threshold for using the simple Taylor approximation.

In e_log.c, there was just a off-by-1 (1 ulp) error in the comment
about the threshold. The precision of the threshold is unimportant,
but the magic numbers in the code are easier to understand when the
threshold is described precisely.

In e_logf.c, mistranslation of the magic numbers gave an off-by-1
(1 * 16 ulps) error in the intended negative bound for the threshold
and an off-by-7 (7 * 16 ulps) error in the intended positive bound for
the threshold, and the intended bounds were not translated from the
double precision bounds so they were unnecessarily small by a factor
of about 2048.

The optimization of using the simple Taylor approximation for args
near a power of 2 is dubious since it only applies to a relatively
small proportion of args, but if it is done then doing it 2048 times
as often _may_ be more efficient. (My benchmarks show unexplained
dependencies on the data that increase with further optimizations
in this area.)

Fixed tanh(-0.0) on ia64 and optimizeed tanh(x) for 2**-55 <= |x| <
2**-28 as a side effect, by merging with the float precision version
of tanh() and the double precision version of sinh().

For tiny x, tanh(x) ~= x, and we used the expression x*(one+x) to
return this value (x) and set the inexact flag iff x != 0. This
doesn't work on ia64 since gcc -O does the dubious optimization
x*(one+x) = x+x*x so as to use fma, so the sign of -0.0 was lost.

Instead, handle tiny x in the same as sinh(), although this is imperfect:
- return x directly and set the inexact flag in a less efficient way.
- increased the threshold for non-tinyness from 2**-55 to 2**-28 so that
many more cases are optimized than are pessimized.

Updated some comments and fixed bugs in others (ranges for half-open
intervals mostly had the open end backwards, and there were nearby style
bugs).

Removed the optimized asm versions of scalb() and scalbf(). These
functions are only for compatibility with obsolete standards. They
shouldn't be used, so they shouldn't be optimized. Use the generic
versions instead.

This fixes scalbf() as a side effect. The optimized asm version left
garbage on the FP stack. I fixed the corresponding bug in the optimized
asm scalb() and scalbn() in 1996. NetBSD fixed it in scalb(), scalbn()
and scalbnf() in 1999 but missed fixing it in scalbf(). Then in 2005
the bug was reimplemented in FreeBSD by importing NetBSD's scalbf().

The generic versions have slightly different error handling:
- the asm versions blindly round the second parameter to a (floating
point) integer and proceed, while the generic versions return NaN
if this rounding changes the value. POSIX permits both behaviours
(these functions are XSI extensions and the behaviour for a bogus
non-integral second parameter is unspecified). Apart from this
and the bug in scalbf(), the behaviour of the generic versions seems
to be identical. (I only exhusatively tested
generic_scalbf(1.0F, anyfloat) == asm_scalb(1.0F, anyfloat). This
covers many representative corner cases involving NaNs and Infs but
doesn't test exception flags. The brokenness of scalbf() showed up
as weird behaviour after testing just 7 integer cases sequentially.)

Backed out rev.1.10. It tried to implement ldexpf() as a weak reference
to scalbf(), but ldexpf() cannot be implemented in that way since the
types of the second parameter differ. ldexpf() can be implemented as
a weak or strong reference to scalbnf() (*) but that was already done
long before rev.1.10 was committed. The old implementation uses a
reference, so rev.1.10 had no effect on applications. The C files for
the scalb() family are not used for amd64 or i386, so rev.1.10 had even
less effect for these arches.

(*) scalbnf() raises the radix to the given exponent, while ldexpf()
raises 2 to the given exponent. Thus the functions are equivalent
except possibly for their error handling iff the radix is 2. Standards
more or less require identical error handling. Under FreeBSD, the
functions are equivalent except for more details being missing in
scalbnf()'s man page.

Add symbol versioning to libm.

Oops, on amd64 (and probably on all non-i386 systems), the previous
commit broke the 2**24 cases where |x| > DBL_MAX/2. There are exponent
range problems not just for denormals (underflow) but for large values
(overflow). Doubles have more than enough exponent range to avoid the
problems, but I forgot to convert enough terms to double, so there was
an x+x term which was sometimes evaluated in float precision.

Unfortunately, this is a pessimization with some combinations of systems
and compilers (it makes no difference on Athlon XP's, but on Athlon64's
it gives a 5% pessimization with gcc-3.4 but not with gcc-3.3).

Exlain the problem better in comments.

Use double precision internally to optimize cbrtf(), and change the
algorithm for the second step significantly to also get a perfectly
rounded result in round-to-nearest mode. The resulting optimization
is about 25% on Athlon64's and 30% on Athlon XP's (about 25 cycles
out of 100 on the former).

Using extra precision, we don't need to do anything special to avoid
large rounding errors in the third step (Newton's method), so we can
regroup terms to avoid a division, increase clarity, and increase
opportunities for parallelism. Rearrangement for parallelism loses
the increase in clarity. We end up with the same number of operations
but with a division reduced to a multiplication.

Using specifically double precision, there is enough extra precision
for the third step to give enough precision for perfect rounding to
float precision provided the previous steps are accurate to 16 bits.
(They were accurate to 12 bits, which was almost minimal for imperfect
rounding in the old version but would be more than enough for imperfect
rounding in this version (9 bits would be enough now).) I couldn't
find any significant time optimizations from optimizing the previous
steps, so I decided to optimize for accuracy instead. The second step
needed a division although a previous commit optimized it to use a
polynomial approximation for its main detail, and this division dominated
the time for the second step. Use the same Newton's method for the
second step as for the third step since this is insignificantly slower
than the division plus the polynomial (now that Newton's method only
needs 1 division), significantly more accurate, and simpler. Single
precision would be precise enough for the second step, but doesn't
have enough exponent range to handle denormals without the special
grouping of terms (as in previous versions) that requires another
division, so we use double precision for both the second and third
steps.

Extract the high and low words together. With gcc-3.4 on uniformly
distributed non-large args, this saves about 14 of 134 cycles for
Athlon64s and about 5 of 199 cycles for AthlonXPs.

Moved the check for x == 0 inside the check for subnormals. With
gcc-3.4 on uniformly distributed non-large args, this saves another
5 cycles on Athlon64s and loses 1 cycle on AthlonXPs.

Use INSERT_WORDS() and not SET_HIGH_WORD() when converting the first
approximation from bits to double. With gcc-3.4 on uniformly distributed
non-large args, this saves another 4 cycles on both Athlon64s and and
AthlonXPs.

Accessing doubles as 2 words may be an optimization on old CPUs, but on
current CPUs it tends to cause extra operations and pipeline stalls,
especially for writes, even when only 1 of the words needs to be accessed.

Removed an unused variable.

Use a minimax polynomial approximation instead of a Pade rational
function approximation for the second step. The polynomial has degree
2 for cbrtf() and 4 for cbrt(). These degrees are minimal for the final
accuracy to be essentially the same as before (slightly smaller).
Adjust the rounding between steps 2 and 3 to match. Unfortunately,
for cbrt(), this breaks the claimed accuracy slightly although incorrect
rounding doesn't. Claim less accuracy since its not worth pessimizing
the polynomial or relying on exhaustive testing to get insignificantly
more accuracy.

This saves about 30 cycles on Athlons (mainly by avoiding 2 divisions)
so it gives an overall optimization in the 10-25% range (a larger
percentage for float precision, especially in 32-bit mode, since other
overheads are more dominant for double precision, surprisingly more
in 32-bit mode).

Fixed code to match comments and the algorithm:
- in preparing for the third approximation, actually make t larger in
magnitude than cbrt(x). After chopping, t must be incremented by 2
ulps to make it larger, not 1 ulp since chopping can reduce it by
almost 1 ulp and it might already be up to half a different-sized-ulp
smaller than cbrt(x). I have not found any cases where this is
essential, but the think-time error bound depends on it. The relative
smallness of the different-sized-ulp limited the bug. If there are
cases where this is essential, then the final error bound would be
5/6+epsilon instead of of 4/6+epsilon ulps (still < 1).
- in preparing for the third approximation, round more carefully (but
still sloppily to avoid branches) so that the claimed error bound of
0.667 ulps is satisfied in all cases tested for cbrt() and remains
satisfied in all cases for cbrtf(). There isn't enough spare precision
for very sloppy rounding to work:
- in cbrt(), even with the inadequate increment, the actual error was
0.6685 in some cases, and correcting the increment increased this
a little. The fix uses sloppy rounding to 25 bits instead of very
sloppy rounding to 21 bits, and starts using uint64_t instead of 2
words for bit manipulation so that rounding more bits is not much
costly.
- in cbrtf(), the 0.667 bound was already satisfied even with the
inadequate increment, but change the code to almost match cbrt()
anyway. There is not enough spare precision in the Newton
approximation to double the inadequate increment without exceeding
the 0.667 bound, and no spare precision to avoid this problem as
in cbrt(). The fix is to round using an increment of 2 smaller-ulps
before chopping so that an increment of 1 ulp is enough. In cbrt(),
we essentially do the same, but move the chop point so that the
increment of 1 is not needed.

Fixed comments to match code:
- in cbrt(), the second approximation is good to 25 bits, not quite 26 bits.
- in cbrt(), don't claim that the second approximation may be implemented
in single precision. Single precision cannot handle the full exponent
range without minor but pessimal changes to renormalize, and although
single precision is enough, 25 bit precision is now claimed and used.

Added comments about some of the magic for the error bound 4/6+epsilon.
I still don't understand why it is 4/6+ and not 6/6+ ulps.

Indent comments at the right of code more consistently.

Added comments about the apparently-magic rational function used in
the second step of approximating cbrt(x). It turns out to be neither
very magic not nor very good. It is just the (2,2) Pade approximation
to 1/cbrt(r) at r = 1, arranged in a strange way to use fewer operations
at a cost of replacing 4 multiplications by 1 division, which is an
especially bad tradeoff on machines where some of the multiplications
can be done in parallel. A Remez rational approximation would give
at least 2 more bits of accuracy, but the (2,2) Pade approximation
already gives 6 more bits than needed. (Changed the comment which
essentially says that it gives 3 more bits.)

Lower order Pade approximations are not quite accurate enough for
double precision but are plenty for float precision. A lower order
Remez rational approximation might be enough for double precision too.
However, rational approximations inherently require an extra division,
and polynomial approximations work well for 1/cbrt(r) at r = 1, so I
plan to switch to using the latter. There are some technical
complications that tend to cost a division in another way.

Optimize by not doing excessive conversions for handling the sign bit.
This gives an optimization of between 9 and 22% on Athlons (largest
for cbrt() on amd64 -- from 205 to 159 cycles).

We extracted the sign bit and worked with |x|, and restored the sign
bit as the last step. We avoided branches to a fault by using accesses
to FP values as bits to clear and restore the sign bit. Avoiding
branches is usually good, but the bit access macros are not so good
(especially for setting FP values), and here they always caused pipeline
stalls on Athlons. Even using branches would be faster except on args
that give perfect branch misprediction, since only mispredicted branches
cause stalls, but it possible to avoid touching the sign bit in FP
values at all (except to preserve it in conversions from bits to FP
not related to the sign bit). Do this. The results are identical
except in 2 of the 3 unsupported rounding modes, since all the
approximations use odd rational functions so they work right on strictly
negative values, and the special case of -0 doesn't use an approximation.

Fixed some especially horrible style bugs (indentation that is neither
KNF nor fdlibmNF combined with multiple statements per line).

Added comments about the magic behind
<cbrt(x) in bits> ~= <x in bits>/3 + BIAS.
Keep the large comments only in the double version as usual.

Fixed some style bugs (mainly grammar and spelling errors in comments).

Fixed the unexpectedly large maximum error after the previous commit.
It was because I forgot to translate the part of the double precision
algorithm that chops t so that t*t is exact. Now the maximum error
is the same as for double precision (almost exactly 2.0/3 ulps).

Fixed all 502518670 errors of more than 1 ulp for cbrtf() on amd64.
The maximum error was 3.56 ulps.

The bug was another translation error. The double precision version
has a comment saying "new cbrt to 23 bits, may be implemented in
precision". This means exactly what it says -- that the 23 bit second
approximation for the double precision cbrt() may be implemented in
single (i.e., float) precision. It doesn't mean what the translation
assumed -- that this approximation, when implemented in float precision,
is good enough for the the final approximation in float precision.
First, float precision needs a 24 bit approximation. The "23 bit"
approximation is actually good to 24 bits on float precision args, but
only if it is evaluated in double precision. Second, the algorithm
requires a cleanup step to ensure its error bound.

In float precision, any reasonable algorithm works for the cleanup
step. Use the same algorithm as for double precision, although this
is much more than enough and is a significant pessimization, and don't
optimize or simplify anything using double precision to implement the
float case, so that the whole double precision algorithm can be verified
in float precision. A maximum error of 0.667 ulps is claimed for cbrt()
and the max for cbrtf() using the same algorithm shouldn't be different,
but the actual max for cbrtf() on amd64 is now 0.9834 ulps. (On i386
-O1 the max is 0.5006 (down from < 0.7) due to extra precision.)

Fixed some magic numbers.

The threshold for not being tiny was too small. Use the usual 2**-12
threshold. As for sinhf, use a different method (now the same as for
sinhf) to set the inexact flag for tiny nonzero x so that the larger
threshold works, although this method is imperfect. As for sinhf,
this change is not just an optimization, since the general code that
we fell into has accuracy problems even for tiny x. On amd64, avoiding
it fixes tanhf on 2*13495596 args with errors of between 1 and 1.3
ulps and thus reduces the total number of args with errors of >= 1 ulp
from 37533748 to 5271278; the maximum error is unchanged at 2.2 ulps.

The magic number 22 is log(DBL_MAX)/2 plus slop. This is bogus for
float precision. Use 9 (log(FLT_MAX)/2 plus less slop than for
double precision). Unlike for coshf and tanhf, this is just an
optimization, and MAX isn't misspelled EPSILON in the commit log.

I started testing with nonstandard rounding modes, and verified that
the chosen thresholds work for all modes modulo problems not related
to thresholds. The best thresholds are not very dependent on the mode,
at least for tanhf.

"Create" ldexpf for non-i386 architectures.

Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>

Fixed the approximation to pio4. pio4_hi must be pio2_hi/2 since it
shares its low half with pio2_hi. pio2_hi is rounded down although
rounding to nearest would be a tiny bit better, so pio4_hi must be
rounded down too. It was rounded to nearest, which happens to be
different in float precision but the same in double precision.

This fixes about 13.5 million errors of more than 1 ulp in asinf().
The largest error was 2.81 ulps on amd64 and 2.57 ulps on i386 -O1.
Now the largest error is 0.93 ulps on amd65 and 0.67 ulps on i386 -O1.

For log1pf(), fixed the approximations to sqrt(2), sqrt(2)-1 and
sqrt(2)/2-1. For log1p(), fixed the approximation to sqrt(2)/2-1.

The end result is to fix an error of 1.293 ulps in
log1pf(0.41421395540 (hex 0x3ed413da))
and an error of 1.783 ulps in
log1p(-0.292893409729003961761) (hex 0x12bec4 00000001)).
The former was the only error of > 1 ulp for log1pf() and the latter
is the only such error that I know of for log1p().

The approximations don't need to be very accurate, but the last 2 need
to be related to the first and be rounded up a little (even more than
1 ulp for sqrt(2)/2-1) for the following implementation-detail reason:
when the arg (x) is not between (the approximations to) sqrt(2)/2-1
and sqrt(2)-1, we commit to using a correction term, but we only
actually use it if 1+x is between sqrt(2)/2 and sqrt(2) according to
the first approximation. Thus we must ensure that
!(sqrt(2)/2-1 < x < sqrt(2)-1) implies !(sqrt(2)/2 < x+1 < sqrt(2)),
where all the sqrt(2)'s are really slightly different approximations
to sqrt(2) and some of the "<"'s are really "<="'s. This was not done.

In log1pf(), the last 2 approximations were rounded up by about 6 ulps
more than needed relative to a good approximation to sqrt(2), but the
actual approximation to sqrt(2) was off by 3 ulps. The approximation
to sqrt(2)-1 ended up being 4 ulps too small, so the algoritm was
broken in 4 cases. The result happened to be broken in 1 case. This
is fixed by using a natural approximation to sqrt(2) and derived
approximations for the others.

In logf(), all the approximations made sense, but the approximation
to sqrt(2)/2-1 was 2 ulps too small (a tiny amount, since we compare
with a granularity of 2**32 ulps), so the algorithm was broken in 2
cases. The result was broken in 1 case. This is fixed by rounding
up the approximation to sqrt(2)/2-1 by 2**32 ulps, so 2**32 cases are
now handled a little differently (still correctly according to my
assertion that the approximations don't need to be very accurate, but
this has not been checked).

Use the usual volatile hack to trick gcc into clipping any extra precision
on assignment.

Extra precision on i386's broke hi+lo decomposition in the usual way.
It caused all except 1 of the 62343 errors of more than 1 ulp for
log1pf() on i386's with gcc -O [-fno-float-store].

Fixed fdlibm[+cygnus] logbf() and logb() on denormals. Adjustment
according to the highest nonzero bit in a denormal was missing.

fdlibm ilogbf() and ilogb() have always had the adjustment, but only
use a small part of their method for handling denormals; use the
normalization method in log[f]() for the main part.

Restored removal of the special handling needed for a result of +-0.
It was lost in rev.1.9. The log message for rev.1.9 says that the
special case of +-0 is handled twice, but it was only handled once,
so it became unhandled, and this happened to break half of the cases
that return +-0:
- round-towards-minus-infinity: 0 < x < 1: result was -0 not 0
- round-to-nearest: -0.5 <= x < 0: result was 0 not -0
- round-towards-plus-infinity: -1 < x < 0: result was 0 not -0
- round-towards-zero: -1 < x < 0: result was 0 not -0

Simplified the fix in rev.1.3. Instead of using long double for
TWO52[sx] to trick gcc into correctly converting TWO52[sx]+x to double
on assignment to "double w", force a correct assignment by assigning
to *(double *)&w. This is cleaner and avoids the double rounding
problem on machines that evaluate double expressions in double
precision. It is not necessary to convert w-TWO52[sx] to double
precision on return as implied in the comment in rev.1.3, since
the difference is exact.

Fixed rint(x) in the following cases:
(1) In round-to-nearest mode, on all machines, fdlibm rint() never
worked for |x| = n+0.75 where n is an even integer between 262144
and 524286 inclusive (2*131072 cases). To avoid double rounding
on some machines, we begin by adjusting x to a value with the 0.25
bit not set, essentially by moving the 0.25 bit to a lower bit
where it works well enough as a guard, but we botched the adjustment
when log2(|x|) == 18 (2*2**52 cases) and ended up just clearing
the 0.25 bit then. Most subcases still worked accidentally since
another lower bit serves as a guard. The case of odd n worked
accidentally because the rounding goes the right way then. However,
for even n, after mangling n+0.75 to 0.5, rounding gives n but the
correct result is n+1.
(2) In round-towards-minus-infinity mode, on all machines, fdlibm rint()
never for x = n+0.25 where n is any integer between -524287 and
-262144 inclusive (262144 cases). In these cases, after mangling
n+0.25 to n, rounding gives n but the correct result is n-1.
(3) In round-towards-plus-infinity mode, on all machines, fdlibm rint()
never for x = n+0.25 where n is any integer between 262144 and
524287 inclusive (262144 cases). In these cases, after mangling
n+0.25 to n, rounding gives n but the correct result is n+1.

A variant of this bug was fixed for the float case in rev.1.9 of s_rintf.c,
but the analysis there is incomplete (it only mentions (1)) and the fix
is buggy.

Example of the problem with double rounding: rint(1.375) on a machine
which evaluates double expressions with just 1 bit of extra precision
and is in round-to-nearest mode. We evaluate the result using
(double)(2**52 + 1.375) - 2**52. Evaluating 2**52 + 1.375 in (53+1) bit
prcision gives 2**52 + 1.5 (first rounding). (Second) rounding of this
to double gives 2**52 + 2.0. Subtracting 2**52 from this gives 2.0 but
we want 1.0. Evaluating 2**52 + 1.375 in double precision would have
given the desired intermediate result of 2**52 + 1.0.

The double rounding problem is relatively rare, so the botched adjustment
can be fixed for most machines by removing the entire adjustment. This
would be a wrong fix (using it is 1 of the bugs in rev.1.9 of s_rintf.c)
since fdlibm is supposed to be generic, but it works in the following cases:
- on all machines that evaluate double expressions in double precision,
provided either long double has the same precision as double (alpha,
and i386's with precision forced to double) or my earlier fix to use
a long double 2**52 is modified to avoid using long double precision.
- on all machines that evaluate double expressions in many more than 11
bits of extra precision. The 1 bit of extra precision in the example
is the worst case. With N bits of extra precision, it sufices to
adjust the bit N bits below the 0.5 bit. For N >= about 52 there is
no such bit so the adjustment is both impossible and unnecessary. The
fix in rev.1.9 of s_rintf.c apparently depends on corresponding magic
in float precision: on all supported machines N is either 0 or >= 24,
so double rounding doesn't occur in practice.
- on all machines that don't use fdlibm rint*() (i386's).
So under FreeBSD, the double rounding problem only affects amd64 now, but
should only affect i386 in future (when double expressions are evaluated
in long double precision).

Fixed roundf(). The following cases never worked in FreeBSD:
- in round-towards-minus-infinity mode, on all machines, roundf(x) never
worked for 0 < |x| < 0.5 (2*0x3effffff cases in all, or almost half of
float space). It was -0 for 0 < x < 0.5 and 0 for -0.5 < x < 0, but
should be 0 and -0, respectively. This is because t = ceilf(|x|) = 1
for these args, and when we adjust t from 1 to 0 by subtracting 1, we
get -0 in this rounding mode, but we want and expected to get 0.
- in round-towards-minus-infinity, round towards zero and round-to-nearest
modes, on machines that evaluate float expressions in float precision
(most machines except i386's), roundf(x) never worked for |x| =
<float value immediately below 0.5> (2 cases in all). It was +-1 but
should have been +-0. This is because t = ceilf(|x|) = 1 for these
args, and when we try to classify |x| by subtracting it from 1 we
get an unexpected rounding error -- the result is 0.5 after rounding
to float in all 3 rounding modes, so we we have forgotten the
difference between |x| and 0.5 and end up returning the same value
as for +-0.5.

The fix is to use floorf() instead of ceilf() and to add 1 instead of
-1 in the adjustment. With floorf() all the expressions used are
always evaluated exactly so there are no rounding problems, and with
adjustments of +1 we don't go near -0 when adjusting.

Attempted to fix round() and roundl() by cloning the fix for roundf().
This has only been tested for round(), only on args representable as
floats. Double expressions are evaluated in double precision even on
i386's, so round(0.5-epsilon) was broken even on i386's. roundl()
must be completely broken on i386's since long double precision is not
really supported. There seem to be no other dependencies on the
precision.

Rearranged the polynomial evaluation to reduce dependencies, as in
k_tanf.c but with different details.

The polynomial is odd with degree 13 for tanf() and odd with degree
9 for sinf(), so the details are not very different for sinf() -- the
term with the x**11 and x**13 coefficients goes awaym and (mysteriously)
it helps to do the evaluation of w = z*z early although moving it later
was a key optimization for tanf(). The details are different but simpler
for cosf() because the polynomial is even and of lower degree.

On Athlons, for uniformly distributed args in [-2pi, 2pi], this gives
an optimization of about 4 cycles (10%) in most cases (13% for sinf()
on AXP, but 0% for cosf() with gcc-3.3 -O1 on AXP). The best case
(sinf() with gcc-3.4 -O1 -fcaller-saves on A64) now takes 33-39 cycles
(was 37-45 cycles). Hardware sinf takes 74-129 cycles. Despite
being fine tuned for Athlons, the optimization is even larger on
some other arches (about 15% on ia64 (pluto2) and 20% on alpha (beast)
with gcc -O2 -fomit-frame-pointer).

Fixed cosf(x) when x is a "negative" NaNs. I broke this in rev.1.10.
cosf(x) is supposed to return something like x when x is a NaN, and
we actually fairly consistently return x-x which is normally very like
x (on i386 and and it is x if x is a quiet NaN and x with the quiet bit
set if x is a signaling NaN. Rev.1.10 broke this by normalising x to
fabsf(x). It's not clear if fabsf(x) is should preserve x if x is a NaN,
but it actually clears the sign bit, and other parts of the code depended
on this.

The bugs can be fixed by saving x before normalizing it, and using the
saved x only for NaNs, and using uint32_t instead of int32_t for ix
so that negative NaNs are not misclassified even if fabsf() doesn't
clear their sign bit, but gcc pessimizes the saving very well, especially
on Athlon XPs (it generates extra loads and stores, and mixes use of
the SSE and i387, and this somehow messes up pipelines). Normalizing
x is not a very good optimization anyway, so stop doing it. (It adds
latency to the FPU pipelines, but in previous versions it helped except
for |x| <= 3pi/4 by simplifying the integer pipelines.) Use the same
organization as in s_sinf.c and s_tanf.c with some branches reordered.
These changes combined recover most of the performance of the unfixed
version on A64 but still lose 10% on AXP with gcc-3.4 -O1 but not with
gcc-3.3 -O1.

Fixed the hi+lo approximation to log(2). The normal 17+24 bit decomposition
that was used doesn't work normally here, since we want to be able to
multiply `hi' by the exponent of x _exactly_, and the exponent of x has
more than 7 significant bits for most denormal x's, so the multiplication
was not always exact despite a cloned comment claiming that it was. (The
comment is correct in the double precision case -- with the normal 33+53
bit decomposition the exponent can have 20 significant bits and the extra
bit for denormals is only the 11th.)

Fixing this had little or no effect for denormals (I think because
more precision is inherently lost for denormals than is lost by roundoff
errors in the multiplication).

The fix is to reduce the precision of the decomposition to 16+24 bits.
Due to 2 bugs in the old deomposition and numerical accidents, reducing
the precision actually increased the precision of hi+lo. The old hi+lo
had about 39 bits instead of at least 41 like it should have had.
There were off-by-1-bit errors in each of hi and lo, apparently due
to mistranslation from the double precision hi and lo. The correct
16 bit hi happens to give about 19 bits of precision, so the correct
hi+lo gives about 43 bits instead of at least 40. The end result is
that expf() is now perfectly rounded (to nearest) except in 52561 cases
instead of except in 67027 cases, and the maximum error is 0.5013 ulps
instead of 0.5023 ulps.

Rearranged the polynomial evaluation some more to reduce dependencies.
Instead of echoing the code in a comment, try to describe why we split
up the evaluation in a special way.

The new optimization is mostly to move the evaluation of w = z*z later
so that everything else (except z = x*x) doesn't have to wait for w.
On Athlons, FP multiplication has a latency of 4 cycles so this
optimization saves 4 cycles per call provided no new dependencies are
introduced. Tweaking the other terms in to reduce dependencies saves
a couple more cycles in some cases (more on AXP than on A64; up to 8
cycles out of 56 altogether in some cases). The previous version had
a similar optimization for s = z*x. Special optimizations like these
probably have a larger effect than the simple 2-way vectorization
permitted (but not activated by gcc) in the old version, since 2-way
vectorization is not enough and the polynomial's degree is so small
in the float case that non-vectorizable dependencies dominate.

On an AXP, tanf() on uniformly distributed args in [-2pi, 2pi] now
takes 34-55 cycles (was 39-59 cycles).

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.

Exploit skew-symmetry to avoid an operation: -sin(x-A) = sin(A-x). This
gives a tiny but hopefully always free optimization in the 2 quadrants
to which it applies. On Athlons, it reduces maximum latency by 4 cycles
in these quadrants but has usually has a smaller effect on total time
(typically ~2 cycles (~5%), but sometimes 8 cycles when the compiler
generates poor code).

Try to use the "proximity" (~) operator consistently in comments
(x ~<= a, not x <= ~a). This got messed up in some places when the
comments were moved from e_rem_pio2f.c.

Added my (non-)copyright.

Changed spelling of the request-to-inline macro name to match the change
of the function name.

Added my (non-)copyright.

In k_tanf.c, added the first set of redundant parentheses to control
grouping which was claimed to be added in the previous commit.

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

Minor cleanups and optimizations:

- Remove dead code that I forgot to remove in the previous commit.

- Calculate the sum of the lower terms of the polynomial (divided by
x**5) in a single expression (sum of odd terms) + (sum of even terms)
with parentheses to control grouping. This is clearer and happens to
give better instruction scheduling for a tiny optimization (an
average of about ~0.5 cycles/call on Athlons).

- Calculate the final sum in a single expression with parentheses to
control grouping too. Change the grouping from
first_term + (second_term + sum_of_lower_terms) to
(first_term + second_term) + sum_of_lower_terms. Normally the first
grouping must be used for accuracy, but extra precision makes any
grouping give a correct result so we can group for efficiency. This
is a larger optimization (average 3-4 cycles/call or 5%).

- Use parentheses to indicate that the C order of left to right evaluation
is what is wanted (for efficiency) in a multiplication too.

The old fdlibm code has several optimizations related to these. 2
involve doing an extra operation that can be done almost in parallel
on some superscalar machines but are pessimizations on sequential
machines. Others involve statement ordering or expression grouping.
All of these except the ordering for the combining the sums of the odd
and even terms seem to be ideal for Athlons, but parallelism is still
limited so all of these optimizations combined together with the ones
in this commit save only ~6-8 cycles (~10%).

On an AXP, tanf() on uniformly distributed args in [-2pi, 2pi] now
takes 39-59 cycles. I don't know of any more optimizations for tanf()
short of writing it all in asm with very MD instruction scheduling.
Hardware fsin takes 122-138 cycles. Most of the optimizations for
tanf() don't work very well for tan[l](). fdlibm tan() now takes
145-365 cycles.

s/5.5/6.0/ in HISTORY section.

Discussed with: ru

Optimized by eliminating the special case for 0.67434 <= |x| < pi/4.

A single polynomial approximation for tan(x) works in infinite precision
up to |x| < pi/2, but in finite precision, to restrict the accumulated
roundoff error to < 1 ulp, |x| must be restricted to less than about
sqrt(0.5/((1.5+1.5)/3)) ~= 0.707. We restricted it a bit more to
give a safety margin including some slop for optimizations. Now that
we use double precision for the calculations, the accumulated roundoff
error is in double-precision ulps so it can easily be made almost 2**29
times smaller than a single-precision ulp. Near x = pi/4 its maximum
is about 0.5+(1.5+1.5)*x**2/3 ~= 1.117 double-precision ulps.

The minimax polynomial needs to be different to work for the larger
interval. I didn't increase its degree the old degree is just large
enough to keep the final error less than 1 ulp and increasing the
degree would be a pessimization. The maximum error is now ~0.80
ulps instead of ~0.53 ulps.

The speedup from this optimization for uniformly distributed args in
[-2pi, 2pi] is 28-43% on athlons, depending on how badly gcc selected
and scheduled the instructions in the old version. The old version
has some int-to-float conversions that are apparently difficult to schedule
well, but gcc-3.3 somehow did everything ~10 cycles or ~10% faster than
gcc-3.4, with the difference especially large on AXPs. On A64s, the
problem seems to be related to documented penalties for moving single
precision data to undead xmm registers. With this version, the speed
is cycles is almost independent of the athlon and gcc version despite
the large differences in instruction selection to use the FPU on AXPs
and SSE on A64s.

Use only double precision for "kernel" tanf (except for returning float).
This is a minor interface change. The function is renamed from
__kernel_tanf() to __kernel_tandf() so that misues of it will cause
link errors and not crashes.

This version is a routine translation with no special optimizations
for accuracy or efficiency. It gives an unimportant increase in
accuracy, from ~0.9 ulps to 0.5285 ulps. Almost all of the error is
from the minimax polynomial (~0.03 ulps and the final rounding step
(< 0.5 ulps). It gives strange differences in efficiency in the -5
to +10% range, with -O1 fairly consistently becoming faster and -O2
slower on AXP and A64 with gcc-3.3 and gcc-3.4.

Simplified setiing up args for __kernel_rem_pio2(). We already have x
with a 24-bit fraction, so we don't need a loop to split it into up to
3 terms with 24-bit fractions.

Quick fix for stack buffer overrun in rev.1.13. Oops. The prec == 1
arg to __kernel_rem_pio2() gives 53-bit (double) precision, not single
precision and/or the array dimension like I thought. prec == 2 is
used in e_rem_pio2.c for double precision although it is documented
to be for 64-bit (extended) precision, and I just reduced it by 1
thinking that this would give the value suitable for 24-bit (float)
precision. Reducing it 1 more to the documented value for float
precision doesn't actually work (it gives errors of ~0.75 ulps in the
reduced arg, but errors of much less than 0.5 ulps are needed; the bug
seems to be in kernel_rem_pio2.c). Keep using a value 1 larger than
the documented value but supply an array large enough hold the extra
unused result from this.

The bug can also be fixed quickly by increasing init_jk[0] in
k_rem_pio2.c from 2 to 3. This gives behaviour identical to using
prec == 1 except it doesn't create the extra result. It isn't clear
how the precision bug affects higher precisions. 113-bit (quad) is
the largest precision, so there is no way to use a large precision
to fix it.

Mess up the "kernel" float trig function .c files with ifdefs so that
they can be #included in other .c files to give inline functions, and
use them to inline the functions in most callers (not in e_lgammaf_r.c).
__kernel_tanf() is too large and complicated for gcc to inline very well.

An athlons, this gives a speed increase under favourable pipeline
conditions of about 10% overall (larger for AXP, smaller for A64).
E.g., on AXP, sinf() on uniformly distributed args in [-2Pi, 2Pi]
now takes 30-56 cycles; it used to take 45-61 cycles; hardware fsin
takes 65-129.

Use double precision to simplify and optimize a long division.

On athlons, this gives a speedup of 10-20% for tanf() on uniformly
distributed args in [-2Pi, 2Pi]. (It only directly applies for 43%
of the args and gives a 16-20% speedup for these (more for AXP than
A64) and this gives an overall speedup of 10-12% which is all that it
should; however, it gives an overall speedup of 17-20% with gcc-3.3
on AXP-A64 by mysteriously effected cases where it isn't executed.)

I originally intended to use double precision for all internals of
float trig functions and will probably still do this, but benchmarking
showed that converting to double precision and back is a pessimization
in cases where a simple float precision calculation works, so it may
be optimal to switch precisions only when using extra precision is
much simpler.

Restored a cleanup in rev.1.9 tthat was lost in rev.1.10.

Moved all the optimizations for |x| <= 9pi/2 from
__ieee754_rem_pio2f() to its 3 callers and manually inline them.

On Athlons, with favourable compiler flags and optimizations and
favourable pipeline conditions, this gives a speedup of 30-40 cycles
for cosf(), sinf() and tanf() on the range pi/4 < |x| <= 9pi/4, so
thes functions are now signifcantly faster than the hardware trig
functions in many cases. E.g., in a benchmark with uniformly distributed
x in [-2pi, 2pi], A64 hardware fcos took 72-129 cycles and cosf() took
37-55 cycles. Out-of-order execution is needed to get both of these
times. The optimizations in this commit apparently work more by
removing 1 serialization point than by reducing latency.

Removed an unused declaration which was so old that it wasn't a prototype
and thus just broke building at any nonzero WARNS level.

Fixed nearby style bugs.

-mdoc sweep.

Minor cleanups:

s_cosf.c and s_sinf.c:
Use a non-bogus magic constant for the threshold of pi/4. It was 2 ulps
smaller than pi/4 rounded down, but its value is not critical so it should
be the result of natural rounding.

s_cosf.c and s_tanf.c:
Use a literal 0.0 instead of an unnecessary variable initialized to
[(float)]0.0. Let the function prototype convert to 0.0F.

Improved wording in some comments.

Attempted to improve indentation of comments.

Rearranged the the optimizations for special cases to reduce the average
number of branches.

Use a non-bogus magic constant for the threshold of pi/4. It was 2 ulps
smaller than pi/4 rounded down, but its value is not critical so it should
be the result of natural rounding. Use "<=" comparisons with rounded-
down thresholds for all small multiples of pi/4.

Cleaned up previous commit:
- use static const variables instead of expressions for multiples of pi/2
to ensure that they are evaluated at compile time. gcc currently
evaluates them at compile time but C99 compilers are not required
to do so. We want compile time evaluation for optimization and don't
care about side effects.
- use M_PI_2 instead of a magic constant for pi/2. We need magic constants
related to pi/2 elsewhere but not here since we just want pi/2 rounded
to double and even prefer it to be rounded in the default rounding mode.
We can depend on the cmpiler being C99ish enough to round M_PI_2 correctly
just as much as we depended on it handling hex constants correctly. This
also fixes a harmless rounding error in the hex constant.
- keep using expressions n*<value for pi/2> in the initializers for the
static const variables. 2*M_PI_2 and 4*M_PI_2 are obviously rounded in
the same way as the corresponding infinite precision expressions for
multiples of pi/2, and 3*M_PI_2 happens to be rounded like this, so we
don't need magic constants for the multiples.
- fixed and/or updated some comments.

Fixed some magic numbers.

The threshold for not being tiny was too small. Use the usual 2**-12
threshold. This change is not just an optimization, since the general
code that we fell into has accuracy problems even for tiny x. Avoiding
it fixes 2*1366 args with errors of more than 1 ulp, with a maximum
error of 1.167 ulps.

The magic number 22 is log(DBL_EPSILON)/2 plus slop. This is bogus
for float precision. Use 9 (~log(FLT_EPSILON)/2 plus less slop than
for double precision). The code for handling the interval
[2**-28, 9_was_22] has accuracy problems even for [9, 22], so this
change happens to fix errors of more than 1 ulp in about 2*17000
cases. It leaves such errors in about 2*1074000 cases, with a max
error of 1.242 ulps.

The threshold for switching from returning exp(x)/2 to returning
exp(x/2)^2/2 was a little smaller than necessary. As for coshf(),
This was not quite harmless since the exp(x/2)^2/2 case is inaccurate,
and fixing it avoids accuracy problems in 2*6 cases, leaving problems
in 2*19997 cases.

Fixed naming errors in pseudo-code in comments.

Fixed some magic numbers.

The threshold for not being tiny was confusing and too small. Use the
usual 2**-12 threshold and simplify the algorithm slightly so that
this threshold works (now use the threshold for sinhf() instead of one
for 1+expm1()). This is just a small optimization.

The magic number 22 is log(DBL_EPSILON)/2 plus slop. This is bogus
for float precision. Use 9 (~log(FLT_EPSILON)/2 plus less slop than
for double precision).

The threshold for switching from returning exp(x)/2 to returning
exp(x/2)^2/2 was a little smaller than necessary. This was not quite
harmless since the exp(x/2)^2/2 case is inaccurate. Fixing it happens
to avoid accuracy problems for 2*6 of the 2*151 args that were handled
by the exp(x)/2 case. This leaves accuracy problems for about 2*19997
args near the overflow threshold (~89); the maximum error there is
2.5029 ulps.

There are also accuracy probles for args in +-[0.5*ln2, 9] -- 2*188885
args with errors of more than 1 ulp, with a maximum error of 1.384 ulps.

Fixed a syntax error and naming errors in pseudo-code in comments.

Imoproved comments for the minimax polynomial.

Removed an unused variable.

Fixed some wrong comments and some nearby misformatting.

Tweaked the minimax polynomial and improved its comments.

Improved comments for the minimax polynomial.

As for the float trig functions, use a minimax polynomial that is
specialized for float precision. The new polynomial has degree 8
instead of 14, and a maximum error of 2**-34.34 (absolute) instead of
2**-30.66. This doesn't affect the final error significantly; the
maximum error was and is about 0.8879 ulps on amd64 -01.

The fdlibm expf() is not used on i386's (the "optimized" asm version
is used), but probably should be since it was already significantly
faster than the asm version on athlons. The asm version has the
advantage of being more accurate, so keep using it for now.

As for __kernel_cosf() and __kernel_sinf(), use a fairly optimal minimax
polynomial for __kernel_tanf(). The old one was the double-precision
polynomial with coefficients truncated to float. Truncation is not
a good way to convert minimax polynomials to lower precision. Optimize
for efficiency and use the lowest-degree polynomial that gives a
relative error of less than 1 ulp. It has degree 13 instead of 27,
and happens to be 2.5 times more accurate (in infinite precision) than
the old polynomial (the maximum error is 0.017 ulps instead of 0.041
ulps).

Unlike for cosf and sinf, the old accuracy was close to being inadequate
-- the polynomial for double precision has a max error of 0.014 ulps
and nearly this small an error is needed. The new accuracy is also a
bit small, but exhaustive checking shows that even the old accuracy
was enough. The increased accuracy reduces the maximum relative error
in the final result on amd64 -O1 from 0.9588 ulps to 0.9044 ulps.

Detach k_rem_pio2f.c from the build since it is now unused. It is a libm
internal so this shouldn't cause version problems.

Use a 53-bit approximation to pi/2 instead of a 33+53 bit one for the
special case pi/4 <= |x| < 3*pi/4. This gives a tiny optimization (it
saves 2 subtractions, which are scheduled well so they take a whole 1
cycle extra on an AthlonXP), and simplifies the code so that the
following optimization is not so ugly.

Optimize for the range 3*pi/4 < |x| < 9*Pi/2 in the same way. On
Athlon{XP,64} systems, this gives a 25-40% optimization (depending a
lot on CFLAGS) for the cosf() and sinf() consumers on this range.
Relative to i387 hardware fcos and fsin, it makes the software versions
faster in most cases instead of slower in most cases. The relative
optimization is smaller for tanf() the inefficient part is elsewhere.

The 53-bit approximation to pi/2 is good enough for pi/4 <= |x| <
3*pi/4 because after losing up to 24 bits to subtraction, we still
have 29 bits of precision and only need 25 bits. Even with only 5
extra bits, it is possible to get perfectly rounded results starting
with the reduced x, since if x is nearly a multiple of pi/2 then x is
not near a half-way case and if x is not nearly a multiple of pi/2
then we don't lose many bits. With our intentionally imperfect rounding
we get the same results for cosf(), sinf() and tanf() as without this
optimization.

The logb() functions are not just ieee754 "test" functions, but are
standard in C99 and POSIX.1-2001+. They are also not deprecated, since
apart from being standard they can handle special args slightly better
than the ilogb() functions.

Move their documentation to ilogb.3. Try to use consistent and improved
wording for both sets of functions. All of ieee854, C99 and POSIX
have better wording and more details for special args.

Add history for the logb() functions and ilogbl(). Fix history for
ilogb().

Moved the optimization for tiny x from __kernel_tan[f](x) to tan[f](x)
so that it can be faster for tiny x and avoided for reduced x.

This improves things a little differently than for cosine and sine.
We still need to reclassify x in the "kernel" functions, but we get
an extra optimization for tiny x, and an overall optimization since
tiny reduced x rarely happens. We also get optimizations for space
and style. A large block of poorly duplicated code to fix a special
case is no longer needed. This supersedes the fixes in k_sin.c revs
1.9 and 1.11 and k_sinf.c 1.8 and 1.10.

Fixed wrong constant for the cutoff for "tiny" in tanf(). It was
2**-28, but should be almost the same as the cutoff in sinf() (2**-12).
The incorrect cutoff protected us from the bugs fixed in k_sinf.c 1.8
and 1.10, except 4 cases of reduced args passed the cutoff and needed
special handling in theory although not in practice. Now we essentially
use a cutoff of 0 for the case of reduced args, so we now have 0 special
args instead of 4.

This change makes no difference to the results for sinf() (since it
only changes the algorithm for the 4 special args and the results for
those happen not to change), but it changes lots of results for sin().
Exhaustive testing is impossible for sin(), but exhaustive testing
for sinf() (relative to a version with the old algorithm and a fixed
cutoff) shows that the changes in the error are either reductions or
from 0.5-epsilon ulps to 0.5+epsilon ulps. The new method just uses
some extra terms in approximations so it tends to give more accurate
results, and there are apparently no problems from having extra
accuracy. On amd64 with -O1, on all float args the error range in ulps
is reduced from (0.500, 0.665] to [0.335, 0.500) in 24168 cases and
increased from 0.500-epsilon to 0.500+epsilon in 24 cases. Non-
exhaustive testing by ucbtest shows no differences.

Updated the comment about the optimization for tiny x (the previous
commit moved it). This includes a comment that the "kernel" sine no
longer works on arg -0, so callers must now handle this case. The kernel
sine still works on all other tiny args; without the optimization it is
just a little slower on these args. I intended it to keep working on
all tiny args, but that seems to be impossible without losing efficiency
or accuracy. (sin(x) ~ x * (1 + S1*x**2 + ...) would preserve -0, but
the approximation must be written as x + S1*x**3 + ... for accuracy.)

Removed dead code for handling tan[f]() on odd multiples of pi/2. This
case never occurs since pi/2 is irrational so no multiple of it can
be represented as a float and we have precise arg reduction so we never
end up with a remainder of 0 in the "kernel" function unless the
original arg is 0.

If this case occurs, then we would now fall through to general code
that returns +-Inf (depending on the sign of the reduced arg) instead
of forcing +Inf. The correct handling would be to return NaN since
we would have lost so much precision that the correct result can be
anything _except_ +-Inf.

Don't reindent the else clause left over from this, although it was already
bogusly indented ("if (foo) return; else ..." just marches the indentation
to the right), since it will be removed too.

Index: k_tan.c
===================================================================
RCS file: /home/ncvs/src/lib/msun/src/k_tan.c,v
retrieving revision 1.10
diff -r1.10 k_tan.c
88,90c88
< if (((ix | low) | (iy + 1)) == 0)
< return one / fabs(x);
< else {
---
> {

Fixed some of the silliness related to rev.1.8. In 1.8, "double" in
a declaration was not translated to "float" although bit fiddling on
double variables was translated. This resulted in garbage being put
into the low word of one of the doubles instead of non-garbage being
put into the only word of the intended float. This had no effect on
any result because:
- with doubles, the algorithm for calculating -1/(x+y) is unnecessarily
complicated. Just returning -1/((double)x+y) would work, and the
misdeclaration gave something like that except for messing up some
low bits with the bit fiddling.
- doubles have plenty of bits to spare so messing up some of the low
bits is unlikely to matter.
- due to other bugs, the buggy code is reached for a whole 4 args out
of all 2**32 float args. The bug fixed by 1.8 only affects a small
percentage of cases and a small percentage of 4 is 0. The 4 args
happen to cause no problems without 1.8, so they are even less likely
to be affected by the bug in 1.8 than average args; in fact, neither
1.8 nor this commit makes any difference to the result for these 4
args (and thus for all args).

Corrections to the log message in 1.8: the bug only applies to tan()
and not tanf(), not because the float type can't represent numbers
large enough to trigger the problem (e.g., the example in the fdlibm-5.3
readme which is > 1.0e269), but because:
- the float type can't represent small enough numbers. For there to be
a possible problem, the original arg for tanf() must lie very near an
odd multiple of pi/2. Doubles can get nearer in absolute units. In
ulps there should be little difference, but ...
- ... the cutoff for "small" numbers is bogus in k_tanf.c. It is still
the double value (2**-28). Since this is 32 times smaller than
FLT_EPSILON and large float values are not very uniformly distributed,
only 6 args other than ones that are initially below the cutoff give
a reduced arg that passes the cutoff (the 4 problem cases mentioned
above and 2 non-problem cases).

Fixing the cutoff makes the bug affect tanf() and much easier to detect
than for tan(). With a cutoff of 2**-12 on amd64 with -O1, 670102
args pass the cutoff; of these, there are 337604 cases where there
might be an error of >= 1 ulp and 5826 cases where there is such an
error; the maximum error is 1.5382 ulps.

The fix in 1.8 works with the reduced cutoff in all cases despite the
bug in it. It changes the result in 84492 cases altogether to fix the
5826 broken cases. Fixing the fix by translating "double" to "float"
changes the result in 42 cases relative to 1.8. In 24 cases the
(absolute) error is increased and in 18 cases it is reduced, but it
remains less than 1 ulp in all cases.

Fixed spelling of remquof() in its prototype.

Fixed some comments added in rev.1.5.

The log message for 1.5 said that some small (one or two ulp) inaccuracies
were fixed, and a comment implied that the critical change is to switch
the rounding mode to to-nearest, with a switch of the precision to
extended at no extra cost. Actually, the errors are very large (ucbtest
finds ones of several hundred ulps), and it is the switch of the
precision that is critical.

Another comment was wrong about NaNs being handled sloppily.

Implement inline functions to give the complex result x+I*y from float
or double args x and y. x+I*y cannot be used directly yet due to compiler
bugs.

Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>

Use double precision to simplify and optimize arg reduction for small
and medium size args too: instead of conditionally subtracting a float
17+24, 17+17+24 or 17+17+17+24 bit approximation to pi/2, always
subtract a double 33+53 bit one. The float version is now closer to
the double version than to old versions of itself -- it uses the same
33+53 bit approximation as the simplest cases in the double version,
and where the float version had to switch to the slow general case at
|x| == 2^7*pi/2, it now switches at |x| == 2^19*pi/2 the same as the
double version.

This speeds up arg reduction by a factor of 2 for |x| between 3*pi/4 and
2^7*pi/4, and by a factor of 7 for |x| between 2^7*pi/4 and 2^19*pi/4.

Start trying to make the float precision trig functions actually worth
using under FreeBSD. Before this commit, all float precision functions
except exp2f() were implemented using only float precision, apparently
because Cygnus needed this in 1993 for embedded systems with slow or
inefficient double precision. For FreeBSD, except possibly on systems
that do floating point entirely in software (very old i386 and now
arm), this just gives a more complicated implementation, many bugs,
and usually worse performance for float precision than for double
precision. The bugs and worse performance were particulary large in
arg reduction for trig functions. We want to divide by an approximation
to pi/2 which has as many as 1584 bits, so we should use the widest
type that is efficient and/or easy to use, i.e., double. Use fdlibm's
__kernel_rem_pio2() to do this as Sun apparently intended. Cygnus's
k_rem_pio2f.c is now unused. e_rem_pio2f.c still needs to be separate
from e_rem_pio2.c so that it can be optimized for float args. Similarly
for long double precision.

This speeds up cosf(x) on large args by a factor of about 2. Correct
arg reduction on large args is still inherently very slow, so hopefully
these args rarely occur in practice. There is much more efficiency
to be gained by using double precision to speed up arg reduction on
medium and small float args.

Use fairly optimal minimax polynomials for __kernel_cosf() and
__kernel_sinf(). The old ones were the double-precision polynomials
with coefficients truncated to float. Truncation is not a good way
to convert minimax polynomials to lower precision. Optimize for
efficiency and use the lowest-degree polynomials that give a relative
error of less than 1 ulp -- degree 8 instead of 14 for cosf and degree
9 instead of 13 for sinf. For sinf, the degree 8 polynomial happens
to be 6 times more accurate than the old degree 14 one, but this only
gives a tiny amount of extra accuracy in results -- we just need to
use a a degree high enough to give a polynomial whose relative accuracy
in infinite precision (but with float coefficients) is a small fraction
of a float ulp (fdlibm generally uses 1/32 for the small fraction, and
the fraction for our degree 8 polynomial is about 1/600).

The maximum relative errors for cosf() and sinf() are now 0.7719 ulps
and 0.7969 ulps, respectively.

Use a better algorithm for reducing the error in __kernel_cos[f]().
This supersedes the fix for the old algorithm in rev.1.8 of k_cosf.c.

I want this change mainly because it is an optimization. It helps
make software cos[f](x) and sin[f](x) faster than the i387 hardware
versions for small x. It is also a simplification, and reduces the
maximum relative error for cosf() and sinf() on machines like amd64
from about 0.87 ulps to about 0.80 ulps. It was validated for cosf()
and sinf() by exhaustive testing. Exhaustive testing is not possible
for cos() and sin(), but ucbtest reports a similar reduction for the
worst case found by non-exhaustive testing. ucbtest's non-exhaustive
testing seems to be good enough to find problems in algorithms but not
maximum relative errors when there are spikes. E.g., short runs of
it find only 3 ulp error where the i387 hardware cos() has an error
of about 2**40 ulps near pi/2.

More fixes for arg reduction near pi/2 on systems with broken assignment
to floats (mainly i386's). All errors of more than 1 ulp for float
precision trig functions were supposed to have been fixed; however,
compiling with gcc -O2 uncovered 18250 more such errors for cosf(),
with a maximum error of 1.409 ulps.

Use essentially the same fix as in rev.1.8 of k_rem_pio2f.c (access a
non-volatile variable as a volatile). Here the -O1 case apparently
worked because the variable is in a 2-element array and it takes -O2
to mess up such a variable by putting it in a register.

The maximum error for cosf() on i386 with gcc -O2 is now 0.5467 (it
is still 0.5650 with gcc -O1). This shows that -O2 still causes some
extra precision, but the extra precision is now good.

Extra precision is harmful mainly for implementing extra precision in
software. We want to represent x+y as w+r where both "+" operations
are in infinite precision and r is tiny compared with w. There is a
standard algorithm for this (Knuth (1981) 4.2.2 Theorem C), and fdlibm
uses this routinely, but the algorithm requires w and r to have the
same precision as x and y. w is just x+y (calculated in the same
finite precision as x and y), and r is a tiny correction term. The
i386 gcc bugs tend to give extra precision in w, and then using this
extra precision in the calculation of r results in the correction
mostly staying in w and being missing from r. There still tends to
be no problem if the result is a simple expression involving w and r
-- modulo spills, w keeps its extra precision and r remains the right
correction for this wrong w. However, here we want to pass w and r
to extern functions. Extra precision is not retained in function args,
so w gets fixed up, but the change to the tiny r is tinier, so r almost
remains as a wrong correction for the right w.

Moved the optimization for tiny x from __kernel_{cos,sin}[f](x) to
{cos_sin}[f](x) so that x doesn't need to be reclassified in the
"kernel" functions to determine if it is tiny (it still needs to be
reclassified in the cosine case for other reasons that will go away).

This optimization is quite large for exponentially distributed x, since
x is tiny for almost half of the domain, but it is a pessimization for
uniformally distributed x since it takes a little time for all cases
but rarely applies. Arg reduction on exponentially distributed x
rarely gives a tiny x unless the reduction is null, so it is best to
only do the optimization if the initial x is tiny, which is what this
commit arranges. The imediate result is an average optimization of
1.4% relative to the previous version in a case that doesn't favour
the optimization (double cos(x) on all float x) and a large
pessimization for the relatively unimportant cases of lgamma[f][_r](x)
on tiny, negative, exponentially distributed x. The optimization should
be recovered for lgamma*() as part of fixing lgamma*()'s low-quality
arg reduction.

Fixed various wrong constants for the cutoff for "tiny". For cosine,
the cutoff is when x**2/2! == {FLT or DBL}_EPSILON/2. We round down
to an integral power of 2 (and for cos() reduce the power by another
1) because the exact cutoff doesn't matter and would take more work
to determine. For sine, the exact cutoff is larger due to the ration
of terms being x**2/3! instead of x**2/2!, but we use the same cutoff
as for cosine. We now use a cutoff of 2**-27 for double precision and
2**-12 for single precision. 2**-27 was used in all cases but was
misspelled 2**27 in comments. Wrong and sloppy cutoffs just cause
missed optimizations (provided the rounding mode is to nearest --
other modes just aren't supported).

Fixed range reduction for large multiples of pi/2 on systems with
broken assignment to floats (e.g., i386 with gcc -O, but not amd64 or
ia64; i386 with gcc -O0 worked accidentally).

Use an unnamed volatile temporary variable to trick gcc -O into clipping
extra precision on assignment. It's surprising that only 1 place needed
to be changed.

For tanf() on i386 with gcc -O, the bug caused errors > 1 ulp with a
density of 2.3% for args larger in magnitude than 128*pi/2, with a
maximum error of 1.624 ulps.

After this fix, exhaustive testing shows that range reduction for
floats works as intended assuming that it is in within a factor of
about 2^16 of working as intended for doubles. It provides >= 8
extra bits of precision for all ranges. On i386:

range max error in double/single ulps extra precision
----- ------------------------------- ---------------
0 to 3*pi/4 0x000d3132 / 0.0016 9+ bits
3*pi/4 to 128*pi/2 0x00160445 / 0.0027 8+
128*pi/2 to +Inf 0x00000030 / 0.00000009 23+
128*pi/2 up, -O0 before fix 0x00000030 / 0.00000009 23+
128*pi/2 up, -O1 before fix 0x10000000 / 0.5 1

The 23+ bits of extra precision for large multiples corresponds to almost
perfect reduction to a pair of floats (24 extra would be perfect).

After this fix, the maximum relative error (relative to the corresponding
fdlibm double precision function) is < 1 ulp for all basic trig functions
on all 2^32 float args on all machines tested:

amd64 ia64 i386-O0 i386-O1
------ ------ ------ ------
cosf: 0.8681 0.8681 0.7927 0.5650
sinf: 0.8733 0.8610 0.7849 0.5651
tanf: 0.9708 0.9329 0.9329 0.7035

Fixed range reduction near (but not very near) medium-sized multiples
of pi/2 (1 line) and expand a comment about related magic (many lines).

The bug was essentially the same as for the +-pi/2 case (a mistranslated
mask), but was smaller so it only significantly affected multiples
starting near +-13*pi/2. At least on amd64, for cosf() on all 2^32
float args, the bug caused 128 errors of >= 1 ulp, with a maximum error
of 1.2393 ulps.

Fix numerous errors of >= 1 ulp for cosf(x) and sinf(x) (1 line)
and add a comment about related magic (many lines)).

__kernel_cos[f]() needs a trick to reduce the error to below 1 ulp
when |x| >= 0.3 for the range-reduced x. Modulo other bugs, naive
code that doesn't use the trick would have an error of >= 1 ulp
in about 0.00006% of cases when |x| >= 0.3 for the unreduced x,
with a maximum relative error of about 1.03 ulps. Mistransation
of the trick from the double precision case resulted in errors in
about 0.2% of cases, with a maximum relative error of about 1.3 ulps.

The mistranslation involved not doing implicit masking of the 32-bit
float word corresponding to to implicit masking of the lower 32-bit
double word by clearing it.

sinf() uses __kernel_cosf() for half of all cases so its errors from
this bug are similar. tanf() is not affected.

The error bounds in the above and in my other recent commit messages
are for amd64. Extra precision for floats on i386's accidentally masks
this bug, but only if k_cosf.c is compiled with -O. Although the extra
precision helps here, this is accidental and depends on longstanding
gcc precision bugs (not clipping extra precision on assignment...),
and the gcc bugs are mostly avoided by compiling without -O. I now
develop libm mainly on amd64 systems to simplify error detection and
debugging.

Oops, the last-minute optimization in rev.1.8 wasn't a good idea. The
17+17+24 bit pi/2 must only be used when subtraction of the first 2
terms in it from the arg is exact. This happens iff the the arg in
bits is one of the 2**17[-1] values on each side of (float)(pi/2).

Revert to the algorithm in rev.1.7 and only fix its threshold for using
the 3-term pi/2. Use the threshold that maximizes the number of values
for which the 3-term pi/2 is used, subject to not changing the algorithm
for comparing with the threshold. The 3-term pi/2 ends up being used
for about half of its usable range (about 64K values on each side).

Fixed syntax error (a missing brace) in previous commit.

Fixed range reduction near (but not very near) +-pi/2. A bug caused
a maximum error of 2.905 ulps for cosf(), but the algorithm for cosf()
is good for < 1 ulps and happens to give perfect rounding (< 0.5 ulps)
near +-pi/2 except for the bug. The extra relative errors for tanf()
were similar (slightly larger). The bug didn't affect sinf() since
sinf'(+-pi/2) is 0.

For range reduction in ~[-3pi/4, -pi/4] and ~[pi/4, 3pi/4] we must
subtract +-pi/2 and the only complication is that this must be done
in extra precision. We have handy 17+24-bit and 17+17+24-bit
approximations to pi/2. If we always used the former then we would
lose up to 24 bits of accuracy due to cancelation of leading bits, but
we need to keep at least 24 bits plus a guard digit or 2, and should
keep as many guard bits as efficiency permits. So we used the
less-precise pi/2 not very near +-pi/2 and switched to using the
more-precise pi/2 very near +-pi/2. However, we got the threshold for
the switch wrong by allowing 19 bits to cancel, so we ended up with
only 21 or 22 bits of accuracy in some cases, which is even worse than
naively subtracting pi/2 would have done.

Exhaustive checking shows that allowing only 17 bits to cancel (min.
accuracy ~24 bits) is sufficient to reduce the maximum error for cosf()
near +-pi/2 to 0.726 ulps, but allowing only 6 bits to cancel (min.
accuracy ~35-bits) happens to give perfect rounding for cosf() at
little extra cost so we prefer that.

We actually (in effect) allow 0 bits to cancel and always use the
17+17+24-bit pi/2 (min. accuracy ~41 bits). This is simpler and
probably always more efficient too. Classifying args to avoid using
this pi/2 when it is not needed takes several extra integer operations
and a branch, but just using it takes only 1 FP operation.

The patch also fixes misspelling of 17 as 24 in many comments.

For the double-precision version, the magic numbers include 33+53 bits
for the less-precise pi/2 and (53-32-1 = 20) bits being allowed to
cancel, so there are ~33-20 = 13 guard bits. This is sufficient except
probably for perfect rounding. The more-precise pi/2 has 33+33+53
bits and we still waste time classifying args to avoid using it.

The bug is apparently from mistranslation of the magic 32 in 53-32-1.
The number of bits allowed to cancel is not critical and we use 32 for
double precision because it allows efficient classification using a
32-bit comparison. For float precision, we must use an explicit mask,
and there are fewer bits so there is less margin for error in their
allocation. The 32 got reduced to 4 but should have been reduced
almost in proportion to the reduction of mantissa bits.

Fixed aliasing bugs in TRUNC() by using the fdlibm macros for access
to doubles as bits. fdlibm-1.1 had similar aliasing bugs, but these
were fixed by NetBSD or Cygnus before a modified version of fdlibm was
imported in 1994. TRUNC() is only used by tgamma() and some
implementation-detail functions. The aliasing bugs were detected by
compiling with gcc -O2 but don't seem to have broken tgamma() on i386's
or amd64's. They broke my modified version of tgamma().

Moved the definition of TRUNC() to mathimpl.h so that it can be fixed
in one place, although the general version is even slower than necessary
because it has to operate on pointers to volatiles to handle its arg
sometimes being volatile. Inefficiency of the fdlibm macros slows
down libm generally, and tgamma() is a relatively unimportant part of
libm. The macros act as if on 32-bit words in memory, so they are
hard to optimize to direct actions on 64-bit double registers for
(non-i386) machines where this is possible. The optimization is too
hard for gcc on amd64's, and declaring variables as volatile makes it
impossible.

Add a missing ldexpf() alias for amd64.

Noticed by: bz@, tjr@

Bump the shared library version number of all libraries that have not
been bumped since RELENG_5.

Reviewed by: ru
Approved by: re (not needed for commit check but in principle...)

Markup nit.

Approved by: re (blanket)

Fixed compile warning.

Approved by: re (blanket)

Assorted markup fixes.

Approved by: re

Prevent these functions from using stack outside of their frame.

Reported by: Marc Olzheim <marcolz at stack dot nl>
OK'd by: das

Revert the last change, the conversion from long double to double can raise
unwanted underflow exceptions.

Pointed out by: das

Use double additions to raise the inexact exception to work around problems
with long double addition on sparc64.

Fix raising the inexact exception (FE_INEXACT) if the result differs from the
argument.

Noticed by: das

Fix truncl.3 MLINKS

More optimized math functions.

Implement truncl() based on floorl().

Add roundl(), lroundl(), and llroundl().

These files should include s_lround.c instead of s_lrint.c.
This only matters for efficiency, not for correctness.

Fix a (coincidentally harmless) bug.

Fix a long-standing bug in k_rem_pio2(), which led to large errors when
tanf() was called with big arguments close to multiples of pi/2.

Reported by: ucbtest via bde

Build exp2(), exp2f(), and related documentation.

Document exp2() and exp2f(), and make other minor tweaks and updates.

Implement exp2() and exp2f().

Implement and document remquo() and remquof().

Fix the double rounding problem with subnormals, and
remove the XXX comments, which no longer apply.

Add missing prototypes for fma() and fmaf(), and remove an inaccurate
comment.

Make the fenv.h routines work for programs that use SSE for
floating-point arithmetic on i386. Now I'm going to make excuses
for why this code is kinda scary:

- To avoid breaking the ABI with 5.3-RELEASE, we can't change
sizeof(fenv_t). I stuck the saved mxcsr in some discontiguous
reserved bits in the existing structure.

- Attempting to access the mxcsr on older processors results
in an illegal instruction exception, so support for SSE must
be detected at runtime. (The extra baggage is optimized away
if either the application or libm is compiled with -msse{,2}.)

I didn't run tests to ensure that this doesn't SIGILL on older 486's
lacking the cpuid instruction or on other processors lacking SSE.
Results from running the fenv regression test on these processors
would be appreciated. (You'll need to compile the test with
-DNO_STRICT_DFL_ENV.) If you have an 80386, or if your processor
supports SSE but the kernel didn't enable it, then you're probably out
of luck.

Also, I un-inlined some of the functions that grew larger as a result
of this change, moving them from fenv.h to fenv.c.

Spell 'fedisableexcept' correctly.

Document feenableexcept(), fedisableexcept(), and fegetexcept().

Replace fegetmask() and fesetmask() with feenableexcept(),
fedisableexcept(), and fegetexcept(). These two sets of routines
provide the same functionality. I implemented the former as an
undocumented internal interface to make the regression test easier to
write. However, fe(enable|disable|get)except() is already part of
glibc, and I would like to avoid gratuitous differences. The only
major flaw in the glibc API is that there's no good way to report
errors on processors that don't support all the unmasked exceptions.

Replace strong references with weak references. There's no
particularly good reason to do this, except that __strong_reference
does type checking, whereas __weak_reference does not.
On Alpha, the compiler won't accept a 'long double' parameter in
place of a 'double' parameter even thought the two types are
identical.

Remove an obsolete sentence from a comment.

- If z is 0, one of x or y is 0, and the other is infinite, raise
an invalid exception and return an NaN.
- If a long double has 113 bits of precision, implement fma in terms
of simple long double arithmetic instead of complicated double arithmetic.
- If a long double is the same as a double, alias fma as fmal.

Document scalbnl and scalblnl.

Document nextafterl and nexttoward{,f,l}.

Add nexttoward to the list of implemented functions, and explicitly
list the four that are still missing.

Document fmal.

Remove ldexp and ldexpf. The former is in libc, and the latter is
identical to scalbnf, which is now aliased as ldexpf. Note that the
old implementations made the mistake of setting errno and were the
only libm routines to do so.

- Remove s_ldexpf.c (now aliased to scalbn.)
- Add nexttoward{,f,l} and nextafterl. On all platforms,
nexttowardl is an alias for nextafterl.
- Add fmal.
- Add man pages for new routines: fmal, nextafterl,
nexttoward{,f,l}, scalb{,l}nl.

Note that on platforms where long double is the same as double, we
generally just alias the double versions of the routines, since doing
so avoids extra work on the source code level and redundant code in
the binary. In particular:

ldbl53 ldbl64/113
fmal s_fma.c s_fmal.c
ldexpl s_scalbn.c s_scalbnl.c
nextafterl s_nextafter.c s_nextafterl.c
nexttoward s_nextafter.c s_nexttoward.c
nexttowardf s_nexttowardf.c s_nexttowardf.c
nexttowardl s_nextafter.c s_nextafterl.c
scalbnl s_scalbn.c s_scalbnl.c

- Define FP_FAST_FMA for sparc64, since fma() is now implemented using
sparc64's 128-bit long doubles.
- Define FP_FAST_FMAL for ia64.
- Prototypes for fmal, frexpl, ldexpl, nextafterl, nexttoward{,f,l},
scalblnl, and scalbnl.

Alias scalbn as ldexpl and scalbnl on platforms where long double is
the same as double.

- Implement scalblnl.
- In scalbln and scalblnf, check the bounds of the second argument.
This is probably unnecessary, but strictly speaking, we should
report an error if someone tries to compute scalbln(x, INT_MAX + 1ll).

Implement nexttowardf. This is used on both platforms with 11-bit
exponents and platforms with 15-bit exponents for long doubles.

Implement nexttoward and nextafterl; the latter is also known as
nexttowardl. These are not needed on machines where long doubles
look like IEEE-754 doubles, so the implementation only supports
the usual long double formats with 15-bit exponents.

Anything bizarre, such as machines where floating-point and integer
data have different endianness, will cause problems. This is the case
with big endian ia64 according to libc/ia64/_fpmath.h. Please contact
me if you managed to get a machine running this way.

- Try harder to trick gcc into not optimizing away statements
that are intended to raise underflow and inexact exceptions.
- On systems where long double is the same as double, nextafter
should be aliased as nexttoward, nexttowardl, and nextafterl.

Implement frexpl.

Alias frexp as frexpl on platforms where a long double is the same as
a double.

Implement fmal.

- Define the LDBL_PREC to be the number of significant bits in a long
double's mantissa.
- Add an assembly version of fmal.

- Define the LDBL_PREC to be the number of significant bits in a long
double's mantissa.
- Add an assembly version of scalbnl.

Define the LDBL_PREC to be the number of significant bits in a long
double's mantissa.

Add an assembly version of fmal.

Add scalbnl, also known as as ldexpl.

Alias scalbnf as ldexpf. The two are identical in binary
floating-point formats.

Fix a mistake in the exponent range.

Work around a gcc bug. This fixes feholdexcept() et al. at -O1.
Symptoms of the problem included assembler warnings and
nondeterministic runtime behavior when a fe*() call that affects the
fpsr is closely followed by a float point op.

The bug (at least, I think it's a bug) is that gcc does not insert a
break between a volatile asm and a dependent instruction if the
volatile asm came from an inlined function. Volatile asms seem to be
fine in other circumstances, even without -mvolatile-asm-stop, so
perhaps the compiler adds the stop bits before inlining takes place.
The problem does not occur at -O0 because inlining is disabled, and it
doesn't happen at -O2 because -fschedule-insns2 knows better.

Un-document the non-extant exp10() and exp10f() functions.
exp10() was a casualty of the transition away from the VAX.

Revert rev 1.8, which causes small (e.g. 2 ulp) errors for some
inputs. The trouble with replacing two floats with a double is that
the latter has 6 extra bits of precision, which actually hurts
accuracy in many cases. All of the constants are optimal when float
arithmetic is used, and would need to be recomputed to do this right.

Noticed by: bde (ucbtest)

Use hardware instructions for sqrt() and sqrtf().

Use double arithmetic instead of simulating it with two floats. This
results in a performance gain on the order of 10% for amd64 (sledge),
ia64 (pluto1), i386+SSE (Pentium 4), and sparc64 (panther), and a
negligible improvement for i386 without SSE. (The i386 port still
uses the hardware instruction, though.)

Remove the i387 versions of atan(), atan2(), and atan2f().
They are slower than the MI routines on modern hardware,
except for degenerate cases such as the Pentium 4.

PR: 67469

Remove i387 versions of asin() and acos(). Although the hardware
instruction was faster on the 486, it's slower than our MD version on
modern processors.

Determined by: bde
PR: 67469

Remove the float versions of the i387 trig functions obtained from
NetBSD. They're buggy, giving particularly for inputs larger in
magnitude than 2**63.

Noticed by: bde
PR: 67469

Fix a small scripting snafu in the previous revision.

Remove another vestige of support for a non-IEEE libm.

Reduce diffs against vendor source (Sun fdlibm 5.3).

Move machine-dependent crud to its own makefile.

Remove wrappers and other cruft intended to support SVID, mistakes in
C90, and other arcana. Most of these features were never fully
supported or enabled by default.

Ok: bde, stefanf

Typo.

Properly terminate sentence.

- Move the functions presently described in in ieee(3) to their own
manpages. They are not very related, so separating them makes it
easier to add meaningful cross-references and extend some of the
descriptions.
- Move the part of math(3) that discusses IEEE 754 to the ieee(3)
manpage.

Define FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD and _ROUND_MASK to
unbreak the build for arm.

Update comment to reflect the code change in the previous revision.

Noticed by: ceri

Many changes, including the following major ones:
- Rearrange the list of functions into categories.
- Remove the ulps column. It was appropriate for only some
of the functions in the list, and correct for even fewer
of them.
- Add some new paragraphs, and remove some old ones about
NaNs that may do more harm than good.
- Document precisions other than double-precision.

If x == y, return y, not x. C99 (though not IEEE 754) requires that
nextafter(+0.0, -0.0) returns -0.0 and nextafter(-0.0, +0.0) returns +0.0.

Add fma() and fmaf(), which implement a fused multiply-add operation.

Sort sections.

Use the \*(If string provided by mdoc(7), to represent infinity.

Removed redundant .br call.

amd64 assembly versions of sqrt(), lrint(), and llrint() using SSE2.

Most libm routines depend on the rounding mode and/or set exception
flags, so they are not pure. Remove the __pure2 annotation from them.
I believe that the following routines and their float and long double
counterparts are the only ones here that can be __pure2:

copysign is* fabs finite fmax fmin fpclassify ilogb nan signbit

When gcc supports FENV_ACCESS, perhaps there will be a new annotation
that allows the other functions to be considered pure when FENV_ACCESS
is off.

Discussed with: bde

Braino. Revert rev 1.50.

Pointy hat to: das

Remove numerous references to VAX floating-point and the setting of
errno, replacing them with a discussion of IEEE exceptions where
appropriate. Cross-reference fenv(3) whenever exceptions are
mentioned.

Set math_errhandling to MATH_ERREXCEPT. Now that we have fenv.h, we
basically support this, subject to gcc's lack of FENV_ACCESS support.
In any case, the previous setting of math_errhandling to 0 is not
allowed by POSIX.

Remove some #if 0'd code.

Tiny markup nits.

Mark all inline asms that read the floating-point control or status
registers as volatile. Instructions that *wrote* to FP state were
already marked volatile, but apparently gcc has license to move
non-volatile asms past volatile asms. This broke amd64's feupdateenv
at -O2 due to a WAR conflict between fnstsw and fldenv there.

Fixed too many of "the", and enclose multi-word argument in double quotes.

Obtained from: ru

Import the subset of J.T. Conklin's single-precision x86-optimized
math routines that appear to be (a) correct and (b) faster than their
MI counterparts on my Pentium 4.

Obtained from: NetBSD

The isnormal() in rev 1.2 should have been isfinite() so subnormals
round correctly.

Noticed by: stefanf

Things that are broken, unneeded, and unused since 1997 belong in the attic.

Markup nits.

Fixed too many of "the", and enclose multi-word argument in double quotes.

Implement and document ceill().

Bump .Dd for the last commit.

Hook up and document floorl().

Implement floorl().

Whitespace nit.

Add MI implementations of [l]lrint[f]() and [l]lround[f]().

Discussed with: bde

Document [l]lrint[f]() and [l]lround[f]().

Faster lrint() and llrint() implementations for x86.

Mark inline stmxcsr instructions as volatile, since this appears to be
the only way to convince gcc that they read the MXCSR. The volatile
annotation may be needed elsewhere as well.

Scheduled mdoc(7) sweep.

Sanitize the markup, as prompted.

GC unused declaration

Cosmetic changes only:
- style
- remove unused variables
- de-support VAX

Inspired by: bin/42388

More updates for math(3):
- Make some minor rearrangements in the introduction.
- Mention the problem with argument reduction on i386.
- Add recently-implemented functions to the table.
- Un-document the error bounds that only apply to the old 4BSD math
library, and fill in the correct values where I know them. No
attempt has been made to document bounds lower than 1 ulp, although
smaller bounds are usually achievable in round-to-nearest mode.

Add and document ilogbl(), a long double version of ilogb().

Use the FP_ILOG macros from <math.h> rather than hardcoded return values.
Also be prepared for FP_ILOGBNAN != INT_MAX.

Reviewed by: md5

Bump the library version numbers for the following libraries:

/lib/{libm,libreadline}
/usr/lib/{libhistory,libopie,libpcap}

in preparation for doing the same thing to RELENG_5. HUGE amounts of
help for determining what to bump provided by kris.

Discussed on: freebsd-current
Approved by: re (not required for commit but something like this should be)

Further refine some #ifs:
- Simplify the logic by using __GNUC_PREREQ__.
Suggested by stefanf.
- Make math.h compile with old (pre-8.0) versions of icc.
Submitted by sf [sic].

Add man pages for the cimag(), conj() and creal() functions.

Only use rfs and wfs if ARM_HARD_FLOAT is defined, and use stubs if it is not,
in order to unbreak arm make world. The right way to do it with soft floats
will be figured out later.
Discussed with: das

Replace s_isnan.c and s_isnanf.c with the more compact s_isnan.c from
libc. The externally-visible effect of this is to add __isnanl() to
libm, which means that libm.so.2 can once again link against libc.so.4
when LD_BIND_NOW is set. This was broken by the addition of fdiml(),
which calls __isnanl().