#
285830 |
|
23-Jul-2015 |
gjb |
- Copy stable/10@285827 to releng/10.2 in preparation for 10.2-RC1 builds. - Update newvers.sh to reflect RC1. - Update __FreeBSD_version to reflect 10.2. - Update default pkg(8) configuration to use the quarterly branch.[1]
Discussed with: re, portmgr [1] Approved by: re (implicit) Sponsored by: The FreeBSD Foundation |
#
256281 |
|
10-Oct-2013 |
gjb |
Copy head (r256279) to stable/10 as part of the 10.0-RELEASE cycle.
Approved by: re (implicit) Sponsored by: The FreeBSD Foundation
|
#
239195 |
|
11-Aug-2012 |
dim |
Add __always_inline to __ieee754_rem_pio2() and __ieee754_rem_pio2f(), since some older versions of gcc refuse to inline these otherwise.
Requested by: bde MFC after: 1 week
|
#
239192 |
|
11-Aug-2012 |
dim |
Change a few extern inline functions in libm to static inline, since they need to refer to static constants, which C99 does not allow for extern inline functions.
While here, change a comment in e_rem_pio2f.c to mention the correct number of bits.
Reviewed by: bde MFC after: 1 week
|
#
193368 |
|
03-Jun-2009 |
ed |
Use ISO C99 style inline semantics in msun.
Because we use ISO C99 nowadays, we can just get rid of enforcing GNU89-style inlining.
|
#
187128 |
|
13-Jan-2009 |
das |
Use __gnu89_inline so that these files will compile with newer versions of gcc, where the meaning of 'inline' was changed to match C99.
Noticed by: rdivacky
|
#
181405 |
|
08-Aug-2008 |
das |
Remove some unused variables.
Reported by: Intel C Compiler
|
#
176640 |
|
28-Feb-2008 |
bde |
Fix and improve some magic numbers for the "medium size" case.
e_rem_pio2.c: This case goes up to about 2**20pi/2, but the comment about it said that it goes up to about 2**19pi/2.
It went too far above 2**pi/2, giving a multiplier fn with 21 significant bits in some cases. This would be harmful except for a numerical accident. It happens that the terms of the approximation to pi/2, when rounded to 33 bits so that multiplications by 20-bit fn's are exact, happen to be rounded to 32 bits so multiplications by 21-bit fn's are exact too, so the bug only complicates the error analysis (we might lose a bit of accuracy but have bits to spare).
e_rem_pio2f.c: The bogus comment in e_rem_pio2.c was copied and the code was changed to be bug-for-bug compatible with it, except the limit was made 90 ulps smaller than necessary. The approximation to pi/2 was not modified except for discarding some of it.
The same rough error analysis that justifies the limit of 2**20pi/2 for double precision only justifies a limit of 2**18pi/2 for float precision. We depended on exhaustive testing to check the magic numbers for float precision. More exaustive testing shows that we can go up to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision, with a the maximum error for cosf() and sinf() unchanged at 0.5009 ulps despite the maximum error in rem_pio2f being ~0.25 ulps. Implement this.
|
#
176569 |
|
25-Feb-2008 |
bde |
Inline __ieee754__rem_pio2f(). On amd64 (A64) and i386 (A64), this gives an average speedup of about 12 cycles or 17% for 9pi/4 < |x| <= 2**19pi/2 and a smaller speedup for larger x, and a small speeddown for |x| <= 9pi/4 (only 1-2 cycles average, but that is 4%).
Inlining this is less likely to bust caches than inlining the float version since it is much smaller (about 220 bytes text and rodata) and has many fewer branches. However, the float version was already large due to its manual inlining of the branches and also the polynomial evaluations.
|
#
176558 |
|
25-Feb-2008 |
bde |
Use a temporary array instead of the arg array y[] for calling __kernel_rem_pio2(). This simplifies analysis of aliasing and thus results in better code for the usual case where __kernel_rem_pio2() is not called. In particular, when __ieee854_rem_pio2[f]() is inlined, it normally results in y[] being returned in registers. I couldn't get this to work using the restrict qualifier.
In float precision, this saves 2-3% in most cases on amd64 and i386 (A64) despite it not being inlined in float precision yet. In double precision, this has high variance, with an average gain of 2% for amd64 and 0.7% for i386 (but a much larger gain for usual cases) and some losses.
|
#
176552 |
|
25-Feb-2008 |
bde |
Change __ieee754_rem_pio2f() to return double instead of float so that this function and its callers cosf(), sinf() and tanf() don't waste time converting values from doubles to floats and back for |x| > 9pi/4. All these functions were optimized a few years ago to mostly use doubles internally and across the __kernel*() interfaces but not across the __ieee754_rem_pio2f() interface.
This saves about 40 cycles in cosf(), sinf() and tanf() for |x| > 9pi/4 on amd64 (A64), and about 20 cycles on i386 (A64) (except for cosf() and sinf() in the upper range). 40 cycles is about 35% for |x| < 9pi/4 <= 2**19pi/2 and about 5% for |x| > 2**19pi/2. The saving is much larger on amd64 than on i386 since the conversions are not easy to optimize except on i386 where some of them are automatic and others are optimized invalidly. amd64 is still about 10% slower in cosf() and tanf() in the lower range due to conversion overhead.
This also gives a tiny speedup for |x| <= 9pi/4 on amd64 (by simplifying the code). It also avoids compiler bugs and/or additional slowness in the conversions on (not yet supported) machines where double_t != double.
|
#
176550 |
|
25-Feb-2008 |
bde |
Fix some off-by-1 errors.
e_rem_pio2.c: Float and double precision didn't work because init_jk[] was 1 too small. It needs to be 2 larger than you might expect, and 1 larger than it was for these precisions, since its test for recomputing needs a margin of 47 bits (almost 2 24-bit units).
init_jk[] seems to be barely enough for extended and quad precisions. This hasn't been completely verified. Callers now get about 24 bits of extra precision for float, and about 19 for double, but only about 8 for extended and quad. 8 is not enough for callers that want to produce extra-precision results, but current callers have rounding errors of at least 0.8 ulps, so another 1/2**8 ulps of error from the reduction won't affect them much.
Add a comment about some of the magic for init_jk[].
e_rem_pio2.c: Double precision worked in practice because of a compensating off-by-1 error here. Extended precision was asked for, and it executed exactly the same code as the unbroken double precision.
e_rem_pio2f.c: Float precision worked in practice because of a compensating off-by-1 error here. Double precision was asked for, and was almost needed, since the cosf() and sinf() callers want to produce extra-precision results, at least internally so that their error is only 0.5009 ulps. However, the extra precision provided by unbroken float precision is enough, and the double-precision code has extra overheads, so the off-by-1 error cost about 5% in efficiency on amd64 and i386.
|
#
176476 |
|
23-Feb-2008 |
bde |
Optimize the 9pi/2 < |x| <= 2**19pi/2 case some more by avoiding an fabs(), a conditional branch, and sign adjustments of 3 variables for x < 0 when the branch is taken. In double precision, even when the branch is perfectly predicted, this saves about 10 cycles or 10% on amd64 (A64) and i386 (A64) for the negative half of the range, but makes little difference for the positive half of the range. In float precision, it also saves about 4 cycles for the positive half of the range on i386, and many more cycles in both halves on amd64 (28 in the negative half and 11 in the positive half for tanf), but the amd64 times for float precision are anomalously slow so the larger improvement is only a side effect.
Previous commits arranged for the x < 0 case to be handled simply: - one part of the rounding method uses the magic number 0x1.8p52 instead of the usual 0x1.0p52. The latter is required for large |x|, but it doesn't work for negative x and we don't need it for large |x|. - another part of the rounding method no longer needs to add `half'. It would have needed to add -half for negative x. - removing the "quick check no cancellation" in the double precision case removed the need to take the absolute value of the quadrant number.
Add my noncopyright in e_rem_pio2.c
|
#
176467 |
|
22-Feb-2008 |
bde |
Avoid using FP-to-integer conversion for !(amd64 || i386) too. Use the FP-to-FP method to round to an integer on all arches, and convert this to an int using FP-to-integer conversion iff irint() is not available. This is cleaner and works well on at least ia64, where it saves 20-30 cycles or about 10% on average for 9Pi/4 < |x| <= 32pi/2 (should be similar up to 2**19pi/2, but I only tested the smaller range).
After the previous commit to e_rem_pio2.c removed the "quick check no cancellation" non-optimization, the result of the FP-to-integer conversion is not needed so early, so using irint() became a much smaller optimization than when it was committed.
An earlier commit message said that cos, cosf, sin and sinf were equally fast on amd64 and i386 except for cos and sin on i386. Actually, cos and sin on amd64 are equally fast to cosf and sinf on i386 (~88 cycles), while cosf and sinf on amd64 are not quite equally slow to cos and sin on i386 (average 115 cycles with more variance).
|
#
176465 |
|
22-Feb-2008 |
bde |
Optimize the 9pi/2 < |x| <= 2**19pi/2 case on amd64 and i386 by avoiding the the double to int conversion operation which is very slow on these arches. Assume that the current rounding mode is the default of round-to-nearest and use rounding operations in this mode instead of faking this mode using the round-towards-zero mode for conversion to int. Round the double to an integer as a double first and as an int second since the double result is needed much earler.
Double rounding isn't a problem since we only need a rough approximation. We didn't support other current rounding modes and produce much larger errors than before if called in a non-default mode.
This saves an average about 10 cycles on amd64 (A64) and about 25 on i386 (A64) for x in the above range. In some cases the saving is over 25%. Most cases with |x| < 1000pi now take about 88 cycles for cos and sin (with certain CFLAGS, etc.), except on i386 where cos and sin (but not cosf and sinf) are much slower at 111 and 121 cycles respectivly due to the compiler only optimizing well for float precision. A64 hardware cos and sin are slower at 105 cycles on i386 and 110 cycles on amd64.
|
#
176410 |
|
19-Feb-2008 |
bde |
Merge cosmetic changes from e_rem_pio2.c 1.10 (convert to __FBSDID(); fix indentation and return type of __ieee754_rem_pio2()).
Remove unused variables.
|
#
176356 |
|
17-Feb-2008 |
das |
Add more pi for long doubles. Also, avoid storing multiple copies of the pi/2 array, as it is unlikely to vary, except in Indiana.
|
#
152707 |
|
23-Nov-2005 |
bde |
Simplified setiing up args for __kernel_rem_pio2(). We already have x with a 24-bit fraction, so we don't need a loop to split it into up to 3 terms with 24-bit fractions.
|
#
152706 |
|
23-Nov-2005 |
bde |
Quick fix for stack buffer overrun in rev.1.13. Oops. The prec == 1 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.
|
#
152596 |
|
19-Nov-2005 |
bde |
Moved all the optimizations for |x| <= 9pi/2 from __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.
|
#
152535 |
|
17-Nov-2005 |
bde |
Rearranged the the optimizations for special cases to reduce the average 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.
|
#
152132 |
|
06-Nov-2005 |
bde |
Use a 53-bit approximation to pi/2 instead of a 33+53 bit one for the special case pi/4 <= |x| < 3*pi/4. This gives a tiny optimization (it saves 2 subtractions, which are scheduled well so they take a whole 1 cycle extra on an AthlonXP), and simplifies the code so that the following optimization is not so ugly.
Optimize for the range 3*pi/4 < |x| < 9*Pi/2 in the same way. On Athlon{XP,64} systems, this gives a 25-40% optimization (depending a lot on CFLAGS) for the cosf() and sinf() consumers on this range. Relative to i387 hardware fcos and fsin, it makes the software versions faster in most cases instead of slower in most cases. The relative optimization is smaller for tanf() the inefficient part is elsewhere.
The 53-bit approximation to pi/2 is good enough for pi/4 <= |x| < 3*pi/4 because after losing up to 24 bits to subtraction, we still have 29 bits of precision and only need 25 bits. Even with only 5 extra bits, it is possible to get perfectly rounded results starting with the reduced x, since if x is nearly a multiple of pi/2 then x is not near a half-way case and if x is not nearly a multiple of pi/2 then we don't lose many bits. With our intentionally imperfect rounding we get the same results for cosf(), sinf() and tanf() as without this optimization.
|
#
151864 |
|
29-Oct-2005 |
bde |
Use double precision to simplify and optimize arg reduction for small and medium size args too: instead of conditionally subtracting a float 17+24, 17+17+24 or 17+17+17+24 bit approximation to pi/2, always subtract a double 33+53 bit one. The float version is now closer to the double version than to old versions of itself -- it uses the same 33+53 bit approximation as the simplest cases in the double version, and where the float version had to switch to the slow general case at |x| == 2^7*pi/2, it now switches at |x| == 2^19*pi/2 the same as the double version.
This speeds up arg reduction by a factor of 2 for |x| between 3*pi/4 and 2^7*pi/4, and by a factor of 7 for |x| between 2^7*pi/4 and 2^19*pi/4.
|
#
151855 |
|
29-Oct-2005 |
bde |
Start trying to make the float precision trig functions actually worth using under FreeBSD. Before this commit, all float precision functions except exp2f() were implemented using only float precision, apparently because Cygnus needed this in 1993 for embedded systems with slow or inefficient double precision. For FreeBSD, except possibly on systems that do floating point entirely in software (very old i386 and now arm), this just gives a more complicated implementation, many bugs, and usually worse performance for float precision than for double precision. The bugs and worse performance were particulary large in arg reduction for trig functions. We want to divide by an approximation to pi/2 which has as many as 1584 bits, so we should use the widest type that is efficient and/or easy to use, i.e., double. Use fdlibm's __kernel_rem_pio2() to do this as Sun apparently intended. Cygnus's k_rem_pio2f.c is now unused. e_rem_pio2f.c still needs to be separate from e_rem_pio2.c so that it can be optimized for float args. Similarly for long double precision.
This speeds up cosf(x) on large args by a factor of about 2. Correct arg reduction on large args is still inherently very slow, so hopefully these args rarely occur in practice. There is much more efficiency to be gained by using double precision to speed up arg reduction on medium and small float args.
|
#
151648 |
|
25-Oct-2005 |
bde |
More fixes for arg reduction near pi/2 on systems with broken assignment to floats (mainly i386's). All errors of more than 1 ulp for float precision trig functions were supposed to have been fixed; however, compiling with gcc -O2 uncovered 18250 more such errors for cosf(), with a maximum error of 1.409 ulps.
Use essentially the same fix as in rev.1.8 of k_rem_pio2f.c (access a non-volatile variable as a volatile). Here the -O1 case apparently worked because the variable is in a 2-element array and it takes -O2 to mess up such a variable by putting it in a register.
The maximum error for cosf() on i386 with gcc -O2 is now 0.5467 (it is still 0.5650 with gcc -O1). This shows that -O2 still causes some extra precision, but the extra precision is now good.
Extra precision is harmful mainly for implementing extra precision in software. We want to represent x+y as w+r where both "+" operations are in infinite precision and r is tiny compared with w. There is a standard algorithm for this (Knuth (1981) 4.2.2 Theorem C), and fdlibm uses this routinely, but the algorithm requires w and r to have the same precision as x and y. w is just x+y (calculated in the same finite precision as x and y), and r is a tiny correction term. The i386 gcc bugs tend to give extra precision in w, and then using this extra precision in the calculation of r results in the correction mostly staying in w and being missing from r. There still tends to be no problem if the result is a simple expression involving w and r -- modulo spills, w keeps its extra precision and r remains the right correction for this wrong w. However, here we want to pass w and r to extern functions. Extra precision is not retained in function args, so w gets fixed up, but the change to the tiny r is tinier, so r almost remains as a wrong correction for the right w.
|
#
151221 |
|
10-Oct-2005 |
bde |
Fixed range reduction near (but not very near) medium-sized multiples of pi/2 (1 line) and expand a comment about related magic (many lines).
The bug was essentially the same as for the +-pi/2 case (a mistranslated mask), but was smaller so it only significantly affected multiples starting near +-13*pi/2. At least on amd64, for cosf() on all 2^32 float args, the bug caused 128 errors of >= 1 ulp, with a maximum error of 1.2393 ulps.
|
#
151147 |
|
09-Oct-2005 |
bde |
Oops, the last-minute optimization in rev.1.8 wasn't a good idea. The 17+17+24 bit pi/2 must only be used when subtraction of the first 2 terms in it from the arg is exact. This happens iff the the arg in bits is one of the 2**17[-1] values on each side of (float)(pi/2).
Revert to the algorithm in rev.1.7 and only fix its threshold for using the 3-term pi/2. Use the threshold that maximizes the number of values for which the 3-term pi/2 is used, subject to not changing the algorithm for comparing with the threshold. The 3-term pi/2 ends up being used for about half of its usable range (about 64K values on each side).
|
#
151112 |
|
08-Oct-2005 |
bde |
Fixed syntax error (a missing brace) in previous commit.
|
#
151110 |
|
08-Oct-2005 |
bde |
Fixed range reduction near (but not very near) +-pi/2. A bug caused a maximum error of 2.905 ulps for cosf(), but the algorithm for cosf() is good for < 1 ulps and happens to give perfect rounding (< 0.5 ulps) near +-pi/2 except for the bug. The extra relative errors for tanf() were similar (slightly larger). The bug didn't affect sinf() since sinf'(+-pi/2) is 0.
For range reduction in ~[-3pi/4, -pi/4] and ~[pi/4, 3pi/4] we must subtract +-pi/2 and the only complication is that this must be done in extra precision. We have handy 17+24-bit and 17+17+24-bit approximations to pi/2. If we always used the former then we would lose up to 24 bits of accuracy due to cancelation of leading bits, but we need to keep at least 24 bits plus a guard digit or 2, and should keep as many guard bits as efficiency permits. So we used the less-precise pi/2 not very near +-pi/2 and switched to using the more-precise pi/2 very near +-pi/2. However, we got the threshold for the switch wrong by allowing 19 bits to cancel, so we ended up with only 21 or 22 bits of accuracy in some cases, which is even worse than naively subtracting pi/2 would have done.
Exhaustive checking shows that allowing only 17 bits to cancel (min. accuracy ~24 bits) is sufficient to reduce the maximum error for cosf() near +-pi/2 to 0.726 ulps, but allowing only 6 bits to cancel (min. accuracy ~35-bits) happens to give perfect rounding for cosf() at little extra cost so we prefer that.
We actually (in effect) allow 0 bits to cancel and always use the 17+17+24-bit pi/2 (min. accuracy ~41 bits). This is simpler and probably always more efficient too. Classifying args to avoid using this pi/2 when it is not needed takes several extra integer operations and a branch, but just using it takes only 1 FP operation.
The patch also fixes misspelling of 17 as 24 in many comments.
For the double-precision version, the magic numbers include 33+53 bits for the less-precise pi/2 and (53-32-1 = 20) bits being allowed to cancel, so there are ~33-20 = 13 guard bits. This is sufficient except probably for perfect rounding. The more-precise pi/2 has 33+33+53 bits and we still waste time classifying args to avoid using it.
The bug is apparently from mistranslation of the magic 32 in 53-32-1. The number of bits allowed to cancel is not critical and we use 32 for double precision because it allows efficient classification using a 32-bit comparison. For float precision, we must use an explicit mask, and there are fewer bits so there is less margin for error in their allocation. The 32 got reduced to 4 but should have been reduced almost in proportion to the reduction of mantissa bits.
|
#
97409 |
|
28-May-2002 |
alfred |
Assume __STDC__, remove non-__STDC__ code.
Reviewed by: md5
|
#
50476 |
|
28-Aug-1999 |
peter |
$Id$ -> $FreeBSD$
|
#
22993 |
|
22-Feb-1997 |
peter |
Revert $FreeBSD$ to $Id$
|
#
21673 |
|
14-Jan-1997 |
jkh |
Make the long-awaited change from $Id$ to $FreeBSD$
This will make a number of things easier in the future, as well as (finally!) avoiding the Id-smashing problem which has plagued developers for so long.
Boy, I'm glad we're not using sup anymore. This update would have been insane otherwise.
|
#
8870 |
|
30-May-1995 |
rgrimes |
Remove trailing whitespace.
|
#
7659 |
|
07-Apr-1995 |
bde |
Submitted by: J.T. Conklin <jtc@wimsey.com>
Second part of update to fdlibm 5.2: speed up argument reduction for trig functions in the case pi/4 < |x| < 3pi/4.
Remove unused static constants ("one").
|
#
2117 |
|
19-Aug-1994 |
jkh |
This commit was generated by cvs2svn to compensate for changes in r2116, which included commits to RCS files with non-trunk default branches.
|
#
2116 |
|
19-Aug-1994 |
jkh |
J.T. Conklin's latest version of the Sun math library.
-- Begin comments from J.T. Conklin: The most significant improvement is the addition of "float" versions of the math functions that take float arguments, return floats, and do all operations in floating point. This doesn't help (performance) much on the i386, but they are still nice to have.
The float versions were orginally done by Cygnus' Ian Taylor when fdlibm was integrated into the libm we support for embedded systems. I gave Ian a copy of my libm as a starting point since I had already fixed a lot of bugs & problems in Sun's original code. After he was done, I cleaned it up a bit and integrated the changes back into my libm. -- End comments
Reviewed by: jkh Submitted by: jtc
|