Commit Graph

325 Commits

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

Unlike for cosf and sinf, the old accuracy was close to being inadequate
-- the polynomial for double precision has a max error of 0.014 ulps
and nearly this small an error is needed.  The new accuracy is also a
bit small, but exhaustive checking shows that even the old accuracy
was enough.  The increased accuracy reduces the maximum relative error
in the final result on amd64 -O1 from 0.9588 ulps to 0.9044 ulps.
2005-11-10 17:43:49 +00:00
bde
35f17c1d45 Detach k_rem_pio2f.c from the build since it is now unused. It is a libm
internal so this shouldn't cause version problems.
2005-11-06 17:59:40 +00:00
bde
e016ebc9a1 Use a 53-bit approximation to pi/2 instead of a 33+53 bit one for the
special case pi/4 <= |x| < 3*pi/4.  This gives a tiny optimization (it
saves 2 subtractions, which are scheduled well so they take a whole 1
cycle extra on an AthlonXP), and simplifies the code so that the
following optimization is not so ugly.

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

The 53-bit approximation to pi/2 is good enough for pi/4 <= |x| <
3*pi/4 because after losing up to 24 bits to subtraction, we still
have 29 bits of precision and only need 25 bits.  Even with only 5
extra bits, it is possible to get perfectly rounded results starting
with the reduced x, since if x is nearly a multiple of pi/2 then x is
not near a half-way case and if x is not nearly a multiple of pi/2
then we don't lose many bits.  With our intentionally imperfect rounding
we get the same results for cosf(), sinf() and tanf() as without this
optimization.
2005-11-06 17:48:02 +00:00
bde
0ec5232d0c The logb() functions are not just ieee754 "test" functions, but are
standard in C99 and POSIX.1-2001+.  They are also not deprecated, since
apart from being standard they can handle special args slightly better
than the ilogb() functions.

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

Add history for the logb() functions and ilogbl().  Fix history for
ilogb().
2005-11-06 12:18:27 +00:00
bde
ea9959fde3 Moved the optimization for tiny x from __kernel_tan[f](x) to tan[f](x)
so that it can be faster for tiny x and avoided for reduced x.

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

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

This change makes no difference to the results for sinf() (since it
only changes the algorithm for the 4 special args and the results for
those happen not to change), but it changes lots of results for sin().
Exhaustive testing is impossible for sin(), but exhaustive testing
for sinf() (relative to a version with the old algorithm and a fixed
cutoff) shows that the changes in the error are either reductions or
from 0.5-epsilon ulps to 0.5+epsilon ulps.  The new method just uses
some extra terms in approximations so it tends to give more accurate
results, and there are apparently no problems from having extra
accuracy.  On amd64 with -O1, on all float args the error range in ulps
is reduced from (0.500, 0.665] to [0.335, 0.500) in 24168 cases and
increased from 0.500-epsilon to 0.500+epsilon in 24 cases.  Non-
exhaustive testing by ucbtest shows no differences.
2005-11-02 14:01:45 +00:00
bde
728b935c7f Updated the comment about the optimization for tiny x (the previous
commit moved it).  This includes a comment that the "kernel" sine no
longer works on arg -0, so callers must now handle this case.  The kernel
sine still works on all other tiny args; without the optimization it is
just a little slower on these args.  I intended it to keep working on
all tiny args, but that seems to be impossible without losing efficiency
or accuracy.  (sin(x) ~ x * (1 + S1*x**2 + ...) would preserve -0, but
the approximation must be written as x + S1*x**3 + ... for accuracy.)
2005-11-02 13:06:49 +00:00
bde
481c63491c Removed dead code for handling tan[f]() on odd multiples of pi/2. This
case never occurs since pi/2 is irrational so no multiple of it can
be represented as a float and we have precise arg reduction so we never
end up with a remainder of 0 in the "kernel" function unless the
original arg is 0.

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

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

