s_cosf.c and s_sinf.c:
Use a non-bogus magic constant for the threshold of pi/4. It was 2 ulps
smaller than pi/4 rounded down, but its value is not critical so it should
be the result of natural rounding.
s_cosf.c and s_tanf.c:
Use a literal 0.0 instead of an unnecessary variable initialized to
[(float)]0.0. Let the function prototype convert to 0.0F.
Improved wording in some comments.
Attempted to improve indentation of comments.
number of branches.
Use a non-bogus magic constant for the threshold of pi/4. It was 2 ulps
smaller than pi/4 rounded down, but its value is not critical so it should
be the result of natural rounding. Use "<=" comparisons with rounded-
down thresholds for all small multiples of pi/4.
Cleaned up previous commit:
- use static const variables instead of expressions for multiples of pi/2
to ensure that they are evaluated at compile time. gcc currently
evaluates them at compile time but C99 compilers are not required
to do so. We want compile time evaluation for optimization and don't
care about side effects.
- use M_PI_2 instead of a magic constant for pi/2. We need magic constants
related to pi/2 elsewhere but not here since we just want pi/2 rounded
to double and even prefer it to be rounded in the default rounding mode.
We can depend on the cmpiler being C99ish enough to round M_PI_2 correctly
just as much as we depended on it handling hex constants correctly. This
also fixes a harmless rounding error in the hex constant.
- keep using expressions n*<value for pi/2> in the initializers for the
static const variables. 2*M_PI_2 and 4*M_PI_2 are obviously rounded in
the same way as the corresponding infinite precision expressions for
multiples of pi/2, and 3*M_PI_2 happens to be rounded like this, so we
don't need magic constants for the multiples.
- fixed and/or updated some comments.
define also, but res_config.h was not included into libc/net/name6.c.
So getipnodebyaddr() ignored the multiple PTRs.
PR: kern/88241
Submitted by: Dan Lukes <dan__at__obluda.cz>
MFC after: 3 days
The threshold for not being tiny was too small. Use the usual 2**-12
threshold. This change is not just an optimization, since the general
code that we fell into has accuracy problems even for tiny x. Avoiding
it fixes 2*1366 args with errors of more than 1 ulp, with a maximum
error of 1.167 ulps.
The magic number 22 is log(DBL_EPSILON)/2 plus slop. This is bogus
for float precision. Use 9 (~log(FLT_EPSILON)/2 plus less slop than
for double precision). The code for handling the interval
[2**-28, 9_was_22] has accuracy problems even for [9, 22], so this
change happens to fix errors of more than 1 ulp in about 2*17000
cases. It leaves such errors in about 2*1074000 cases, with a max
error of 1.242 ulps.
The threshold for switching from returning exp(x)/2 to returning
exp(x/2)^2/2 was a little smaller than necessary. As for coshf(),
This was not quite harmless since the exp(x/2)^2/2 case is inaccurate,
and fixing it avoids accuracy problems in 2*6 cases, leaving problems
in 2*19997 cases.
Fixed naming errors in pseudo-code in comments.
The threshold for not being tiny was confusing and too small. Use the
usual 2**-12 threshold and simplify the algorithm slightly so that
this threshold works (now use the threshold for sinhf() instead of one
for 1+expm1()). This is just a small optimization.
The magic number 22 is log(DBL_EPSILON)/2 plus slop. This is bogus
for float precision. Use 9 (~log(FLT_EPSILON)/2 plus less slop than
for double precision).
The threshold for switching from returning exp(x)/2 to returning
exp(x/2)^2/2 was a little smaller than necessary. This was not quite
harmless since the exp(x/2)^2/2 case is inaccurate. Fixing it happens
to avoid accuracy problems for 2*6 of the 2*151 args that were handled
by the exp(x)/2 case. This leaves accuracy problems for about 2*19997
args near the overflow threshold (~89); the maximum error there is
2.5029 ulps.
There are also accuracy probles for args in +-[0.5*ln2, 9] -- 2*188885
args with errors of more than 1 ulp, with a maximum error of 1.384 ulps.
Fixed a syntax error and naming errors in pseudo-code in comments.
specialized for float precision. The new polynomial has degree 8
instead of 14, and a maximum error of 2**-34.34 (absolute) instead of
2**-30.66. This doesn't affect the final error significantly; the
maximum error was and is about 0.8879 ulps on amd64 -01.
The fdlibm expf() is not used on i386's (the "optimized" asm version
is used), but probably should be since it was already significantly
faster than the asm version on athlons. The asm version has the
advantage of being more accurate, so keep using it for now.
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.
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.
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().
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.
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.)
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 {
---
> {
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.
rewritten, now timers created with same sigev_notify_attributes will
run in same thread, this allows user to organize which timers can
run in same thread to save some thread resource.
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.
