Commit Graph

536 Commits

Author SHA1 Message Date
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