summaryrefslogtreecommitdiffstats
path: root/lib/msun
Commit message (Collapse)AuthorAgeFilesLines
* Fix some style bogosity from fdlibm.das2008-08-032-12/+12
|
* Minor improvements:das2008-08-033-20/+20
| | | | | | | - Improve the order of some tests. - Fix style. Submitted by: bde
* A few minor corrections, including some from bde:das2008-08-023-17/+16
| | | | | | | - When y/x is huge, it's faster and more accurate to return pi/2 instead of pi - pi/2. - There's no need for 3 lines of bit fiddling to compute -z. - Fix a comment.
* On i386, gcc truncates long double constants to double precisiondas2008-08-024-8/+127
| | | | | | | | | | | | | | | | at compile time regardless of the dynamic precision, and there's no way to disable this misfeature at compile time. Hence, it's impossible to generate the appropriate tables of constants for the long double inverse trig functions in a straightforward way on i386; this change hacks around the problem by encoding the underlying bits in the table. Note that these functions won't pass the regression test on i386, even with the FPU set to extended precision, because the regression test is similarly damaged by gcc. However, the tests all pass when compiled with a modified version of gcc. Reported by: bde
* Fix some problems with asinf(), acosf(), atanf(), and atan2f():das2008-08-014-71/+42
| | | | | | | | | | | | | | | | | | | | | | | | - Adjust several constants for float precision. Some thresholds that were appropriate for double precision were never changed when these routines were converted to float precision. This has an impact on performance but not accuracy. (Submitted by bde.) - Reduce the degrees of the polynomials used. A smaller degree suffices for float precision. - In asinf(), use double arithmetic in part of the calculation to avoid a corner case and some complicated arithmetic involving a division and some buggy constants. This improves performance and accuracy. Max error (ulps): asinf acosf atanf before 0.925 0.782 0.852 after 0.743 0.804 0.852 As bde points out, it's cheaper for asin*() and acos*() to use polynomials instead of rational functions, but that's a task for another day.
* Add implementations of acosl(), asinl(), atanl(), atan2l(),das2008-07-3120-69/+900
| | | | | | | and cargl(). Reviewed by: bde sparc64 testing resources from: remko
* Set WARNS=1.das2008-07-311-0/+1
| | | | | I believe I've committed all the bits necessary to make this compile on all supported architectures. :crosses fingers:
* The high part of the mantissa is 64 bits on sparc64.das2008-07-312-2/+2
|
* As in other parts of libm, mark a few constants as volatile to preventdas2008-07-316-12/+20
| | | | | | spurious optimizations. gcc doesn't support FENV_ACCESS, so when it folds constants, it assumes that the rounding mode is always the default and floating point exceptions never matter.
* Sort the .PATH entries to give a more reasonable order of precedence:das2008-07-181-4/+5
| | | | | | | | | | | | | 1. architecture-specific files 2. long double format-specific files 3. bsdsrc 4. src 5. man The original order was virtually the opposite of this. This should not cause any functional changes at this time. The difference is only significant when one wants to override, say, a generic foo.c with a more specialized foo.c (as opposed to foo.S).
* Fix a typo in the cosl() prototype.das2008-06-281-2/+2
|
* Implement fmodl.das2008-06-195-25/+177
| | | | Document fmodl and fix some errors in the fmod manpage.
* Symbol.map is handled by cpp, so use C-style commentsgonzo2008-05-031-1/+3
| | | | Approved by: cognet (mentor)
* Add mips support to libm, from mips2-jnpr perforce branch.imp2008-04-264-0/+259
|
* Fix some corner cases:das2008-04-032-10/+20
| | | | | | | | | - fma(x, y, z) returns z, not NaN, if z is infinite, x and y are finite, x*y overflows, and x*y and z have opposite signs. - fma(x, y, z) doesn't generate an overflow, underflow, or inexact exception if z is NaN or infinite, as per IEEE 754R. - If the rounding mode is set to FE_DOWNWARD, fma(1.0, 0.0, -0.0) is -0.0, not +0.0.
* Remove a (bogus) remnant of debugging this on sparc64.das2008-03-311-1/+1
|
* Add assembly versions of remquol() and remainderl().das2008-03-306-4/+228
|
* Hook remquol() and remainderl() up to the build.das2008-03-304-12/+33
|
* Implement remainderl() as a wrapper around remquol(). The extra workdas2008-03-302-0/+44
| | | | remquol() performs to compute the quotient is negligible.
* Implement remquol() based on remquo().das2008-03-302-0/+183
|
* Implement csqrtl().das2008-03-305-13/+135
|
* Hook hypotl() and cabsl() up to the build.das2008-03-303-6/+7
|
* Document hypotl().das2008-03-301-12/+17
| | | | Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
* Alias hypotl() and cabsl() for platforms where long double is the samedas2008-03-302-2/+12
| | | | as double.
* Implement cabsl() in terms of hypotl().das2008-03-301-0/+20
| | | | Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
* Implement hypotl(). This is bde's conversion of fdlibm hypot(), with minordas2008-03-301-0/+139
| | | | fixes for ld128 by me.
* Use fabs[f]() instead of bit fiddling for setting absolute values.bde2008-03-302-4/+4
| | | | | | | | | | | This makes little difference in float precision, but in double precision gives a speedup of about 30% on amd64 (A64 CPU) and i386 (A64). This depends on fabs[f]() being inline and efficient. The bit fiddling (or any use of SET_HIGH_WORD(), which libm does too much because it was best on old 32-bit machines) always causes packing overheads and sometimes causes stalls in the packing, since it operates on only part of a variable in the double precision case. It apparently did cause stalls in a critical path here.
* Use the expression fabs(x+0.0)-fabs(y+0.0) instead ofbde2008-03-302-2/+2
| | | | | | | | fabs(x+0.0)+fabs(y+0.0) when mixing NaNs. This improves consistency of the result by making it harder for the compiler to reorder the operands. (FP addition is not necessarily commutative because the order of operands makes a difference on some machines iff the operands are both NaNs.)
* Fix a missing mask in a hi+lo decomposition. Thus bug made the extrabde2008-03-301-1/+1
| | | | | | precision in software useless, so hypotf() had some errors in the 1-2 ulp range unless there is extra precision in hardware (as happens on i386).
* Include math.h for the fmaf() prototype.das2008-03-291-0/+2
|
* Fix some rather obscene code that has ambiguous if...if...else...das2008-03-296-10/+36
| | | | constructs in it.
* 1 << 47 needs to be written 1ULL << 47.das2008-03-021-1/+1
|
* Hook up sqrtl() to the build.das2008-03-025-13/+24
|
* MD implementations of sqrtl().das2008-03-022-0/+76
|
* MI implementation of sqrtl(). This is very slow and shoulddas2008-03-023-4/+168
| | | | be overridden when hardware sqrt is available.
* Fix and improve some magic numbers for the "medium size" case.bde2008-02-282-4/+4
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Inline __ieee754__rem_pio2f(). On amd64 (A64) and i386 (A64), thisbde2008-02-254-0/+15
| | | | | | | | | | | | | 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.
* Use a temporary array instead of the arg array y[] for callingbde2008-02-252-8/+8
| | | | | | | | | | | | | | __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.
* Change __ieee754_rem_pio2f() to return double instead of float so thatbde2008-02-255-27/+24
| | | | | | | | | | | | | | | | | | | | | | 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.
* Fix some off-by-1 errors.bde2008-02-253-6/+11
| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | 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.
* Let PowerPC world optionally build with -msoft-float. For FPU-less PowerPCraj2008-02-241-0/+5
| | | | | | | | | variations (e500 currently), this provides a gcc-level FPU emulation and is an alternative approach to the recently introduced kernel-level emulation (FPU_EMU). Approved by: cognet (mentor) MFp4: e500
* Optimize the 9pi/2 < |x| <= 2**19pi/2 case some more by avoiding anbde2008-02-232-11/+8
| | | | | | | | | | | | | | | | | | | | | | | | | 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
* Avoid using FP-to-integer conversion for !(amd64 || i386) too. Use thebde2008-02-222-6/+4
| | | | | | | | | | | | | | | | | | | 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).
* Remove the "quick check no cancellation" optimization forbde2008-02-221-12/+1
| | | | | | | | | | | | | | | | | | 9pi/2 < |x| < 32pi/2 since it is only a small or negative optimation and it gets in the way of further optimizations. It did one more branch to avoid some integer operations and to use a different dependency on previous results. The branches are fairly predictable so they are usually not a problem, so whether this is a good optimization depends mainly on the timing for the previous results, which is very machine-dependent. On amd64 (A64), this "optimization" is a pessimization of about 1 cycle or 1%; on ia64, it is an optimization of about 2 cycles or 1%; on i386 (A64), it is an optimization of about 5 cycles or 4%; on i386 (Celeron P2) it is an optimization of about 4 cycles or 3% for cos but a pessimization of about 5 cycles for sin and 1 cycle for tan. I think the new i386 (A64) slowness is due to an pipeline stall due to an avoidable load-store mismatch (so the old timing was better), and the i386 (Celeron) variance is due to its branch predictor not being too good.
* Optimize the 9pi/2 < |x| <= 2**19pi/2 case on amd64 and i386 by avoidingbde2008-02-222-0/+18
| | | | | | | | | | | | | | | | | | | | | | 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.
* Add an irint() function in inline asm for amd64 and i386. irint() isbde2008-02-221-0/+30
| | | | | | | | | | | | | the same as lrint() except it returns int instead of long. Though the extern lrint() is fairly fast on these arches, it still takes about 12 cycles longer than the inline version, and 12 cycles is a lot in applications where [li]rint() is used to avoid slow conversions that are only a couple of times slower. This is only for internal use. The libm versions of *rint*() should also be inline, but that would take would take more header engineering. Implementing irint() instead of lrint() also avoids a conflict with the extern declaration of the latter.
* Optimize the conversion to bits a little (by about 11 cycles or 16%bde2008-02-221-5/+13
| | | | | | | | | | | | on i386 (A64), 5 cycles on amd64 (A64), and 3 cycles on ia64). gcc tends to generate very bad code for accessing floating point values as bits except when the integer accesses have the same width as the floating point values, and direct accesses to bit-fields (as is common only for long double precision) always gives such accesses. Use the expsign access method, which is good for 80-bit long doubles and hopefully no worse for 128-bit long doubles. Now the generated code is less bad. There is still unnecessary copying of the arg on amd64 and i386 and mysterious extra slowness on amd64.
* Optimize the fixup for +-0 by using better classification for this casebde2008-02-221-2/+4
| | | | | and by using a table lookup to avoid a branch when this case occurs. On i386, this saves 1-4 cycles out of about 64 for non-large args.
* Fix rintl() on signaling NaNs and unsupported formats.bde2008-02-221-5/+3
|
* s/rcsid/__FBSDID/das2008-02-2279-241/+158
|
OpenPOWER on IntegriCloud