Commit Graph

840 Commits

Author SHA1 Message Date
Warner Losh
98fd12e162 Revert r241756 2012-10-22 13:21:11 +00:00
Warner Losh
b3f94d242f Revert r241755 2012-10-22 13:20:31 +00:00
Warner Losh
527a7c56f7 Document the method used to compute expf. Taken from exp, with
changes to reflect differences in computation between the two.
2012-10-19 22:47:44 +00:00
Warner Losh
aa616aa22d Document the methods used to compute logf. Taken and edited from the
double version, with adaptations for the differences between it and
the float version.
2012-10-19 22:46:48 +00:00
Joel Dahl
284eab32c0 mdoc: avoid nested displays. 2012-10-14 14:45:54 +00:00
Steve Kargl
f81d134e7e * Update the comment that explains the choice of values in the
table and the requirement on trailing zero bits.

* Remove the __aligned() compiler directives as these were found
  to have a negative effect on the produced code.

Submitted by:	bde
Approved by:	das (mentor)
2012-10-13 19:53:11 +00:00
Steve Kargl
a077586c53 * src/math_private.h:
. Change the API for the LD80C by removing the explicit passing
    of the sign bit.  The sign can be determined from the last
    parameter of the macro.
  . On i386, load long double by bit manipulations to work around
    at least a gcc compiler issue.  On non-i386 ld80 architectures,
    use a simple assignment.

* ld80/s_expl.c:
  . Update the only consumer of LD80C.

Submitted by:	bde
Approved by:	das (mentor)
2012-09-29 16:40:12 +00:00
Steve Kargl
dba466c344 * ld80/s_expl.c:
. Fix the threshold for expl(x) where |x| is small.
  . Also update the previously incorrect comment to match the
    new threshold.

* ld128/s_expl.c:
  . Re-order logic in exceptional cases to match the logic used in
    other long double functions.
  . Fix the threshold for expl(x) where is |x| is small.
  . Also update the previously incorrect comment to match the
    new threshold.

Submitted by:	bde
Approved by:	das (mentor)
2012-09-23 18:32:03 +00:00
Steve Kargl
724c1ee29f Fix whitespace issue.
Approved by:	das (mentor, implicit)
2012-09-23 18:13:46 +00:00
Steve Kargl
8f647ffd7f * ld80/s_expl.c:
. Guard a comment from reformatting by indent(1).
  . Re-order variables in declarations to alphabetical order.
  . Remove a banal comment.

* ld128/s_expl.c:
  . Add a comment to point to ld80/s_expl.c for implementation details.
  . Move the #define of INTERVAL to reduce the diff with ld80/s_expl.c.
  . twom10000 does not need to be volatile, so move its declaration.
  . Re-order variables in declarations to alphabetical order.
  . Add a comment that describes the argument reduction.
  . Remove the same banal comment found in ld80/s_expl.c.

Reviewed by:	bde
Approved by:	das (mentor)
2012-09-23 18:06:27 +00:00
Steve Kargl
c1a077829a * Update the lookup table to use 53-bit high and low values.
Also, update the comment to describe the choice of using
  a high and low decomposition of 2^(i/INTERNVAL) for
  0 <= i <= INTERVAL in preparation for an implementation of
  expm1l.

* Move the #define of INTERVAL above the comment, because the
  comment refers to INTERVAL.

Reviewed by:	bde
Approved by:	das (mentor)
2012-09-23 17:36:01 +00:00
Steve Kargl
340076f0f6 * Use ENTERI() and RETURNI() to toggle the rounding precision if
necessary, so that cosl(), sinl() and tanl() work on i386 even
  for naive callers.

Suggested by:	bde
Reviewed by:	bde
Approved by: 	das (mentor)
2012-09-22 15:38:29 +00:00
Steve Kargl
e6f9129aa0 * Make STRICT_ALIGN() work for doubles as well as for floats. This
only affects i386.  The double case was intentionally left broken
  as an optimization, but we are getting closer to supporting
  applications and/or kernels that change the (FreeBSD i386) default
  rounding precision from FP_PD to FP_PE and never change it back,
  and this requires the STRICT_ALIGN()s that were added to support
  FP_PE to actually work in all precisions.

* Remove an extraneous semicolon at the end of a macro that was
  supposed to be function-like.

