Commit Graph

400 Commits

Author SHA1 Message Date
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
Warner Losh
ee7093a640 Remove California Regent's clause 3, per letter 2007-01-09 01:02:06 +00:00
David Schultz
9abb1ff616 Implement modfl(). 2007-01-07 07:54:21 +00:00
David Schultz
8185b32b5a Fix a problem relating to fesetenv() clobbering i387 register stack.
Details: As a side-effect of restoring a saved FP environment,
fesetenv() overwrites the tag word, which indicates which i387
registers are in use.  Normally this isn't a problem because
the calling convention requires the register stack to be empty
on function entry and exit.  However, fesetenv() is inlined, so we
need to tell gcc explicitly that the i387 registers get clobbered.

PR:	85101
2007-01-06 21:46:23 +00:00
David Schultz
93e0663877 Fix a cut-and-paste-o. 2007-01-06 21:23:20 +00:00
David Schultz
829d55ac9c Correctly handle NaN. 2007-01-06 21:22:57 +00:00
David Schultz
9fa229fc8d Correctly handle inf/nan. This routine is currently unused because we
seem to have assembly versions for all architectures, but it can't
hurt to fix it.
2007-01-06 21:22:38 +00:00
David Schultz
6642c3fa74 Remove modf from libm's symbol map. It's actually in libc for
hysterical raisins.
2007-01-06 21:18:17 +00:00
David Schultz
3cb636ce18 Remove an unneeded fnstcw instruction.
Noticed by:	bde
2007-01-05 07:15:26 +00:00
David Schultz
cc0d85b680 Remove a note pertaining to the Alpha. 2007-01-05 07:14:26 +00:00
Bruce Evans
fae6222bdb Moved __BEGIN_DECLS up a little so that it covers __test_sse() and C++
isn't broken,

PR:		104425
2006-10-14 20:35:56 +00:00
Ruslan Ermilov
2b46c64c9c Remove alpha left-overs. 2006-08-22 08:03:01 +00:00
Bruce Evans
d79d610d9c Fixed the threshold for using the simple Taylor approximation.
In e_log.c, there was just a off-by-1 (1 ulp) error in the comment
about the threshold.  The precision of the threshold is unimportant,
but the magic numbers in the code are easier to understand when the
threshold is described precisely.

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

The optimization of using the simple Taylor approximation for args
near a power of 2 is dubious since it only applies to a relatively
small proportion of args, but if it is done then doing it 2048 times
as often _may_ be more efficient.  (My benchmarks show unexplained
dependencies on the data that increase with further optimizations
in this area.)
2006-07-07 04:33:08 +00:00
Bruce Evans
fe72622ebe Fixed tanh(-0.0) on ia64 and optimizeed tanh(x) for 2**-55 <= |x| <
2**-28 as a side effect, by merging with the float precision version
of tanh() and the double precision version of sinh().

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

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

Updated some comments and fixed bugs in others (ranges for half-open
intervals mostly had the open end backwards, and there were nearby style
bugs).
2006-07-05 22:59:33 +00:00
Bruce Evans
3454a5a101 Removed the optimized asm versions of scalb() and scalbf(). These
functions are only for compatibility with obsolete standards.  They
shouldn't be used, so they shouldn't be optimized.  Use the generic
versions instead.

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

The generic versions have slightly different error handling:
- the asm versions blindly round the second parameter to a (floating
  point) integer and proceed, while the generic versions return NaN
  if this rounding changes the value.  POSIX permits both behaviours
  (these functions are XSI extensions and the behaviour for a bogus
  non-integral second parameter is unspecified).   Apart from this
  and the bug in scalbf(), the behaviour of the generic versions seems
  to be identical.  (I only exhusatively tested
  generic_scalbf(1.0F, anyfloat) == asm_scalb(1.0F, anyfloat).  This
  covers many representative corner cases involving NaNs and Infs but
  doesn't test exception flags.  The brokenness of scalbf() showed up
  as weird behaviour after testing just 7 integer cases sequentially.)
2006-07-05 20:06:42 +00:00
Bruce Evans
8eca9455de Backed out rev.1.10. It tried to implement ldexpf() as a weak reference
to scalbf(), but ldexpf() cannot be implemented in that way since the
types of the second parameter differ.  ldexpf() can be implemented as
a weak or strong reference to scalbnf() (*) but that was already done
long before rev.1.10 was committed.  The old implementation uses a
reference, so rev.1.10 had no effect on applications.  The C files for
the scalb() family are not used for amd64 or i386, so rev.1.10 had even
less effect for these arches.