Index: k_tan.c
===================================================================
RCS file: /home/ncvs/src/lib/msun/src/k_tan.c,v
retrieving revision 1.10
diff -r1.10 k_tan.c
88,90c88
< 			if (((ix | low) | (iy + 1)) == 0)
< 				return one / fabs(x);
< 			else {
---
> 			{
2005-11-02 06:45:21 +00:00
bde
d568fc134a Fixed some of the silliness related to rev.1.8. In 1.8, "double" in
a declaration was not translated to "float" although bit fiddling on
double variables was translated.  This resulted in garbage being put
into the low word of one of the doubles instead of non-garbage being
put into the only word of the intended float.  This had no effect on
any result because:
- with doubles, the algorithm for calculating -1/(x+y) is unnecessarily
  complicated.  Just returning -1/((double)x+y) would work, and the
  misdeclaration gave something like that except for messing up some
  low bits with the bit fiddling.
- doubles have plenty of bits to spare so messing up some of the low
  bits is unlikely to matter.
- due to other bugs, the buggy code is reached for a whole 4 args out
  of all 2**32 float args.  The bug fixed by 1.8 only affects a small
  percentage of cases and a small percentage of 4 is 0.  The 4 args
  happen to cause no problems without 1.8, so they are even less likely
  to be affected by the bug in 1.8 than average args; in fact, neither
  1.8 nor this commit makes any difference to the result for these 4
  args (and thus for all args).

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

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

The fix in 1.8 works with the reduced cutoff in all cases despite the
bug in it.  It changes the result in 84492 cases altogether to fix the
5826 broken cases.  Fixing the fix by translating "double" to "float"
changes the result in 42 cases relative to 1.8.  In 24 cases the
(absolute) error is increased and in 18 cases it is reduced, but it
remains less than 1 ulp in all cases.
2005-11-02 05:37:31 +00:00
bde
bce05c8c60 Fixed spelling of remquof() in its prototype. 2005-10-30 12:34:58 +00:00
bde
eb7e930697 Fixed some comments added in rev.1.5.
The log message for 1.5 said that some small (one or two ulp) inaccuracies
were fixed, and a comment implied that the critical change is to switch
the rounding mode to to-nearest, with a switch of the precision to
extended at no extra cost.  Actually, the errors are very large (ucbtest
finds ones of several hundred ulps), and it is the switch of the
precision that is critical.

Another comment was wrong about NaNs being handled sloppily.
2005-10-30 12:21:02 +00:00
bde
26610cfe9b Implement inline functions to give the complex result x+I*y from float
or double args x and y.  x+I*y cannot be used directly yet due to compiler
bugs.

Submitted by:	Steve Kargl <sgk@troutmask.apl.washington.edu>
2005-10-29 17:14:11 +00:00
bde
bbfb40721e Use double precision to simplify and optimize arg reduction for small
and medium size args too: instead of conditionally subtracting a float
17+24, 17+17+24 or 17+17+17+24 bit approximation to pi/2, always
subtract a double 33+53 bit one.  The float version is now closer to
the double version than to old versions of itself -- it uses the same
33+53 bit approximation as the simplest cases in the double version,
and where the float version had to switch to the slow general case at
|x| == 2^7*pi/2, it now switches at |x| == 2^19*pi/2 the same as the
double version.

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

This speeds up cosf(x) on large args by a factor of about 2.  Correct
arg reduction on large args is still inherently very slow, so hopefully
these args rarely occur in practice.  There is much more efficiency
to be gained by using double precision to speed up arg reduction on
medium and small float args.
2005-10-29 08:15:29 +00:00
bde
8e62cdabe0 Use fairly optimal minimax polynomials for __kernel_cosf() and
__kernel_sinf().  The old ones were the double-precision polynomials
with coefficients truncated to float.  Truncation is not a good way
to convert minimax polynomials to lower precision.  Optimize for
efficiency and use the lowest-degree polynomials that give a relative
error of less than 1 ulp -- degree 8 instead of 14 for cosf and degree
9 instead of 13 for sinf.  For sinf, the degree 8 polynomial happens
to be 6 times more accurate than the old degree 14 one, but this only
gives a tiny amount of extra accuracy in results -- we just need to
use a a degree high enough to give a polynomial whose relative accuracy
in infinite precision (but with float coefficients) is a small fraction
of a float ulp (fdlibm generally uses 1/32 for the small fraction, and
the fraction for our degree 8 polynomial is about 1/600).

The maximum relative errors for cosf() and sinf() are now 0.7719 ulps
and 0.7969 ulps, respectively.
2005-10-28 13:36:58 +00:00
bde
96c89ee304 Use a better algorithm for reducing the error in __kernel_cos[f]().
This supersedes the fix for the old algorithm in rev.1.8 of k_cosf.c.

I want this change mainly because it is an optimization.  It helps
make software cos[f](x) and sin[f](x) faster than the i387 hardware
versions for small x.  It is also a simplification, and reduces the
maximum relative error for cosf() and sinf() on machines like amd64
from about 0.87 ulps to about 0.80 ulps.  It was validated for cosf()
and sinf() by exhaustive testing.  Exhaustive testing is not possible
for cos() and sin(), but ucbtest reports a similar reduction for the
worst case found by non-exhaustive testing.  ucbtest's non-exhaustive
testing seems to be good enough to find problems in algorithms but not
maximum relative errors when there are spikes.  E.g., short runs of
it find only 3 ulp error where the i387 hardware cos() has an error
of about 2**40 ulps near pi/2.
2005-10-26 12:36:18 +00:00
bde
d6cdac5f7a More fixes for arg reduction near pi/2 on systems with broken assignment
to floats (mainly i386's).  All errors of more than 1 ulp for float
precision trig functions were supposed to have been fixed; however,
compiling with gcc -O2 uncovered 18250 more such errors for cosf(),
with a maximum error of 1.409 ulps.

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

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

Extra precision is harmful mainly for implementing extra precision in
software.  We want to represent x+y as w+r where both "+" operations
are in infinite precision and r is tiny compared with w.  There is a
standard algorithm for this (Knuth (1981) 4.2.2 Theorem C), and fdlibm
uses this routinely, but the algorithm requires w and r to have the
same precision as x and y.  w is just x+y (calculated in the same
finite precision as x and y), and r is a tiny correction term.  The
i386 gcc bugs tend to give extra precision in w, and then using this
extra precision in the calculation of r results in the correction
mostly staying in w and being missing from r.  There still tends to
be no problem if the result is a simple expression involving w and r
-- modulo spills, w keeps its extra precision and r remains the right
correction for this wrong w.  However, here we want to pass w and r
to extern functions.  Extra precision is not retained in function args,
so w gets fixed up, but the change to the tiny r is tinier, so r almost
remains as a wrong correction for the right w.
2005-10-25 12:13:37 +00:00
bde
5931c79161 Moved the optimization for tiny x from __kernel_{cos,sin}[f](x) to
{cos_sin}[f](x) so that x doesn't need to be reclassified in the
"kernel" functions to determine if it is tiny (it still needs to be
reclassified in the cosine case for other reasons that will go away).

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

Fixed various wrong constants for the cutoff for "tiny".  For cosine,
the cutoff is when x**2/2! == {FLT or DBL}_EPSILON/2.  We round down
to an integral power of 2 (and for cos() reduce the power by another
1) because the exact cutoff doesn't matter and would take more work
to determine.  For sine, the exact cutoff is larger due to the ration
of terms being x**2/3! instead of x**2/2!, but we use the same cutoff
as for cosine.  We now use a cutoff of 2**-27 for double precision and
2**-12 for single precision.  2**-27 was used in all cases but was
misspelled 2**27 in comments.  Wrong and sloppy cutoffs just cause
missed optimizations (provided the rounding mode is to nearest --
other modes just aren't supported).
2005-10-24 14:08:36 +00:00
bde
86f27343be Fixed range reduction for large multiples of pi/2 on systems with
broken assignment to floats (e.g., i386 with gcc -O, but not amd64 or
ia64; i386 with gcc -O0 worked accidentally).

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

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

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

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

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

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

          amd64     ia64      i386-O0   i386-O1
	  ------    ------    ------    ------
cosf:     0.8681    0.8681    0.7927    0.5650
sinf:     0.8733    0.8610    0.7849    0.5651
tanf:     0.9708    0.9329    0.9329    0.7035
2005-10-11 07:56:05 +00:00
bde
32945bd185 Fixed range reduction near (but not very near) medium-sized multiples
of pi/2 (1 line) and expand a comment about related magic (many lines).

The bug was essentially the same as for the +-pi/2 case (a mistranslated
mask), but was smaller so it only significantly affected multiples
starting near +-13*pi/2.  At least on amd64, for cosf() on all 2^32
float args, the bug caused 128 errors of >= 1 ulp, with a maximum error
of 1.2393 ulps.
2005-10-10 20:02:02 +00:00
bde
6210e62129 Fix numerous errors of >= 1 ulp for cosf(x) and sinf(x) (1 line)
and add a comment about related magic (many lines)).

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

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

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

The error bounds in the above and in my other recent commit messages
are for amd64.  Extra precision for floats on i386's accidentally masks
this bug, but only if k_cosf.c is compiled with -O.  Although the extra
precision helps here, this is accidental and depends on longstanding
gcc precision bugs (not clipping extra precision on assignment...),
and the gcc bugs are mostly avoided by compiling without -O.  I now
develop libm mainly on amd64 systems to simplify error detection and
debugging.
2005-10-09 21:07:23 +00:00
bde
485c06b5bb Oops, the last-minute optimization in rev.1.8 wasn't a good idea. The
17+17+24 bit pi/2 must only be used when subtraction of the first 2
terms in it from the arg is exact.  This happens iff the the arg in
bits is one of the 2**17[-1] values on each side of (float)(pi/2).

Revert to the algorithm in rev.1.7 and only fix its threshold for using
the 3-term pi/2.  Use the threshold that maximizes the number of values
for which the 3-term pi/2 is used, subject to not changing the algorithm
for comparing with the threshold.  The 3-term pi/2 ends up being used
for about half of its usable range (about 64K values on each side).
2005-10-09 04:29:08 +00:00
bde
39c85cfe13 Fixed syntax error (a missing brace) in previous commit. 2005-10-08 22:55:36 +00:00
bde
13f78201bb Fixed range reduction near (but not very near) +-pi/2. A bug caused
a maximum error of 2.905 ulps for cosf(), but the algorithm for cosf()
is good for < 1 ulps and happens to give perfect rounding (< 0.5 ulps)
near +-pi/2 except for the bug.  The extra relative errors for tanf()
were similar (slightly larger).  The bug didn't affect sinf() since
sinf'(+-pi/2) is 0.

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

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

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

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

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

The bug is apparently from mistranslation of the magic 32 in 53-32-1.
The number of bits allowed to cancel is not critical and we use 32 for
double precision because it allows efficient classification using a
32-bit comparison.  For float precision, we must use an explicit mask,
and there are fewer bits so there is less margin for error in their
allocation.  The 32 got reduced to 4 but should have been reduced
almost in proportion to the reduction of mantissa bits.
2005-10-08 22:43:55 +00:00
bde
41d865435a Fixed aliasing bugs in TRUNC() by using the fdlibm macros for access
to doubles as bits.  fdlibm-1.1 had similar aliasing bugs, but these
were fixed by NetBSD or Cygnus before a modified version of fdlibm was
imported in 1994.  TRUNC() is only used by tgamma() and some
implementation-detail functions.  The aliasing bugs were detected by
compiling with gcc -O2 but don't seem to have broken tgamma() on i386's
or amd64's.  They broke my modified version of tgamma().

Moved the definition of TRUNC() to mathimpl.h so that it can be fixed
in one place, although the general version is even slower than necessary
because it has to operate on pointers to volatiles to handle its arg
sometimes being volatile.  Inefficiency of the fdlibm macros slows
down libm generally, and tgamma() is a relatively unimportant part of
libm.  The macros act as if on 32-bit words in memory, so they are
hard to optimize to direct actions on 64-bit double registers for
(non-i386) machines where this is possible.  The optimization is too
hard for gcc on amd64's, and declaring variables as volatile makes it
impossible.
2005-09-19 11:28:19 +00:00
das
665ea151e9 Add a missing ldexpf() alias for amd64.
Noticed by:	bz@, tjr@
2005-09-12 20:54:00 +00:00
kensmith
f97f77429f Bump the shared library version number of all libraries that have not
been bumped since RELENG_5.

Reviewed by:	ru
Approved by:	re (not needed for commit check but in principle...)
2005-07-22 17:19:05 +00:00
ru
0c80b11e62 Markup nit.
Approved by:	re (blanket)
2005-06-16 21:56:03 +00:00
ru
3f3ef36f49 Fixed compile warning.
Approved by:	re (blanket)
2005-06-16 21:55:45 +00:00
ru
38fc91ca96 Assorted markup fixes.
Approved by:	re
2005-06-15 19:04:04 +00:00
deischen
5d3cf26519 Prevent these functions from using stack outside of their frame.
Reported by:	Marc Olzheim <marcolz at stack dot nl>
OK'd by:	das
2005-05-06 15:44:20 +00:00
stefanf
16dd1d18f5 Revert the last change, the conversion from long double to double can raise
unwanted underflow exceptions.

Pointed out by:	das
2005-04-28 19:45:55 +00:00
stefanf
13322eaf9b Use double additions to raise the inexact exception to work around problems
with long double addition on sparc64.
2005-04-22 09:57:55 +00:00
stefanf
77782516eb Fix raising the inexact exception (FE_INEXACT) if the result differs from the
argument.

Noticed by:	das
2005-04-22 08:30:33 +00:00
ache
60e7539065 Fix truncl.3 MLINKS 2005-04-17 19:57:52 +00:00
das
9c49c2a65a More optimized math functions. 2005-04-16 21:12:55 +00:00
das
a9fd105354 Implement truncl() based on floorl(). 2005-04-16 21:12:47 +00:00
das
15bd306d7a Add roundl(), lroundl(), and llroundl(). 2005-04-08 01:24:08 +00:00
das
4df744a471 These files should include s_lround.c instead of s_lrint.c.
This only matters for efficiency, not for correctness.
2005-04-08 00:52:27 +00:00
das
9977034544 Fix a (coincidentally harmless) bug. 2005-04-08 00:52:16 +00:00
das
bb5d9cb768 Fix a long-standing bug in k_rem_pio2(), which led to large errors when
tanf() was called with big arguments close to multiples of pi/2.

Reported by:	ucbtest via bde
2005-04-05 23:27:47 +00:00
das
1db6f984bc Build exp2(), exp2f(), and related documentation. 2005-04-05 02:57:39 +00:00
das
1324b71e9d Document exp2() and exp2f(), and make other minor tweaks and updates. 2005-04-05 02:57:28 +00:00
das
bf28283937 Implement exp2() and exp2f(). 2005-04-05 02:57:15 +00:00
das
da9b203aaf Implement and document remquo() and remquof(). 2005-03-25 04:40:44 +00:00
das
79b831e3a1 Fix the double rounding problem with subnormals, and
remove the XXX comments, which no longer apply.
2005-03-18 02:27:59 +00:00
das
9f8ee2b273 Add missing prototypes for fma() and fmaf(), and remove an inaccurate
comment.
2005-03-18 01:47:42 +00:00
das
fdf53809bb Make the fenv.h routines work for programs that use SSE for
floating-point arithmetic on i386.  Now I'm going to make excuses
for why this code is kinda scary:

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

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

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

Also, I un-inlined some of the functions that grew larger as a result
of this change, moving them from fenv.h to fenv.c.
2005-03-17 22:21:46 +00:00
das
5b7d321e53 Spell 'fedisableexcept' correctly. 2005-03-16 22:34:14 +00:00
das
bfdcd78bfc Document feenableexcept(), fedisableexcept(), and fegetexcept(). 2005-03-16 19:04:28 +00:00
das
6448887f3b Replace fegetmask() and fesetmask() with feenableexcept(),
fedisableexcept(), and fegetexcept().  These two sets of routines
provide the same functionality.  I implemented the former as an
undocumented internal interface to make the regression test easier to
write.  However, fe(enable|disable|get)except() is already part of
glibc, and I would like to avoid gratuitous differences.  The only
major flaw in the glibc API is that there's no good way to report
errors on processors that don't support all the unmasked exceptions.
2005-03-16 19:03:46 +00:00