Submitted by:	bde
Approved by:	das (mentor)
2012-09-22 15:19:11 +00:00
Tijl Coosemans
71dad5d6ad Optimise i387 trigonometric functions. Replace "andw 0x400,%ax \ jnz" with
"sahf \ jp", "fprem1" with "fprem" and "fstsw %ax" with "fnstsw %ax".
2012-09-16 16:58:49 +00:00
Joel Dahl
ff34c32ac4 Minor mdoc fix. 2012-09-11 17:40:06 +00:00
Marcel Moolenaar
7750ad47a9 Sync FreeBSD's bmake branch with Juniper's internal bmake branch.
Requested by: Simon Gerraty <sjg@juniper.net>
2012-08-22 19:25:57 +00:00
Dimitry Andric
f049a6cb4a 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
MFC after:	1 week
2012-08-11 15:47:22 +00:00
Dimitry Andric
2b795b2921 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 after:	1 week
2012-08-11 11:13:48 +00:00
Steve Kargl
82b89f4850 ieeefp.h is only needed on i386 class hardware.
Submitted by:	bde
Approved by:	das (pre-approved)
2012-07-30 21:58:28 +00:00
Steve Kargl
ca50c4b871 Whitespace.
Submitted by:	bde
Approved by:	das (pre-approved)
2012-07-30 21:55:49 +00:00
Steve Kargl
8345cbd275 Replace the macro name NUM with INTERVALS. This change provides
compatibility with the INTERVALS macro used in the soon-to-be-commmitted
expm1l() and someday-to-be-committed log*l() functions.

Add a comment into ld128/s_expl.c noting at gcc issue that was
deleted when rewriting ld80/e_expl.c as ld128/s_expl.c.

Requested by:	bde
Approved by:	das (mentor)
2012-07-26 04:05:08 +00:00
Steve Kargl
f7cfe68f59 * ld80/expl.c:
. Remove a few #ifdefs that should have been removed in the initial
    commit.
  . Sort fpmath.h to its rightful place.

* ld128/s_expl.c:
  . Replace EXPMASK with its actual value.
  . Sort fpmath.h to its rightful place.

Requested by:	bde
Approved by:	das (mentor)
2012-07-26 03:59:33 +00:00
Steve Kargl
3458a81e65 Replace code that toggles between 53 and 64 bits on i386
class hardware with the ENTERI and RETURNI macros, which
are now available in math_private.h.

Suggested by:	bde
Approved by: das (mentor)
2012-07-26 03:50:24 +00:00
Steve Kargl
8a489dc248 Hook ld80/s_expl.c or ld128/s_expl.c into the building of libm.
PR: standards/152415
Approved by: das (mentor)
2012-07-23 19:23:49 +00:00
Steve Kargl
b83ccea32c Compute the exponential of x for Intel 80-bit format and IEEE 128-bit
format.  These implementations are based on

PTP Tang, "Table-driven implementation of the exponential function
in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15,
144-157 (1989).

PR: standards/152415
Submitted by: kargl
Reviewed by: bde, das
Approved by: das (mentor)
2012-07-23 19:13:55 +00:00
David Chisnall
8f62d384b9 Allow inclusion of libc++ <cmath> to work after including math.h
Submitted by:	Yamaya Takashi
Reviewed by:	das
MFC after:	1 week
2012-05-27 12:54:41 +00:00
Glen Barber
a3fb6da9ba General mdoc(7) and typo fixes.
PR:		167734
Submitted by:	Nobuyuki Koganemaru (kogane!jp.freebsd.org)
MFC after:	3 days
2012-05-11 20:06:46 +00:00
Dag-Erling Smørgrav
084f679591 I stopped using my middle name years ago. 2012-04-25 18:07:35 +00:00
Joel Dahl
748611c9c0 mdoc: fix column names, indentation, column separation within each row, and
quotation. Also make sure we have the same amount of columns in each row as
the number of columns we specify in the head arguments.

Reviewed by:	brueffer
2012-04-07 09:05:30 +00:00
David Schultz
1cbd288942 Fix a bug in remquo{,f,l}, in which the quotient didn't always have the
correct sign when the remainder was 0.

Fix a separate bug in remquo alone, in which the remainder and
quotient were both off by a bit in certain cases involving subnormal
remainders.

The bugs affected all platforms except amd64 and i386, on which the
routines are implemented in assembly.

PR:		166463
Submitted by:	Ilya Burylov
MFC after:	2 weeks
2012-04-07 03:59:12 +00:00
Joel Dahl
41949a1ed5 Remove superfluous paragraph macro. 2012-03-25 12:13:24 +00:00
David Schultz
6e3ab4e1b3 Fix a small nit noted by bde: exp_x should be of type float, not double. 2012-01-20 07:02:42 +00:00
David Schultz
b7e83566eb Add a change I missed in r230367 (don't inline arm's fenv.h functions). 2012-01-20 07:01:58 +00:00
David Schultz
7082a2cf72 Don't inline fenv.h functions on arm for now. Inlining makes sense:
the function bodies require only 2 to 10 instructions.  However, it
leads to application binaries that refer to a private ABI, namely, the
softfloat innards in libc.  This could complicate future changes in
the implementation of the floating-point emulation layer, so it seems
best to have programs refer to the official fe* entry points in libm.
2012-01-20 06:54:30 +00:00
David Schultz
7bde21b1f9 Add an implementation of fenv.h intended for platforms that lack an FPU and
use softfloat.