(*) scalbnf() raises the radix to the given exponent, while ldexpf()
raises 2 to the given exponent.  Thus the functions are equivalent
except possibly for their error handling iff the radix is 2.  Standards
more or less require identical error handling.  Under FreeBSD, the
functions are equivalent except for more details being missing in
scalbnf()'s man page.
2006-07-05 02:16:29 +00:00
Daniel Eischen
d7eda46253 Add symbol versioning to libm. 2006-03-27 23:59:45 +00:00
Bruce Evans
fd2891004d Oops, on amd64 (and probably on all non-i386 systems), the previous
commit broke the 2**24 cases where |x| > DBL_MAX/2.  There are exponent
range problems not just for denormals (underflow) but for large values
(overflow).  Doubles have more than enough exponent range to avoid the
problems, but I forgot to convert enough terms to double, so there was
an x+x term which was sometimes evaluated in float precision.

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

Exlain the problem better in comments.
2006-01-05 09:18:48 +00:00
Bruce Evans
4bb9780353 Use double precision internally to optimize cbrtf(), and change the
algorithm for the second step significantly to also get a perfectly
rounded result in round-to-nearest mode.  The resulting optimization
is about 25% on Athlon64's and 30% on Athlon XP's (about 25 cycles
out of 100 on the former).

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

Using specifically double precision, there is enough extra precision
for the third step to give enough precision for perfect rounding to
float precision provided the previous steps are accurate to 16 bits.
(They were accurate to 12 bits, which was almost minimal for imperfect
rounding in the old version but would be more than enough for imperfect
rounding in this version (9 bits would be enough now).)  I couldn't
find any significant time optimizations from optimizing the previous
steps, so I decided to optimize for accuracy instead.  The second step
needed a division although a previous commit optimized it to use a
polynomial approximation for its main detail, and this division dominated
the time for the second step.  Use the same Newton's method for the
second step as for the third step since this is insignificantly slower
than the division plus the polynomial (now that Newton's method only
needs 1 division), significantly more accurate, and simpler.  Single
precision would be precise enough for the second step, but doesn't
have enough exponent range to handle denormals without the special
grouping of terms (as in previous versions) that requires another
division, so we use double precision for both the second and third
steps.
2006-01-05 07:57:31 +00:00
Bruce Evans
5776f433ab Extract the high and low words together. With gcc-3.4 on uniformly
distributed non-large args, this saves about 14 of 134 cycles for
Athlon64s and about 5 of 199 cycles for AthlonXPs.

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

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

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

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

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

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

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

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

Lower order Pade approximations are not quite accurate enough for
double precision but are plenty for float precision.  A lower order
Remez rational approximation might be enough for double precision too.
However, rational approximations inherently require an extra division,
and polynomial approximations work well for 1/cbrt(r) at r = 1, so I
plan to switch to using the latter.  There are some technical
complications that tend to cost a division in another way.
2005-12-15 16:23:22 +00:00
Bruce Evans
ec761d7501 Optimize by not doing excessive conversions for handling the sign bit.
This gives an optimization of between 9 and 22% on Athlons (largest
for cbrt() on amd64 -- from 205 to 159 cycles).

We extracted the sign bit and worked with |x|, and restored the sign
bit as the last step.  We avoided branches to a fault by using accesses
to FP values as bits to clear and restore the sign bit.  Avoiding
branches is usually good, but the bit access macros are not so good
(especially for setting FP values), and here they always caused pipeline
stalls on Athlons.  Even using branches would be faster except on args
that give perfect branch misprediction, since only mispredicted branches
cause stalls, but it possible to avoid touching the sign bit in FP
values at all (except to preserve it in conversions from bits to FP
not related to the sign bit).  Do this.  The results are identical
except in 2 of the 3 unsupported rounding modes, since all the
approximations use odd rational functions so they work right on strictly
negative values, and the special case of -0 doesn't use an approximation.
2005-12-13 20:17:23 +00:00
Bruce Evans
7d5a4821ba Fixed some especially horrible style bugs (indentation that is neither
KNF nor fdlibmNF combined with multiple statements per line).
2005-12-13 18:22:00 +00:00
Bruce Evans
af7f99131d Added comments about the magic behind
<cbrt(x) in bits> ~= <x in bits>/3 + BIAS.
Keep the large comments only in the double version as usual.

