summaryrefslogtreecommitdiffstats
path: root/lib/msun/src
Commit message (Collapse)AuthorAgeFilesLines
...
* Fixed some especially horrible style bugs (indentation that is neitherbde2005-12-132-12/+14
| | | | KNF nor fdlibmNF combined with multiple statements per line).
* Added comments about the magic behindbde2005-12-112-14/+25
| | | | | | | <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).
* Fixed the unexpectedly large maximum error after the previous commit.bde2005-12-111-2/+2
| | | | | | 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).
* Fixed all 502518670 errors of more than 1 ulp for cbrtf() on amd64.bde2005-12-111-1/+13
| | | | | | | | | | | | | | | | | | | | | | | | | | 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.)
* Fixed some magic numbers.bde2005-12-111-8/+8
| | | | | | | | | | | | | | | | | | | | | | 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.
* "Create" ldexpf for non-i386 architectures.obrien2005-12-061-0/+2
| | | | Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
* Fixed the approximation to pio4. pio4_hi must be pio2_hi/2 since itbde2005-12-041-1/+1
| | | | | | | | | | | 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.
* For log1pf(), fixed the approximations to sqrt(2), sqrt(2)-1 andbde2005-12-042-8/+22
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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).
* Use the usual volatile hack to trick gcc into clipping any extra precisionbde2005-12-041-1/+1
| | | | | | | | 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].
* Fixed fdlibm[+cygnus] logbf() and logb() on denormals. Adjustmentbde2005-12-032-8/+20
| | | | | | | | 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.
* Restored removal of the special handling needed for a result of +-0.bde2005-12-031-0/+2
| | | | | | | | | | | 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
* Simplified the fix in rev.1.3. Instead of using long double forbde2005-12-031-8/+2
| | | | | | | | | | 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.
* Fixed rint(x) in the following cases:bde2005-12-031-0/+9
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | (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).
* Fixed roundf(). The following cases never worked in FreeBSD:bde2005-12-023-18/+18
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | - 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.
* Rearranged the polynomial evaluation to reduce dependencies, as inbde2005-11-302-9/+13
| | | | | | | | | | | | | | | | | | | | 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).
* Fixed cosf(x) when x is a "negative" NaNs. I broke this in rev.1.10.bde2005-11-301-11/+19
| | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fixed the hi+lo approximation to log(2). The normal 17+24 bit decompositionbde2005-11-301-4/+4
| | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Rearranged the polynomial evaluation some more to reduce dependencies.bde2005-11-281-8/+20
| | | | | | | | | | | | | | | | | | | | | Instead of echoing the code in a comment, try to describe why we split up the evaluation in a special way. The new optimization is mostly to move the evaluation of w = z*z later so that everything else (except z = x*x) doesn't have to wait for w. On Athlons, FP multiplication has a latency of 4 cycles so this optimization saves 4 cycles per call provided no new dependencies are introduced. Tweaking the other terms in to reduce dependencies saves a couple more cycles in some cases (more on AXP than on A64; up to 8 cycles out of 56 altogether in some cases). The previous version had a similar optimization for s = z*x. Special optimizations like these probably have a larger effect than the simple 2-way vectorization permitted (but not activated by gcc) in the old version, since 2-way vectorization is not enough and the polynomial's degree is so small in the float case that non-vectorizable dependencies dominate. On an AXP, tanf() on uniformly distributed args in [-2pi, 2pi] now takes 34-55 cycles (was 39-59 cycles).
* Fixed about 50 million errors of infinity ulps and about 3 million errorsbde2005-11-281-1/+1
| | | | | | | | | | | | | | | | | | | | | | | | | | of between 1.0 and 1.8509 ulps for lgammaf(x) with x between -2**-21 and -2**-70. As usual, the cutoff for tiny args was not correctly translated to float precision. It was 2**-70 but 2**-21 works. Not as usual, having a too-small threshold was worse than a pessimization. It was just a pessimization for (positive) args between 2**-70 and 2**-21, but for the first ~50 million (negative) args below -2**-70, the general code overflowed and gave a result of infinity instead of correct (finite) results near 70*log(2). For the remaining ~361 million negative args above -2**21, the general code gave almost-acceptable errors (lgamma[f]() is not very accurate in general) but the pessimization was larger than for misclassified tiny positive args. Now the max error for lgammaf(x) with |x| < 2**-21 is 0.7885 ulps, and speed and accuracy are almost the same for positive and negative args in this range. The maximum error overall is still infinity ulps. A cutoff of 2**-70 is probably wastefully small for the double precision case. Smaller cutoffs can be used to reduce the max error to nearly 0.5 ulps for tiny args, but this is useless since the general algrorithm for nearly-tiny args is not nearly that accurate -- it has a max error of about 1 ulp.
* Exploit skew-symmetry to avoid an operation: -sin(x-A) = sin(A-x). Thisbde2005-11-282-4/+4
| | | | | | | | gives a tiny but hopefully always free optimization in the 2 quadrants to which it applies. On Athlons, it reduces maximum latency by 4 cycles in these quadrants but has usually has a smaller effect on total time (typically ~2 cycles (~5%), but sometimes 8 cycles when the compiler generates poor code).
* Try to use the "proximity" (~) operator consistently in commentsbde2005-11-282-8/+10
| | | | | | | (x ~<= a, not x <= ~a). This got messed up in some places when the comments were moved from e_rem_pio2f.c. Added my (non-)copyright.
* Changed spelling of the request-to-inline macro name to match the changebde2005-11-282-4/+6
| | | | | | | | | of the function name. Added my (non-)copyright. In k_tanf.c, added the first set of redundant parentheses to control grouping which was claimed to be added in the previous commit.
* Use only double precision for "kernel" cosf and sinf (except forbde2005-11-286-71/+44
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | returning float). The functions are renamed from __kernel_{cos,sin}f() to __kernel_{cos,sin}df() so that misuses of them will cause link errors and not crashes. This version is an almost-routine translation with no special optimizations for accuracy or efficiency. The not-quite-routine part is that in __kernel_cosf(), regenerating the minimax polynomial with double precision coefficients gives a coefficient for the x**2 term that is not quite -0.5, so the literal 0.5 in the code and the related `hz' variable need to be modified; also, the special code for reducing the error in 1.0-x**2*0.5 is no longer needed, so it is convenient to adjust all the logic for the x**2 term a little. Note that without extra precision, it would be very bad to use a coefficient of other than -0.5 for the x**2 term -- the old version depends on multiplication by -0.5 being infinitely precise so as not to need even more special code for reducing the error in 1-x**2*0.5. This gives an unimportant increase in accuracy, from ~0.8 to ~0.501 ulps. Almost all of the error is from the final rounding step, since the choice of the minimax polynomials so that their contribution to the error is a bit less than 0.5 ulps just happens to give contributions that are significantly less (~.001 ulps). An Athlons, for uniformly distributed args in [-2pi, 2pi], this gives overall speed increases in the 10-20% range, despite giving a speed decrease of typically 19% (from 31 cycles up to 37) for sinf() on args in [-pi/4, pi/4].
* Minor cleanups and optimizations:bde2005-11-241-11/+5
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | - Remove dead code that I forgot to remove in the previous commit. - Calculate the sum of the lower terms of the polynomial (divided by x**5) in a single expression (sum of odd terms) + (sum of even terms) with parentheses to control grouping. This is clearer and happens to give better instruction scheduling for a tiny optimization (an average of about ~0.5 cycles/call on Athlons). - Calculate the final sum in a single expression with parentheses to control grouping too. Change the grouping from first_term + (second_term + sum_of_lower_terms) to (first_term + second_term) + sum_of_lower_terms. Normally the first grouping must be used for accuracy, but extra precision makes any grouping give a correct result so we can group for efficiency. This is a larger optimization (average 3-4 cycles/call or 5%). - Use parentheses to indicate that the C order of left to right evaluation is what is wanted (for efficiency) in a multiplication too. The old fdlibm code has several optimizations related to these. 2 involve doing an extra operation that can be done almost in parallel on some superscalar machines but are pessimizations on sequential machines. Others involve statement ordering or expression grouping. All of these except the ordering for the combining the sums of the odd and even terms seem to be ideal for Athlons, but parallelism is still limited so all of these optimizations combined together with the ones in this commit save only ~6-8 cycles (~10%). On an AXP, tanf() on uniformly distributed args in [-2pi, 2pi] now takes 39-59 cycles. I don't know of any more optimizations for tanf() short of writing it all in asm with very MD instruction scheduling. Hardware fsin takes 122-138 cycles. Most of the optimizations for tanf() don't work very well for tan[l](). fdlibm tan() now takes 145-365 cycles.
* Optimized by eliminating the special case for 0.67434 <= |x| < pi/4.bde2005-11-241-16/+7
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | A single polynomial approximation for tan(x) works in infinite precision up to |x| < pi/2, but in finite precision, to restrict the accumulated roundoff error to < 1 ulp, |x| must be restricted to less than about sqrt(0.5/((1.5+1.5)/3)) ~= 0.707. We restricted it a bit more to give a safety margin including some slop for optimizations. Now that we use double precision for the calculations, the accumulated roundoff error is in double-precision ulps so it can easily be made almost 2**29 times smaller than a single-precision ulp. Near x = pi/4 its maximum is about 0.5+(1.5+1.5)*x**2/3 ~= 1.117 double-precision ulps. The minimax polynomial needs to be different to work for the larger interval. I didn't increase its degree the old degree is just large enough to keep the final error less than 1 ulp and increasing the degree would be a pessimization. The maximum error is now ~0.80 ulps instead of ~0.53 ulps. The speedup from this optimization for uniformly distributed args in [-2pi, 2pi] is 28-43% on athlons, depending on how badly gcc selected and scheduled the instructions in the old version. The old version has some int-to-float conversions that are apparently difficult to schedule well, but gcc-3.3 somehow did everything ~10 cycles or ~10% faster than gcc-3.4, with the difference especially large on AXPs. On A64s, the problem seems to be related to documented penalties for moving single precision data to undead xmm registers. With this version, the speed is cycles is almost independent of the athlon and gcc version despite the large differences in instruction selection to use the FPU on AXPs and SSE on A64s.
* Use only double precision for "kernel" tanf (except for returning float).bde2005-11-233-29/+20
| | | | | | | | | | | | | | This is a minor interface change. The function is renamed from __kernel_tanf() to __kernel_tandf() so that misues of it will cause link errors and not crashes. This version is a routine translation with no special optimizations for accuracy or efficiency. It gives an unimportant increase in accuracy, from ~0.9 ulps to 0.5285 ulps. Almost all of the error is from the minimax polynomial (~0.03 ulps and the final rounding step (< 0.5 ulps). It gives strange differences in efficiency in the -5 to +10% range, with -O1 fairly consistently becoming faster and -O2 slower on AXP and A64 with gcc-3.3 and gcc-3.4.
* Simplified setiing up args for __kernel_rem_pio2(). We already have xbde2005-11-231-17/+9
| | | | | with a 24-bit fraction, so we don't need a loop to split it into up to 3 terms with 24-bit fractions.
* Quick fix for stack buffer overrun in rev.1.13. Oops. The prec == 1bde2005-11-231-4/+4
| | | | | | | | | | | | | | | | | | | | | arg to __kernel_rem_pio2() gives 53-bit (double) precision, not single precision and/or the array dimension like I thought. prec == 2 is used in e_rem_pio2.c for double precision although it is documented to be for 64-bit (extended) precision, and I just reduced it by 1 thinking that this would give the value suitable for 24-bit (float) precision. Reducing it 1 more to the documented value for float precision doesn't actually work (it gives errors of ~0.75 ulps in the reduced arg, but errors of much less than 0.5 ulps are needed; the bug seems to be in kernel_rem_pio2.c). Keep using a value 1 larger than the documented value but supply an array large enough hold the extra unused result from this. The bug can also be fixed quickly by increasing init_jk[0] in k_rem_pio2.c from 2 to 3. This gives behaviour identical to using prec == 1 except it doesn't create the extra result. It isn't clear how the precision bug affects higher precisions. 113-bit (quad) is the largest precision, so there is no way to use a large precision to fix it.
* Mess up the "kernel" float trig function .c files with ifdefs so thatbde2005-11-216-0/+25
| | | | | | | | | | | | they can be #included in other .c files to give inline functions, and use them to inline the functions in most callers (not in e_lgammaf_r.c). __kernel_tanf() is too large and complicated for gcc to inline very well. An athlons, this gives a speed increase under favourable pipeline conditions of about 10% overall (larger for AXP, smaller for A64). E.g., on AXP, sinf() on uniformly distributed args in [-2Pi, 2Pi] now takes 30-56 cycles; it used to take 45-61 cycles; hardware fsin takes 65-129.
* Use double precision to simplify and optimize a long division.bde2005-11-211-15/+1
| | | | | | | | | | | | | | | | On athlons, this gives a speedup of 10-20% for tanf() on uniformly distributed args in [-2Pi, 2Pi]. (It only directly applies for 43% of the args and gives a 16-20% speedup for these (more for AXP than A64) and this gives an overall speedup of 10-12% which is all that it should; however, it gives an overall speedup of 17-20% with gcc-3.3 on AXP-A64 by mysteriously effected cases where it isn't executed.) I originally intended to use double precision for all internals of float trig functions and will probably still do this, but benchmarking showed that converting to double precision and back is a pessimization in cases where a simple float precision calculation works, so it may be optimal to switch precisions only when using extra precision is much simpler.
* Restored a cleanup in rev.1.9 tthat was lost in rev.1.10.bde2005-11-201-2/+2
|
* Moved all the optimizations for |x| <= 9pi/2 frombde2005-11-194-67/+105
| | | | | | | | | | | | | | __ieee754_rem_pio2f() to its 3 callers and manually inline them. On Athlons, with favourable compiler flags and optimizations and favourable pipeline conditions, this gives a speedup of 30-40 cycles for cosf(), sinf() and tanf() on the range pi/4 < |x| <= 9pi/4, so thes functions are now signifcantly faster than the hardware trig functions in many cases. E.g., in a benchmark with uniformly distributed x in [-2pi, 2pi], A64 hardware fcos took 72-129 cycles and cosf() took 37-55 cycles. Out-of-order execution is needed to get both of these times. The optimizations in this commit apparently work more by removing 1 serialization point than by reducing latency.
* Minor cleanups:bde2005-11-173-24/+21
| | | | | | | | | | | | | | | 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.
* Rearranged the the optimizations for special cases to reduce the averagebde2005-11-171-42/+53
| | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fixed some magic numbers.bde2005-11-131-6/+6
| | | | | | | | | | | | | | | | | | | | | | | | 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.
* Fixed some magic numbers.bde2005-11-131-6/+6
| | | | | | | | | | | | | | | | | | | | | | | | 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.
* Imoproved comments for the minimax polynomial.bde2005-11-121-10/+11
| | | | | | Removed an unused variable. Fixed some wrong comments and some nearby misformatting.
* Tweaked the minimax polynomial and improved its comments.bde2005-11-121-5/+5
|
* Improved comments for the minimax polynomial.bde2005-11-121-4/+4
|
* As for the float trig functions, use a minimax polynomial that isbde2005-11-121-9/+7
| | | | | | | | | | | | 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.
* As for __kernel_cosf() and __kernel_sinf(), use a fairly optimal minimaxbde2005-11-101-17/+11
| | | | | | | | | | | | | | | | | | 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.
* Use a 53-bit approximation to pi/2 instead of a 33+53 bit one for thebde2005-11-061-9/+39
| | | | | | | | | | | | | | | | | | | | | | | | 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.
* Moved the optimization for tiny x from __kernel_tan[f](x) to tan[f](x)bde2005-11-024-46/+14
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Updated the comment about the optimization for tiny x (the previousbde2005-11-021-2/+4
| | | | | | | | | | 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.)
* Removed dead code for handling tan[f]() on odd multiples of pi/2. Thisbde2005-11-022-6/+2
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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 { --- > {
* Fixed some of the silliness related to rev.1.8. In 1.8, "double" inbde2005-11-021-1/+1
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Implement inline functions to give the complex result x+I*y from floatbde2005-10-291-0/+42
| | | | | | | 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>
* Use double precision to simplify and optimize arg reduction for smallbde2005-10-291-97/+22
| | | | | | | | | | | | | | 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.
* Start trying to make the float precision trig functions actually worthbde2005-10-291-30/+28
| | | | | | | | | | | | | | | | | | | | | | | | 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.
* Use fairly optimal minimax polynomials for __kernel_cosf() andbde2005-10-282-16/+15
| | | | | | | | | | | | | | | | | | __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.
OpenPOWER on IntegriCloud