Thanks to Ian Lepore for testing and debugging this patch.  The fenv
regression tests pass (at least for Ian's arm chip) with this change.
2012-01-16 04:09:17 +00:00
Ulrich Spörlein
bf3f9db657 Convert files to UTF-8 and add some copyright markers where missing. 2012-01-07 16:13:56 +00:00
David Chisnall
e481b86384 Expose the unimplemented libm functions in the math.h header. This allows C++'s <cmath> to work without the compiler complaining that the C++ versions are calling implicitly-declared functions. You will still get a linker error when they are called. OpenBSD 5.0 claims to fully implement the C99 <math.h> stuff, so might be worth investigating...
Reviewed by:	das
Approved by:	dim (mentor)
2011-11-12 19:55:48 +00:00
David Schultz
1f52bcdde4 Minor corrections and clarifications regarding exceptions. 2011-10-21 14:23:59 +00:00
David Schultz
0c7e4d5fc1 Fix a regression introduced in r226371: When the high part of x*y
exactly cancels with z, return the low part of x*y instead of
discarding it.
2011-10-21 06:30:43 +00:00
David Schultz
bc23acdc32 Fix a corner case: tan(large + Inf i) == NaN + NaN i. 2011-10-21 06:30:16 +00:00
David Schultz
c6df46bafb Improved handling of large x in ccosh{,f}():
- Handle cases where exp(x) would overflow, but ccosh(x) ~= exp(x) / 2
  shouldn't.
- Use the ccosh(x) ~= exp(x) / 2 approximation to simplify the calculation
  when x is large.

Similarly for csinh().  Also fixed the return value of csinh(-Inf +- 0i).
2011-10-21 06:29:32 +00:00
David Schultz
d4657ac799 Use __ldexp_exp() to simplify things and improve accuracy for x near
the overflow threshold.
2011-10-21 06:28:47 +00:00
David Schultz
12188b77a2 The cexp() and {,c}{cos,sin}h functions all need to be able to compute
exp(x) scaled down by some factor, and the challenge is doing this
accurately when exp(x) would overflow.  This change replaces all of
the tricks we've been using with common __ldexp_exp() and
__ldexp_cexp() routines that handle all the scaling.

bde plans to improve on this further by moving the guts of exp() into
k_exp.c and handling the scaling in a more direct manner.  But the
current approach is simple and adequate for now.
2011-10-21 06:27:56 +00:00
David Schultz
f2ea2b9d27 Use STRICT_ASSIGN() to ensure that the compiler doesn't screw things
up by storing x in a wider type than it's supposed to.

Submitted by:	bde
2011-10-21 06:26:38 +00:00
David Schultz
cd24d79843 Per IEEE754r, pow(1, y) is 1 even if y is NaN, and pow(-1, +-Inf) is 1.
MFC after:	2 weeks
2011-10-21 06:26:07 +00:00
David Schultz
741ae1d017 Bugfix: feenableexcept() and fedisableexcept() should just return the
old exception mask, not mask | ~FE_ALL_EXCEPT.

MFC after:	2 weeks
2011-10-21 06:25:31 +00:00
David Schultz
715a43edd0 It's no longer accurate to say that math.h "constitute[s] the C math
library," since complex.h, tgmath.h, and fenv.h are also part of the
math library.  Replace the outdated sentence with some references to
the other parts.
2011-10-17 06:10:32 +00:00
David Schultz
3daee1d6c3 Add c{cos,sin,tan}{,h}{,f} functions. This is joint work with
bde and kargl.
2011-10-17 05:41:03 +00:00
Eitan Adler
36daf0495a - change "is is" to "is" or "it is"
- change "the the" to "the"

Approved by:	lstewart
Approved by:	sahil (mentor)
MFC after:	3 days
2011-10-16 14:30:28 +00:00
David Schultz
5d9fefacf2 Use #include "fenv.h" instead of #include <fenv.h>. This makes it
more convenient to compile the math library by itself.

Requested by:	bde
2011-10-16 05:37:56 +00:00
David Schultz
0ddbe81e4f Fix some non-standard variable declarations. 2011-10-16 05:37:20 +00:00
David Schultz
4ce31585dc Optimize the case of pure imaginary arguments. Calls like this are
common, e.g., in DFT implementations.

Discussed with:	bde, kargl
2011-10-16 05:37:01 +00:00
David Schultz
e62394861d Move the macros GET_LDBL_EXPSIGN() and SET_LDBL_EXPSIGN() into
math_private.h, so they can be used elsewhere in the math library.
2011-10-16 05:36:39 +00:00
David Schultz
eccf8b3a30 Remove an unused variable. 2011-10-16 05:36:23 +00:00
David Schultz
e595c01fd8 Remove some unnecessary initializations.
Obtained from:	DragonFlyBSD
2011-10-15 07:00:28 +00:00
David Schultz
b052ec9065 Various changes to improve the accuracy and speed of log{2,10}{,f}.
- Rename __kernel_log() to k_log1p().
- Move some of the work that was previously done in the kernel log into
  the callers.  This enables further refactoring to improve accuracy or
  speed, although I don't recall the details.
- Use extra precision when adding the final scaling term, which improves
  accuracy.
- Describe and work around compiler problems that break some of the
  multiprecision calculations.

A fix for a small bug is also included:
- Add a special case for log*(1).  This is needed to ensure that log*(1) == +0
  instead of -0, even when the rounding mode is FE_DOWNWARD.

Submitted by:	bde
2011-10-15 05:23:28 +00:00
David Schultz
5ebf26f6a1 Style fixes and updates to comments.
Submitted by:	bde
2011-10-15 05:00:56 +00:00
David Schultz
15d5d9deb1 Don't define FP_FAST_FMA on sparc64; with the recent fixes, fma() is
no longer "fast" on sparc64.  (It really wasn't to begin with, since
the old implementation was using long doubles, and long doubles are
emulated in software on sparc64.)
2011-10-15 04:24:54 +00:00
David Schultz
3d4dfde48a Add INSERT_WORD64 and EXTRACT_WORD64 macros for use in s_fma.c. 2011-10-15 04:22:55 +00:00
David Schultz
90a83ac60a Replace two lines accidentally removed in r226218. Thanks to bde
for noticing this.
2011-10-15 04:17:20 +00:00
David Schultz
ec95083025 Fix a double-rounding bug in fma{,f,l}. The bug would occur in
round-to-nearest mode when the result, rounded to twice machine
precision, was exactly halfway between two machine-precision
values.  The essence of the fix is to simulate a "sticky bit" in
the pathological cases, which is how hardware implementations
break the ties.

MFC after:	1 month
2011-10-15 04:16:58 +00:00
David Schultz
7792754061 Refactor this code by introducing separate functions to handle the
extra-precision add and multiply operations. This simplifies future
work but shouldn't result in any functional change.
2011-10-11 05:17:45 +00:00
David Schultz
d78e594bc9 Provide external definitions of all of the standardized functions in
fenv.h that are currently inlined.

The definitions are provided in fenv.c via 'extern inline'
declaractions.  This assumes the compiler handles 'extern inline' as
specified in C99, which has been true under FreeBSD since 8.0.

The goal is to eventually remove the 'static' keyword from the inline
definitions in fenv.h, so that non-inlined references all wind up
pointing to the same external definition like they're supposed to.
I am deferring the second step to provide a window where
newly-compiled apps will still link against old math libraries.
(This isn't supported, but there's no need to cause undue breakage.)

Reviewed by:    stefanf, bde
2011-10-10 15:43:09 +00:00
Steve Kargl
edfca01f71 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)
2011-06-19 17:07:58 +00:00
Ben Laurie
5f301949ef Fix clang warnings.
Approved by:	philip (mentor)
2011-06-18 13:56:33 +00:00
Steve Kargl
9aa461b570 Clean up the unneeded cpp macro INLINE_REM_PIO2L.
Reviewed by:	das
Approved by:	das (mentor)
2011-05-30 19:41:28 +00:00
Steve Kargl
c273267e83 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
2011-04-29 23:13:43 +00:00
Steve Kargl
1cd0ec03d6 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
2011-03-12 19:37:35 +00:00
Steve Kargl
44bc6ce420 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)
2011-03-12 17:03:41 +00:00
Steve Kargl
dfe5233b59 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
2011-03-12 16:50:39 +00:00
David Schultz
d132715684 Add cexp() to the complex(3) manpage. Thanks to bde for pointing out
that I missed this.
2011-03-07 08:54:20 +00:00
David Schultz
d76098d531 Remove part of an uncommitted change that snuck into the last commit. 2011-03-07 08:46:14 +00:00
David Schultz
97a539be1f Convert log10f() to use __kernel_log(), which is more accurate and simpler. 2011-03-07 03:12:08 +00:00
David Schultz
acda0929b2 Convert log10() to use __kernel_log(), which is more accurate and simpler. 2011-03-07 03:11:27 +00:00
David Schultz
f3732b5aaa Add cexp() and cexpf().
Reviewed by:	bde (earlier version)
2011-03-07 03:09:24 +00:00
Rebecca Cran
6bccea7c2b Fix typos - remove duplicate "the".
PR:	bin/154928
Submitted by:	Eitan Adler <lists at eitanadler.com>
MFC after: 	3 days
2011-02-21 09:01:34 +00:00
Murray Stokely
11725f2142 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
2011-02-20 05:29:00 +00:00
David Schultz
b5209b6228 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
2011-02-10 07:38:38 +00:00
David Schultz
5ffd745ec2 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
2011-02-10 07:38:13 +00:00
David Schultz
e044d80d08 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
2011-02-10 07:37:50 +00:00
David Schultz
b775d18789 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
2011-02-10 07:37:29 +00:00
Konstantin Belousov
9f2e8bdff3 Remove duplicate .note.GNU-stack section declaration.
Reported by:	arundel
2011-02-04 21:54:06 +00:00
Konstantin Belousov
8997563c9a Add section .note.GNU-stack for assembly files used by 386 and amd64. 2011-01-07 16:13:12 +00:00
David Schultz
f353f24867 Another minor nit: Make sure the constant here is a float so the compiler
doesn't promote the entire expression to double.
2010-12-07 03:29:36 +00:00
David Schultz
63687c8b08 Fix various nits in style and comments that were pointed out by bde.
Code changes verified with md5.
2010-12-07 02:19:15 +00:00
David Schultz
177668d11f Add log2() and log2f(). 2010-12-05 22:11:22 +00:00
David Schultz
e7780530fa 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).
2010-12-05 22:11:03 +00:00
David Schultz
9bb4c87bb1 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
2010-12-03 00:05:49 +00:00
Ulrich Spörlein
e61ffaea2a 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
2010-11-13 10:54:10 +00:00
Dimitry Andric
b3b74eb4e1 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)
2010-10-01 20:14:36 +00:00
Warner Losh
fca5b4f4d9 This is exactly the same as the .else, so remove it. 2010-09-13 04:23:23 +00:00
Warner Losh
165721e3b1 MFtbemd: Move to using MACHINE_CPUARCH, now that it is safe. 2010-09-13 01:44:56 +00:00
Nathan Whitehorn
b12277d1d4 Repair some build breakage introduced in r211725 and garbage collect some
code made obsolete in the same commit.
2010-08-28 15:03:11 +00:00
Warner Losh
25faff346c MFtbemd:
Prefer MACHNE_CPUARCH to MACHINE_ARCH in most contexts where you want
to test of all the CPUs of a given family conform.
2010-08-23 22:24:11 +00:00
Joel Dahl
c2025a7660 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>
2010-08-16 15:18:30 +00:00
Nathan Whitehorn
3c766d1c8f powerpc64 floating-point is identical to powerpc, so use the same
code on both architectures.
2010-07-10 14:40:57 +00:00
David Schultz
efd0f253c2 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>
2010-06-12 17:32:05 +00:00
Ulrich Spörlein
df26753b39 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.
2010-06-09 07:31:32 +00:00
Ulrich Spörlein
e3303e900e mdoc: spell macros correctly, there's no need for the backslash escape 2010-06-02 10:20:38 +00:00
Ulrich Spörlein
25c8b07fc2 mdoc: Garbage collect unused/unneeded macros 2010-05-27 13:56:27 +00:00
Ulrich Spörlein
0b31f1f731 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
2010-05-13 12:08:11 +00:00
Ulrich Spörlein
0afc94c17a 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
2010-05-13 12:07:55 +00:00
Ulrich Spörlein
63d46d1d5e Fix several typos in macros or macro misusage.
Found by:	make manlint
Reviewed by:	ru
Approved by:	philip (mentor)
2010-03-12 10:01:06 +00:00
Konstantin Belousov
60d818ef9c 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
2010-02-03 20:23:47 +00:00
Ed Schouten
a4f326ddd5 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/
2009-06-11 13:59:51 +00:00
Ed Schouten
b492f2899f Use ISO C99 style inline semantics in msun.
Because we use ISO C99 nowadays, we can just get rid of enforcing
GNU89-style inlining.
2009-06-03 08:16:34 +00:00
Attilio Rao
9235ed7199 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
2009-05-25 14:37:10 +00:00
David Schultz
fbaabc11ae Namespace: scalb() is withdrawn from POSIX. 2009-03-14 18:58:53 +00:00
David Schultz
24863729f9 Eliminate __real__ and __imag__ gccisms. 2009-03-14 18:24:15 +00:00
David Schultz
ea8e257984 C99 TC2 now wants FP_FAST_FMA* to be defined to 1, if the macros are
defined at all. See also: defect report #223.
2009-02-07 05:41:24 +00:00
David Schultz
4630140ce2 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
2009-01-13 05:13:20 +00:00
David Schultz
4984f138f8 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
2009-01-08 06:12:03 +00:00
Marcel Moolenaar
74aed9855d 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
2008-12-23 22:20:59 +00:00
David Schultz
941ab616c7 Remove some unused variables.
Reported by:	Intel C Compiler
2008-08-08 00:21:27 +00:00
David Schultz
65faab72b3 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.
2008-08-08 00:15:16 +00:00
David Schultz
5e9470f1d5 Implement cproj{,f,l}(). 2008-08-07 15:07:48 +00:00
David Schultz
ad5e21e24b 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.
2008-08-07 14:39:56 +00:00
David Schultz
e142bd0ec9 Fix some style bogosity from fdlibm. 2008-08-03 17:49:05 +00:00
David Schultz
e034558322 Minor improvements:
- Improve the order of some tests.
- Fix style.