Fixed some style bugs (mainly grammar and spelling errors in comments).
2005-12-11 19:51:30 +00:00
Bruce Evans
288a8c86cb Fixed the unexpectedly large maximum error after the previous commit.
It was because I forgot to translate the part of the double precision
algorithm that chops t so that t*t is exact.  Now the maximum error
is the same as for double precision (almost exactly 2.0/3 ulps).
2005-12-11 17:58:14 +00:00
Bruce Evans
6de073b4ef Fixed all 502518670 errors of more than 1 ulp for cbrtf() on amd64.
The maximum error was 3.56 ulps.

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

In float precision, any reasonable algorithm works for the cleanup
step.  Use the same algorithm as for double precision, although this
is much more than enough and is a significant pessimization, and don't
optimize or simplify anything using double precision to implement the
float case, so that the whole double precision algorithm can be verified
in float precision.  A maximum error of 0.667 ulps is claimed for cbrt()
and the max for cbrtf() using the same algorithm shouldn't be different,
but the actual max for cbrtf() on amd64 is now 0.9834 ulps.  (On i386
-O1 the max is 0.5006 (down from < 0.7) due to extra precision.)
2005-12-11 13:22:01 +00:00
Bruce Evans
1a787460ba Fixed some magic numbers.
The threshold for not being tiny was too small.  Use the usual 2**-12
threshold.  As for sinhf, use a different method (now the same as for
sinhf) to set the inexact flag for tiny nonzero x so that the larger
threshold works, although this method is imperfect.  As for sinhf,
this change is not just an optimization, since the general code that
we fell into has accuracy problems even for tiny x.  On amd64, avoiding
it fixes tanhf on 2*13495596 args with errors of between 1 and 1.3
ulps and thus reduces the total number of args with errors of >= 1 ulp
from 37533748 to 5271278; the maximum error is unchanged at 2.2 ulps.

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

I started testing with nonstandard rounding modes, and verified that
the chosen thresholds work for all modes modulo problems not related
to thresholds.  The best thresholds are not very dependent on the mode,
at least for tanhf.
2005-12-11 11:40:55 +00:00
David E. O'Brien
9b39b7cba6 "Create" ldexpf for non-i386 architectures.
Submitted by:	Steve Kargl <sgk@troutmask.apl.washington.edu>
2005-12-06 20:12:38 +00:00
Bruce Evans
0f06be5a4d Fixed the approximation to pio4. pio4_hi must be pio2_hi/2 since it
shares its low half with pio2_hi.  pio2_hi is rounded down although
rounding to nearest would be a tiny bit better, so pio4_hi must be
rounded down too.  It was rounded to nearest, which happens to be
different in float precision but the same in double precision.

This fixes about 13.5 million errors of more than 1 ulp in asinf().
The largest error was 2.81 ulps on amd64 and 2.57 ulps on i386 -O1.
Now the largest error is 0.93 ulps on amd65 and 0.67 ulps on i386 -O1.
2005-12-04 13:52:46 +00:00
Bruce Evans
d48ea9753c For log1pf(), fixed the approximations to sqrt(2), sqrt(2)-1 and
sqrt(2)/2-1.  For log1p(), fixed the approximation to sqrt(2)/2-1.

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

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

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

In logf(), all the approximations made sense, but the approximation
to sqrt(2)/2-1 was 2 ulps too small (a tiny amount, since we compare
with a granularity of 2**32 ulps), so the algorithm was broken in 2
cases.  The result was broken in 1 case.  This is fixed by rounding
up the approximation to sqrt(2)/2-1 by 2**32 ulps, so 2**32 cases are
now handled a little differently (still correctly according to my
assertion that the approximations don't need to be very accurate, but
this has not been checked).
2005-12-04 12:30:44 +00:00
Bruce Evans
669152498a Use the usual volatile hack to trick gcc into clipping any extra precision
on assignment.

