summaryrefslogtreecommitdiffstats
path: root/lib/msun/src
Commit message (Collapse)AuthorAgeFilesLines
* Fix a bug where the wrong argument was passed to INSERT_WORDS().das2011-02-101-2/+2
| | | | | | | This bug results in a type mismatch that happens to be harmless because of the way INSERT_WORDS() works. Submitted by: bde
* For small arguments, these functions use simple approximations,das2011-02-105-5/+5
| | | | | | | e.g. cos(small) = 1, sin(small) = small. This commit tightens the thresholds at which the simple approximations are used. Reviewed by: bde
* Fix a bogus threshold that was copied from the double precision version.das2011-02-101-1/+1
| | | | | | | This commit should have no effect on correctness; it merely changes the threshold at which a simpler approximation can be used. Reviewed by: bde
* Another minor nit: Make sure the constant here is a float so the compilerdas2010-12-071-1/+1
| | | | doesn't promote the entire expression to double.
* Fix various nits in style and comments that were pointed out by bde.das2010-12-073-9/+13
| | | | Code changes verified with md5.
* Add log2() and log2f().das2010-12-054-0/+118
|
* Add a "kernel" log function, based on e_log.c, which is useful fordas2010-12-052-0/+171
| | | | | | | | | | | | | | | implementing accurate logarithms in different bases. This is based on an approach bde coded up years ago. This function should always be inlined; it will be used in only a few places, and rudimentary tests show a 40% performance improvement in implementations of log2() and log10() on amd64. The kernel takes a reduced argument x and returns the same polynomial approximation as e_log.c, but omitting the low-order term. The low-order term is much larger than the rest of the approximation, so the caller of the kernel function can scale it to the appropriate base in extra precision and obtain a much more accurate answer than by using log(x)/log(b).
* Fix bug in jn(3) and jnf(3) that led to -inf resultsuqs2010-11-132-2/+12
| | | | | | | | | | | | | | | | | | | | Explanation by Steve: jn[f](n,x) for certain ranges of x uses downward recursion to compute the value of the function. The recursion sequence that is generated is proportional to the actual desired value, so a normalization step is taken. This normalization is j0[f](x) divided by the zeroth sequence member. As Bruce notes, near the zeros of j0[f](x) the computed value can have giga-ULP inaccuracy. I found for the 1st zero of j0f(x) only the leading decimal digit is correct. The solution to the issue is fairly straight forward. The zeros of j0(x) and j1(x) never coincide, so as j0(x) approaches a zero, the normalization constant switches to j1[f](x) divided by the 2nd sequence member. The expectation is that j1[f](x) is a more accurately computed value. PR: bin/144306 Submitted by: Steven G. Kargl <kargl@troutmask.apl.washington.edu> Reviewed by: bde MFC after: 7 days
* Introduce __isnanf() as an alias for isnanf(), and make the isnan()das2010-06-122-2/+5
| | | | | | | | | | | macro expand to __isnanf() instead of isnanf() for float arguments. This change is needed because isnanf() isn't declared in strict POSIX or C99 mode. Compatibility note: Apps using isnan(float) that are compiled after this change won't link against an older libm. Reported by: Florian Forster <octo@verplant.org>
* Use the documented machine constraint for SSE registers.ed2009-06-111-1/+1
| | | | | | | | The amd64-specific bits of msun use an undocumented constraint, which is less likely to be supported by other compilers (such as Clang). Change the code to use a more common machine constraint. Obtained from: /projects/clangbsd/
* Use ISO C99 style inline semantics in msun.ed2009-06-036-15/+30
| | | | | Because we use ISO C99 nowadays, we can just get rid of enforcing GNU89-style inlining.
* Namespace: scalb() is withdrawn from POSIX.das2009-03-141-1/+4
|
* Eliminate __real__ and __imag__ gccisms.das2009-03-144-15/+42
|
* C99 TC2 now wants FP_FAST_FMA* to be defined to 1, if the macros aredas2009-02-071-3/+3
| | | | defined at all. See also: defect report #223.
* Use __gnu89_inline so that these files will compile with newer versionsdas2009-01-135-5/+5
| | | | | | of gcc, where the meaning of 'inline' was changed to match C99. Noticed by: rdivacky
* Fix the types of INFINITY and NAN, which were broken in r131851. Theydas2009-01-081-2/+2
| | | | | | | | should both be floats, not doubles. PR: 127795 Submitted by: Christoph Mallon MFC after: 2 weeks
* Add support for the FPA floating-point format on ARM. Themarcel2008-12-231-2/+12
| | | | | | | | | FPA floating-point format is identical to the VFP format, but is always stored in big-endian. Introduce _IEEE_WORD_ORDER to describe the byte-order of the FP representation. Obtained from: Juniper Networks, Inc
* Remove some unused variables.das2008-08-084-4/+1
| | | | Reported by: Intel C Compiler
* In the linedas2008-08-083-3/+3
| | | | | | | #pragma STDC CX_LIMITED_RANGE ON the "ON" needs to be in caps. gcc doesn't understand this pragma anyway and assumes it is always on in any case, but icc supports it and cares about the case.
* Implement cproj{,f,l}().das2008-08-073-0/+133
|
* Use cpack() and the gcc extension __imag__ to implement cimag() anddas2008-08-076-6/+18
| | | | | | | | | | | | conj() instead of using expressions like z * I. The latter is bad for several reasons: 1. It is implemented using arithmetic, which is unnecessary, and can generate floating point exceptions, contrary to the requirements on these functions. 2. gcc implements complex multiplication using a formula that breaks down for infinities, e.g., it gives INFINITY * I == nan + inf I.
* 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-022-4/+24
| | | | | | | | | | | | | | | | 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-3110-0/+413
| | | | | | | and cargl(). Reviewed by: bde sparc64 testing resources from: remko
* 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.
* Implement fmodl.das2008-06-192-2/+149
| | | | Document fmodl and fix some errors in the fmod manpage.
* 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
|
* Hook remquol() and remainderl() up to the build.das2008-03-301-1/+1
|
* 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-302-0/+113
|
* Hook hypotl() and cabsl() up to the build.das2008-03-301-2/+0
|
* 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.
* 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.
OpenPOWER on IntegriCloud