Submitted by:	bde
2008-08-03 17:39:54 +00:00
David Schultz
9d7d093689 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.
2008-08-02 19:17:00 +00:00
David Schultz
1192a80ed1 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
2008-08-02 03:56:22 +00:00
David Schultz
8862f666ad 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.
2008-08-01 01:24:25 +00:00
David Schultz
17303c626f Add implementations of acosl(), asinl(), atanl(), atan2l(),
and cargl().

Reviewed by:			bde
sparc64 testing resources from:	remko
2008-07-31 22:41:26 +00:00
David Schultz
d84f6378b9 Set WARNS=1.
I believe I've committed all the bits necessary to make this compile
on all supported architectures. :crosses fingers:
2008-07-31 20:11:37 +00:00
David Schultz
71da0e146b The high part of the mantissa is 64 bits on sparc64. 2008-07-31 20:09:47 +00:00
David Schultz
16608a810d 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.
2008-07-31 19:57:50 +00:00
David Schultz
6684d0ace3 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).
2008-07-18 02:18:34 +00:00
David Schultz
018b56d61f Fix a typo in the cosl() prototype. 2008-06-28 01:43:24 +00:00
David Schultz
f2dc0faad0 Implement fmodl.
Document fmodl and fix some errors in the fmod manpage.
2008-06-19 22:39:53 +00:00
Oleksandr Tymoshenko
41ea62f986 Symbol.map is handled by cpp, so use C-style comments
Approved by:	cognet (mentor)
2008-05-03 21:16:08 +00:00
Warner Losh
62c97aefc5 Add mips support to libm, from mips2-jnpr perforce branch. 2008-04-26 12:20:29 +00:00
David Schultz
92a1a6e169 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.
2008-04-03 06:14:51 +00:00
David Schultz
8087c515ab Remove a (bogus) remnant of debugging this on sparc64. 2008-03-31 13:11:45 +00:00
David Schultz
074fb64d9a Add assembly versions of remquol() and remainderl(). 2008-03-30 21:21:53 +00:00
David Schultz
c7392feecc Hook remquol() and remainderl() up to the build. 2008-03-30 20:48:02 +00:00
David Schultz
a2e5f27559 Implement remainderl() as a wrapper around remquol(). The extra work
remquol() performs to compute the quotient is negligible.
2008-03-30 20:47:42 +00:00
David Schultz
cef56f9d6d Implement remquol() based on remquo(). 2008-03-30 20:47:26 +00:00
David Schultz
511dd36b32 Implement csqrtl(). 2008-03-30 20:07:15 +00:00
David Schultz
84c1c0a1ca Hook hypotl() and cabsl() up to the build. 2008-03-30 20:03:46 +00:00
David Schultz
01a13522ad Document hypotl().
Submitted by:	Steve Kargl <sgk@troutmask.apl.washington.edu>
2008-03-30 20:03:29 +00:00
David Schultz
a641fc76eb Alias hypotl() and cabsl() for platforms where long double is the same
as double.
2008-03-30 20:03:06 +00:00
David Schultz
2264157a42 Implement cabsl() in terms of hypotl().
Submitted by:	Steve Kargl <sgk@troutmask.apl.washington.edu>
2008-03-30 20:02:03 +00:00
David Schultz
d23166b015 Implement hypotl(). This is bde's conversion of fdlibm hypot(), with minor
fixes for ld128 by me.
2008-03-30 20:01:50 +00:00
Bruce Evans
42ee187c3c 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.
2008-03-30 18:07:12 +00:00
Bruce Evans
c0c7ddd3a8 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.)
2008-03-30 17:28:27 +00:00
Bruce Evans
f94997c8d7 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).
2008-03-30 17:17:42 +00:00
David Schultz
a1af0d70da Include math.h for the fmaf() prototype. 2008-03-29 16:38:29 +00:00
David Schultz
ee0730e61e Fix some rather obscene code that has ambiguous if...if...else...
constructs in it.
2008-03-29 16:37:59 +00:00
David Schultz
3e13dd37ff 1 << 47 needs to be written 1ULL << 47. 2008-03-02 20:16:55 +00:00
David Schultz
e43c8f6acc Hook up sqrtl() to the build. 2008-03-02 01:48:17 +00:00
David Schultz
c6f56f9f41 MD implementations of sqrtl(). 2008-03-02 01:48:08 +00:00
David Schultz
c6a4447b64 MI implementation of sqrtl(). This is very slow and should
be overridden when hardware sqrt is available.
2008-03-02 01:47:58 +00:00
Bruce Evans
a278d99026 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.
2008-02-28 16:22:36 +00:00
Bruce Evans
e822ea5b2a 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.
2008-02-25 22:19:17 +00:00
Bruce Evans
c32951b16e 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.
2008-02-25 18:28:58 +00:00
Bruce Evans
70d818a20e 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.
2008-02-25 13:33:20 +00:00
Bruce Evans
0d1564b6c7 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.
2008-02-25 11:43:20 +00:00
Rafal Jaworowski
56ae1bed48 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
2008-02-24 19:22:53 +00:00
Bruce Evans
60a50c2585 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
2008-02-23 12:53:21 +00:00
Bruce Evans
dbf10e45c4 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).
2008-02-22 18:43:23 +00:00
Bruce Evans
7c1b5e7953 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.
2008-02-22 17:26:24 +00:00
Bruce Evans
43590b1517 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.
2008-02-22 15:55:14 +00:00
Bruce Evans
0ddfa46b44 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.
2008-02-22 14:11:03 +00:00
Bruce Evans
f839bac29c 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.
2008-02-22 11:59:05 +00:00
Bruce Evans
a7aa8cc980 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.
2008-02-22 10:04:53 +00:00
Bruce Evans
33843eef65 Fix rintl() on signaling NaNs and unsupported formats. 2008-02-22 09:21:14 +00:00
David Schultz
5aa554c7e5 s/rcsid/__FBSDID/ 2008-02-22 02:30:36 +00:00
David Schultz
fab324dfa4 Remove an unused variable. 2008-02-22 02:27:34 +00:00
David Schultz
7cd50f4d94 Eliminate some warnings. 2008-02-22 02:26:51 +00:00
Bruce Evans
f21d26becb 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.
2008-02-19 15:42:46 +00:00
Bruce Evans
9e9d3bc9f1 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).
2008-02-19 15:30:58 +00:00
Bruce Evans
9ce8756044 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().
2008-02-19 12:54:14 +00:00
David Schultz
345241c5e0 Document return values better. 2008-02-18 19:02:49 +00:00
David Schultz
71c11dd528 Add tgammaf() as a simple wrapper around tgamma(). 2008-02-18 17:27:11 +00:00
Bruce Evans
be396b71c1 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.
2008-02-18 15:39:52 +00:00
Bruce Evans
19a9e1bb1c Fix a typo which broke k_tanl.c on !(amd64 || i386). 2008-02-18 14:09:41 +00:00
Bruce Evans
38662c9698 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.
2008-02-18 14:02:12 +00:00
David Schultz
842d1d5c98 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.
2008-02-17 21:53:19 +00:00
David Schultz
234b60cd97 Hook up sinl(), cosl(), and tanl() to the build. 2008-02-17 07:33:51 +00:00
David Schultz
8e77cc6431 Add implementations of sinl(), cosl(), and tanl().
Submitted by:	Steve Kargl <sgk@apl.washington.edu>
2008-02-17 07:33:12 +00:00
David Schultz
f869a8c5f3 Documentation for sinl(), cosl(), and tanl(). 2008-02-17 07:32:44 +00:00
David Schultz
61f955827d 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)
2008-02-17 07:32:31 +00:00
David Schultz
de336b0c5e 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
2008-02-17 07:32:14 +00:00
David Schultz
079299f710 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.
2008-02-17 07:31:59 +00:00
Bruce Evans
63b4a1f80c 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
2008-02-15 07:01:40 +00:00
Bruce Evans
5014f8ded4 Oops, the weak reference for ceill(), floorl() and truncl() was in the
wrong file.  This broke arm and powerpc.