Extra precision on i386's broke hi+lo decomposition in the usual way.
It caused all except 1 of the 62343 errors of more than 1 ulp for
log1pf() on i386's with gcc -O [-fno-float-store].
2005-12-04 08:57:54 +00:00
Bruce Evans
00b1756b1e Fixed fdlibm[+cygnus] logbf() and logb() on denormals. Adjustment
according to the highest nonzero bit in a denormal was missing.

fdlibm ilogbf() and ilogb() have always had the adjustment, but only
use a small part of their method for handling denormals; use the
normalization method in log[f]() for the main part.
2005-12-03 11:57:19 +00:00
Bruce Evans
1186054263 Restored removal of the special handling needed for a result of +-0.
It was lost in rev.1.9.  The log message for rev.1.9 says that the
special case of +-0 is handled twice, but it was only handled once,
so it became unhandled, and this happened to break half of the cases
that return +-0:
- round-towards-minus-infinity:  0   <  x < 1:  result was -0 not  0
- round-to-nearest:             -0.5 <= x < 0:  result was  0 not -0
- round-towards-plus-infinity:  -1   <  x < 0:  result was  0 not -0
- round-towards-zero:           -1   <  x < 0:  result was  0 not -0
2005-12-03 09:00:29 +00:00
Bruce Evans
3fc5a433e9 Simplified the fix in rev.1.3. Instead of using long double for
TWO52[sx] to trick gcc into correctly converting TWO52[sx]+x to double
on assignment to "double w", force a correct assignment by assigning
to *(double *)&w.  This is cleaner and avoids the double rounding
problem on machines that evaluate double expressions in double
precision.  It is not necessary to convert w-TWO52[sx] to double
precision on return as implied in the comment in rev.1.3, since
the difference is exact.
2005-12-03 07:38:35 +00:00
Bruce Evans
7441377544 Fixed rint(x) in the following cases:
(1) In round-to-nearest mode, on all machines, fdlibm rint() never
    worked for |x| = n+0.75 where n is an even integer between 262144
    and 524286 inclusive (2*131072 cases).  To avoid double rounding
    on some machines, we begin by adjusting x to a value with the 0.25
    bit not set, essentially by moving the 0.25 bit to a lower bit
    where it works well enough as a guard, but we botched the adjustment
    when log2(|x|) == 18 (2*2**52 cases) and ended up just clearing
    the 0.25 bit then.  Most subcases still worked accidentally since
    another lower bit serves as a guard.  The case of odd n worked
    accidentally because the rounding goes the right way then.  However,
    for even n, after mangling n+0.75 to 0.5, rounding gives n but the
    correct result is n+1.
(2) In round-towards-minus-infinity mode, on all machines, fdlibm rint()
    never for x = n+0.25 where n is any integer between -524287 and
    -262144 inclusive (262144 cases).  In these cases, after mangling
    n+0.25 to n, rounding gives n but the correct result is n-1.
(3) In round-towards-plus-infinity mode, on all machines, fdlibm rint()
    never for x = n+0.25 where n is any integer between 262144 and
    524287 inclusive (262144 cases).  In these cases, after mangling
    n+0.25 to n, rounding gives n but the correct result is n+1.

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

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

The double rounding problem is relatively rare, so the botched adjustment
can be fixed for most machines by removing the entire adjustment.  This
would be a wrong fix (using it is 1 of the bugs in rev.1.9 of s_rintf.c)
since fdlibm is supposed to be generic, but it works in the following cases:
- on all machines that evaluate double expressions in double precision,
  provided either long double has the same precision as double (alpha,
  and i386's with precision forced to double) or my earlier fix to use
  a long double 2**52 is modified to avoid using long double precision.
- on all machines that evaluate double expressions in many more than 11
  bits of extra precision.  The 1 bit of extra precision in the example
  is the worst case.  With N bits of extra precision, it sufices to
  adjust the bit N bits below the 0.5 bit.  For N >= about 52 there is
  no such bit so the adjustment is both impossible and unnecessary.  The
  fix in rev.1.9 of s_rintf.c apparently depends on corresponding magic
  in float precision: on all supported machines N is either 0 or >= 24,
  so double rounding doesn't occur in practice.
- on all machines that don't use fdlibm rint*() (i386's).
So under FreeBSD, the double rounding problem only affects amd64 now, but
should only affect i386 in future (when double expressions are evaluated
in long double precision).
2005-12-03 07:23:30 +00:00
Bruce Evans
5792e54aa9 Fixed roundf(). The following cases never worked in FreeBSD:
- in round-towards-minus-infinity mode, on all machines, roundf(x) never
  worked for 0 < |x| < 0.5 (2*0x3effffff cases in all, or almost half of
  float space).  It was -0 for 0 < x < 0.5 and 0 for -0.5 < x < 0, but
  should be 0 and -0, respectively.  This is because t = ceilf(|x|) = 1
  for these args, and when we adjust t from 1 to 0 by subtracting 1, we
  get -0 in this rounding mode, but we want and expected to get 0.