Reported by:	grehan
2008-02-14 15:10:34 +00:00
Bruce Evans
3365b45e5e 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().
2008-02-14 13:44:03 +00:00
Bruce Evans
b4437c3d32 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.
2008-02-14 10:23:51 +00:00
Bruce Evans
011cbae1fe 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().
2008-02-14 09:42:24 +00:00
Bruce Evans
e7c95ee5fe s_ceill.c
s_floorl.c
s_truncl.c
2008-02-13 17:38:16 +00:00
Bruce Evans
74d68da630 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.
2008-02-13 16:56:52 +00:00
Bruce Evans
6597187205 Fix the C version of ceill(x) for -1 < x <= -0 in all rounding modes.
The result should be -0, but was +0.
2008-02-13 15:22:53 +00:00
Bruce Evans
f01bfe5c6d 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.
2008-02-13 10:44:44 +00:00
Bruce Evans
828f7b4a82 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.
2008-02-13 08:36:13 +00:00
Bruce Evans
02ef796d23 Use hardware remainder on amd64 since it is 5 to 10 times faster than
software remainder and is already used for remquo().
2008-02-13 06:01:48 +00:00
Bruce Evans
a2ddfa5ea7 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().
2008-02-12 17:11:36 +00:00
Bruce Evans
51f86873af 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.
2008-02-11 05:20:02 +00:00
Bruce Evans
52453261e9 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.
2008-02-09 12:53:15 +00:00
Bruce Evans
6d656800db Fix a comment about coefficients and expand a related one. 2008-02-09 10:36:07 +00:00
Bruce Evans
fbe8fb4d7b 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.
2008-02-08 01:45:52 +00:00
Bruce Evans
aa7c7c47cf 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().
2008-02-08 01:22:13 +00:00
Bruce Evans
a00672cff9 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().
2008-02-07 09:42:19 +00:00
Bruce Evans
a373e66b85 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.
2008-02-07 03:17:05 +00:00
Bruce Evans
ce56838fdc 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.
2008-02-06 06:35:21 +00:00
David Schultz
b134ea7211 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.
2008-01-28 01:19:07 +00:00
Bruce Evans
fc84b771b4 Fix a harmless type error in 1.9. 2008-01-25 21:09:21 +00:00
Bruce Evans
f2a1477818 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.
2008-01-21 13:46:21 +00:00
Bruce Evans
85c309021f 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).
2008-01-20 04:09:44 +00:00
Bruce Evans
a9b721d6b2 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.
2008-01-19 22:50:50 +00:00
Bruce Evans
5b62c3808e 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).
2008-01-19 22:05:14 +00:00
Bruce Evans
684217d889 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.
2008-01-19 21:37:14 +00:00
Bruce Evans
fa7fdac725 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().
2008-01-19 20:02:55 +00:00
Bruce Evans
0814af48f7 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.
2008-01-19 18:13:21 +00:00
Bruce Evans
6a876b92fb 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().
2008-01-19 16:37:57 +00:00
David Schultz
86c2e0c047 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.
2008-01-18 22:19:04 +00:00
David Schultz
3d2cc91218 Hook up exp2l() and related docs to the build. 2008-01-18 21:43:10 +00:00
David Schultz
5526551600 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.
2008-01-18 21:43:00 +00:00
David Schultz
968b39e3b9 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.
2008-01-18 21:42:46 +00:00
Bruce Evans
1880ccbd79 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).
2008-01-17 17:02:11 +00:00
David Schultz
a2d171e440 Optimize this a bit better.
Submitted by:	bde (although these aren't all of his changes)
2008-01-15 23:31:24 +00:00
David Schultz
d3f9671a7d Implement rintl(), nearbyintl(), lrintl(), and llrintl().
Thanks to bde@ for feedback and testing of rintl().
2008-01-14 02:12:07 +00:00
David Schultz
73b2958b94 - 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.
2008-01-11 04:18:25 +00:00
David Schultz
45310fdb5d 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
2007-12-20 03:16:55 +00:00
David Schultz
58c9a67ed7 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
2007-12-20 00:06:33 +00:00
David Schultz
7cd4a83267 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
2007-12-18 23:46:32 +00:00
David Schultz
0ba1fd2f72 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
2007-12-18 01:15:20 +00:00
Bruce Evans
ccef8c4fcb Oops, the previous commit was not needed -- the file was committed but
not checked out due to my checkout error.
2007-12-17 18:21:23 +00:00
Bruce Evans
a18b106ffc 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.
2007-12-17 18:12:06 +00:00
Bruce Evans
9ed67737f2 Don't try to build s_nanl.c before it is committed. 2007-12-17 13:20:38 +00:00
David Schultz
6821aba9e5 Add logbl(3) to libm. 2007-12-17 03:53:38 +00:00
David Schultz
3be0479b4c Document the fact that we have nan(3) now, and make some minor clarifications
in other places.
2007-12-17 01:04:43 +00:00
David Schultz
4b6b574455 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.
2007-12-16 21:19:28 +00:00
David Schultz
3c27af2a44 1. Add csqrt{,f}(3).
2. Put carg{,f}(3) under the FBSD_1.1 namespace where it belongs
   (requested by kan@)
2007-12-15 08:39:03 +00:00
David Schultz
aaf70b2314 Implement and document csqrt(3) and csqrtf(3). 2007-12-15 08:38:44 +00:00
David Schultz
ce448a2e74 Update the standards section, and make a minor clarification about the
return value of sqrt.
2007-12-14 07:53:09 +00:00
David Schultz
80f974729f Typo in previous commit 2007-12-14 03:08:10 +00:00
David Schultz
39ebc398b6 Symbol.map additions for carg and cargf. (They're in C99, so I didn't
add a new version for them.)
2007-12-14 03:06:50 +00:00
David Schultz
9768c3fea8 s/C90/C99/ 2007-12-12 23:50:00 +00:00
David Schultz
367d55260f Add a "STANDARDS" section. 2007-12-12 23:49:40 +00:00
David Schultz
205bd64894 Implement carg(3) and cargf(3).
Rotting in an old src tree since: March 2005
2007-12-12 23:43:51 +00:00
Bruce Evans
b5e547df33 Oops, back out previous commit since it was backwards to a wrong branch. 2007-06-14 05:57:13 +00:00
Bruce Evans
d382c5ebb4 MFC: 1.11: fix the threshold for (not) using the simple Taylor approximation. 2007-06-14 05:51:00 +00:00
Bruce Evans
a8a2e00ebf 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>
2007-06-11 07:48:52 +00:00
Bruce Evans
20a990117d 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.
2007-05-29 07:13:07 +00:00
Daniel Eischen
419ecd5dee Bump library versions in preparation for 7.0.
Ok'd by:	kan
2007-05-21 02:49:08 +00:00
Daniel Eischen
00fb440c1a 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.
2007-05-13 14:12:40 +00:00
Bruce Evans
9698b3b564 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.
2007-05-02 16:54:22 +00:00
Bruce Evans
e95cc9b700 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).
2007-05-02 15:24:49 +00:00
Bruce Evans
dd936b27fc 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.
2007-05-02 13:49:28 +00:00
Daniel Eischen
5f864214bb Use C comments since we now preprocess these files with CPP. 2007-04-29 14:05:22 +00:00