- in round-towards-minus-infinity, round towards zero and round-to-nearest
  modes, on machines that evaluate float expressions in float precision
  (most machines except i386's), roundf(x) never worked for |x| =
  <float value immediately below 0.5> (2 cases in all).  It was +-1 but
  should have been +-0.  This is because t = ceilf(|x|) = 1 for these
  args, and when we try to classify |x| by subtracting it from 1 we
  get an unexpected rounding error -- the result is 0.5 after rounding
  to float in all 3 rounding modes, so we we have forgotten the
  difference between |x| and 0.5 and end up returning the same value
  as for +-0.5.

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

Attempted to fix round() and roundl() by cloning the fix for roundf().
This has only been tested for round(), only on args representable as
floats.  Double expressions are evaluated in double precision even on
i386's, so round(0.5-epsilon) was broken even on i386's.  roundl()
must be completely broken on i386's since long double precision is not
really supported.  There seem to be no other dependencies on the
precision.
2005-12-02 13:45:06 +00:00
Bruce Evans
f4b01a9edf Rearranged the polynomial evaluation to reduce dependencies, as in
k_tanf.c but with different details.

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

On Athlons, for uniformly distributed args in [-2pi, 2pi], this gives
an optimization of about 4 cycles (10%) in most cases (13% for sinf()
on AXP, but 0% for cosf() with gcc-3.3 -O1 on AXP).  The best case
(sinf() with gcc-3.4 -O1 -fcaller-saves on A64) now takes 33-39 cycles
(was 37-45 cycles).  Hardware sinf takes 74-129 cycles.  Despite
being fine tuned for Athlons, the optimization is even larger on
some other arches (about 15% on ia64 (pluto2) and 20% on alpha (beast)
with gcc -O2 -fomit-frame-pointer).
2005-11-30 11:51:17 +00:00
Bruce Evans
8d3b309bad Fixed cosf(x) when x is a "negative" NaNs. I broke this in rev.1.10.
cosf(x) is supposed to return something like x when x is a NaN, and
we actually fairly consistently return x-x which is normally very like
x (on i386 and and it is x if x is a quiet NaN and x with the quiet bit
set if x is a signaling NaN.  Rev.1.10 broke this by normalising x to
fabsf(x).  It's not clear if fabsf(x) is should preserve x if x is a NaN,
but it actually clears the sign bit, and other parts of the code depended
on this.

The bugs can be fixed by saving x before normalizing it, and using the
saved x only for NaNs, and using uint32_t instead of int32_t for ix
so that negative NaNs are not misclassified even if fabsf() doesn't
clear their sign bit, but gcc pessimizes the saving very well, especially
on Athlon XPs (it generates extra loads and stores, and mixes use of
the SSE and i387, and this somehow messes up pipelines).  Normalizing
x is not a very good optimization anyway, so stop doing it.  (It adds
latency to the FPU pipelines, but in previous versions it helped except
for |x| <= 3pi/4 by simplifying the integer pipelines.)  Use the same
organization as in s_sinf.c and s_tanf.c with some branches reordered.
These changes combined recover most of the performance of the unfixed
version on A64 but still lose 10% on AXP with gcc-3.4 -O1 but not with
gcc-3.3 -O1.
2005-11-30 06:47:18 +00:00
Bruce Evans
908801933a Fixed the hi+lo approximation to log(2). The normal 17+24 bit decomposition
that was used doesn't work normally here, since we want to be able to
multiply `hi' by the exponent of x _exactly_, and the exponent of x has
more than 7 significant bits for most denormal x's, so the multiplication
was not always exact despite a cloned comment claiming that it was.  (The
comment is correct in the double precision case -- with the normal 33+53
bit decomposition the exponent can have 20 significant bits and the extra
bit for denormals is only the 11th.)

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

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