355395 |
04-Dec-2019 |
dim |
r355120 | dim | 2019-11-26 23:01:09 +0100 (Tue, 26 Nov 2019) | 32 lines
The fdlibm hypot() implementations shouldn't potentially left-shift negative numbers (invoking undefined behavior)
Summary: Various paths through hypot(x, y) will multiply x and y by a power of two, perform the calculation in a range where IEEE-754 provides greater precision, then undo the multiplication to determine the true result. Undoing that multiplication is implemented as t1*w, where t1=2**k.
2**k is often computed by taking the high word of 1.0, then adding k<<20 (for doubles or long doubles) or k<<23 (for floats) to it, then overwriting that high word. But when k is negative this left-shifts a negative value -- and that's undefined behavior in many editions of C and C++.
This patch should fix all hypot implementations to compute 2**k without triggering this particular bit of undefined behavior.
Test Plan: I've only very lightly tested out the hypot(double, double) change, in SpiderMonkey's JavaScript engine, for consistency with prior behavior. The other functions' changes have more or less only been eyeballed. Careful examination appreciated! Do note, however, that an error in any of these changes would most likely produce a value that is incorrect by a factor of two, so any mistake would most likely be glaring if invoked.
Submitted by: Jeff Walden <jwalden@mit.edu> Obtained from: https://github.com/freebsd/freebsd/pull/414 Reviewed by: dim, lwhsu Differential Revision: https://reviews.freebsd.org/D22354 |
354596 |
10-Nov-2019 |
dim |
MFC r354255:
Add __isnan()/__isnanf() aliases for compatibility with glibc and CUDA
Even though clang comes with a number of internal CUDA wrapper headers, compiling sample CUDA programs will result in errors similar to:
In file included from <built-in>:1: In file included from /usr/lib/clang/9.0.0/include/__clang_cuda_runtime_wrapper.h:204: /usr/home/arr/cuda/var/cuda-repo-10-0-local-10.0.130-410.48/usr/local/cuda-10.0//include/crt/math_functions.hpp:2910:7: error: no matching function for call to '__isnan' if (__isnan(a)) { ^~~~~~~ /usr/lib/clang/9.0.0/include/__clang_cuda_device_functions.h:460:16: note: candidate function not viable: call to __device__ function from __host__ function __DEVICE__ int __isnan(double __a) { return __nv_isnand(__a); } ^
CUDA expects __isnan() and __isnanf() declarations to be available, which are glibc specific extensions, equivalent to the regular isnan() and isnanf().
To provide these, define __isnan() and __isnanf() as aliases of the already existing static inline functions __inline_isnan() and __inline_isnanf() from math.h.
Reported by: arrowd PR: 241550 |
352835 |
28-Sep-2019 |
dim |
MFC r352710:
Do not left-shift a negative number (inducing undefined behavior in C/C++) in exp(3), expf(3), expm1(3) and expm1f(3) during intermediate computations that compute the IEEE-754 bit pattern for |2**k| for integer |k|.
The implementations of exp(3), expf(3), expm1(3) and expm1f(3) need to compute IEEE-754 bit patterns for 2**k in certain places. (k is an integer and 2**k is exactly representable in IEEE-754.)
Currently they do things like 0x3FF0'0000+(k<<20), which is to say they take the bit pattern representing 1 and then add directly to the exponent field to get the desired power of two. This is fine when k is non-negative.
But when k<0 (and certain classes of input trigger this), this left-shifts a negative number -- an operation with undefined behavior in C and C++.
The desired semantics can be achieved by instead adding the possibly-negative k to the IEEE-754 exponent bias to get the desired exponent field, _then_ shifting that into its proper overall position.
(Note that in case of s_expm1.c and s_expm1f.c, there are SET_HIGH_WORD and SET_FLOAT_WORD uses further down in each of these files that perform shift operations involving k, but by these points k's range has been restricted to 2 < k <= 56, and the shift operations under those circumstances can't do anything that would be UB.)
Submitted by: Jeff Walden, https://github.com/jswalden Obtained from: https://github.com/freebsd/freebsd/pull/411 Obtained from: https://github.com/freebsd/freebsd/pull/412 |
343258 |
21-Jan-2019 |
pfg |
MFC r343023: msun: reduce diff between src/e_j0.c and src/e_j0f.c
PR: 229501 |
326830 |
13-Dec-2017 |
dim |
MFC r326748:
Document the existence and precision of the remaining long double functions for which an imprecise stub implementation was added in r255294, namely powl(3) and tgammal(3).
Submitted by: Steve Kargl
MFC r326753:
Correct r326748, indicating that tgammal(3) is mapped to tgamma(3), not to itself.
Noticed by: jilles |
326747 |
10-Dec-2017 |
dim |
MFC r326670:
Remove an unused incude from lib/msun/bsdsrc/b_log.c.
Submitted by: Steve Kargl |
326746 |
10-Dec-2017 |
dim |
MFC r326669:
Remove the sentence in math(3) about some long double math functions not being available.
Submitted by: Steve Kargl |
321421 |
24-Jul-2017 |
ngie |
MFC r290605,r290606:
r290605:
Document powl(3)
PR: 191751
r290606:
Bump .Dd |
319379 |
01-Jun-2017 |
ngie |
MFC r319157:
fma_test: mute a warning about unreachable code on amd64 by restructuring the #ifdef block to only handle the rest of the logic in the loop in the #else case.
CID: 1346844 |
319376 |
01-Jun-2017 |
ngie |
MFC r319158:
logarithm_test: assert that feclearexcept succeeds
This helps ensure that test preconditons are fulfilled.
CID: 1346572 |
319252 |
30-May-2017 |
asomers |
MFC r315032-r315036, r315039, r315041
r315032: Increase WARNS for iconv tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315033: Increase WARNS for nss tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315034: Document that the msun tests require WARNS=0
ATF tests have a default WARNS of 0, unlike other usermode programs. This change is technically a noop, but it documents that the msun tests don't work with any warnings enabled, at least not on all architectures.
Reviewed by: ngie Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315035: Increase WARNS for libcrypt tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie, julian Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315036: Increase WARNS for libmp tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie, julian Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315039: Increase WARNS for libutil tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie, julian Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933
r315041: Increase WARNS for pw tests
ATF tests have a default WARNS of 0, unlike other usermode programs.
Reviewed by: ngie, julian Sponsored by: Spectra Logic Corporation Differential Revision: https://reviews.freebsd.org/D9933 |
318365 |
16-May-2017 |
dim |
MFC r318259:
Silence a -Wunused warning about the junk variable being used to raise an inexact floating point exception. The variable cannot be eliminated, unfortunately, otherwise the desired addition triggering the exception will be emitted neither by clang, nor by gcc.
Reviewed by: Steve Kargl, bde |
314023 |
21-Feb-2017 |
ngie |
MFC r312213,r313713:
r312213:
Turn COMPILER_VERSION/COMPILER_TYPE make check into a compile-time check of the clang version
This works around breakage on ^/stable/10 when running installworld from a ^/stable/10 host where the test wouldn't be compiled on the first go-around and would be missing when make installworld is run.
PR: 208703
r313713:
Handle clang 4.x+ with the compile-time exception added in r312213
It also fails the assertions noted in bug 208703
PR: 208703 PR: 217084 |
313535 |
10-Feb-2017 |
ngie |
MFC r311925,r311968,r311969,r312008:
r311925:
Import testcase updates with code contributed back to NetBSD
This also (inadvertently) contains an update to contrib/netbsd-tests/lib/libc/sys/t_wait.c (new testcases).
In collaboration with: christos@NetBSD.org
r311968:
Fix lib/libc/sys/access_test after r311925
sys/param.h needs to be #included in order for __FreeBSD_version to be checked
r311969:
Remove __HAVE_LONG_DOUBLE #define from t_strtod.c and place it in Makefile
This is to enable support in other testcases
Inspired by lib/msun/tests/Makefile .
r312008:
Upgrade NetBSD tests to 01.11.2017_23.20 snapshot
This contains some new testcases in /usr/tests/...:
- .../lib/libc - .../lib/libthr - .../lib/msun - .../sys/kern
Tested on: amd64, i386 |
313498 |
10-Feb-2017 |
ngie |
MFC r305358,r305449,r305451,r306367,r306397,r309474:
This also contains a merge of ^/projects/netbsd-tests-update-12@r304035 . This change never hit ^/head because bin/cat's behavior was changed (on ^/head) to match NetBSD.
PR: 210607
r305358:
Update contrib/netbsd-tests with new content from NetBSD
This updates the snapshot from 09/30/2014 to 08/11/2016
This brings in a number of new testcases from upstream, most notably:
- bin/cat - lib/libc - lib/msun - lib/libthr - usr.bin/sort
lib/libc/tests/stdio/open_memstream_test.c was moved to lib/libc/tests/stdio/open_memstream2_test.c to accomodate the new open_memstream test from NetBSD.
Tested on: amd64 (VMware fusion VM; various bare metal platforms); i386 (VMware fusion VM); make tinderbox
r305449:
Install h_db to unbreak some of the lib/libc/db testcases after r305358
r305451:
Fix lib/libc/rpc test assumptions added in r305358
- Require root in the tcp/udp subtests (it's needed on FreeBSD when registering services). - Skip the tests if service registration fails.
r306367 (by br):
Allow up to 6 arguments only on MIPS.
r306397 (by br):
Use right piece of code for FreeBSD.
r309474:
Don't build :strvis_locale if VIS_NOLOCALE is undefined
The copy of contrib/libc-vis on ^/stable/10 doesn't contain all of the features in the ^/stable/11 // ^/head version, including VIS_NOLOCALE. The risk is lower in conditionally running the test instead of backporting the newer version of libc-vis |
313488 |
09-Feb-2017 |
ngie |
MFC r289172,r290254:
r289172:
Refactor the test/ Makefiles after recent changes to bsd.test.mk (r289158) and netbsd-tests.test.mk (r289151)
- Eliminate explicit OBJTOP/SRCTOP setting - Convert all ad hoc NetBSD test integration over to netbsd-tests.test.mk - Remove unnecessary TESTSDIR setting - Use SRCTOP where possible for clarity
r290254:
Remove unused variable (SRCDIR) |
309740 |
09-Dec-2016 |
emaste |
MFC r309298: libm: remove duplicate version script entries
These symbols already appear in the common lib/msun/Symbol.map. Duplicate entries produce an error with LLVM's LLD linker.
Sponsored by: The FreeBSD Foundation |
309510 |
03-Dec-2016 |
trasz |
MFC r301761:
Fix frexpl() declaration to not include the field name. |
298491 |
22-Apr-2016 |
ngie |
MFC r297836:
Fix appending -O0 to CFLAGS
The previous method would completely nerf CFLAGS once bsd.progs.mk had recursed into the per-PROG logic and make the CFLAGS for tap testcases to -O0, instead of appending to CFLAGS for all of the tap testcases. |
294243 |
18-Jan-2016 |
ngie |
MFC r292497:
Integrate the remaining tools/regression/lib/msun testcases into the FreeBSD test suite under lib/msun/tests |
293267 |
06-Jan-2016 |
ngie |
MFC r292328:
Integrate a number of testcases from tools/regression/lib/msun into the FreeBSD test suite
There's no functional change with these testcases; they're purposely being left in TAP format for the time being
Other testcases which crash on amd64/i386 as-is have not been integrated yet (they need to be retested on a later version of CURRENT, as I haven't used i386 in some time) |
291815 |
04-Dec-2015 |
bdrewery |
MFC r284289,r288267:
r284289: Misc fixes from projects/bmake r288267: Remove redundant .NOPATH. |
286688 |
12-Aug-2015 |
dim |
MFC r286515:
In libm's exp2(3), avoid left-shifting a negative integer, which is undefined. Replace it with the intended value, in a defined way.
Reviewed by: bde |
284810 |
25-Jun-2015 |
tijl |
MFC r271651, r271719, r272138, r272457, r272845, r275476, r275518, r275614, r275819, r276176, r278154, r278160, r278339, r279127, r279240, r279491, r279493, r279856, r283032, r284423, r284426, r284427, r284428
Merge changes to libm from the past 9 months. This includes improvements to the Bessel functions and adds the C99 function lgammal. |
282854 |
13-May-2015 |
ngie |
MFC r282057,r282092,r282106:
r282057:
Build/install libc, librt, libthr, and msun NetBSD test suites on all architectures
r282092 (by andrew):
Correct the spelling of MACHINE_CPUARCH, MACHINE_CPU is not set on arm64.
r282106:
Remove per-architecture checks for enabling HAVE_FENV_H
The conditional came from NetBSD, where only select architectures have this header/support
All architectures on FreeBSD have the necessary support though, so the conditional's completely unnecessary
make tinderbox done on all architectures (including arm64, where the issue occurred before) this time |
280396 |
23-Mar-2015 |
pfg |
MFC r277948: msun: use previously ignored "in" value.
This fixes evaluation of exceptional values in scalblnl(). While here, simplify the code as suggested by Bruce Evans.
Reported by: clang static analyzer |
276672 |
05-Jan-2015 |
ngie |
MFC r274618,r276521:
MFC Note: text corrected for r274618
r274618:
Add reachover Makefiles for contrib/netbsd-tests/lib/msun
A variant of this code has been tested on amd64/i386 for some time by EMC/Isilon on 10-STABLE/11-CURRENT. It builds on other architectures, but the code will remain off until it's proven it works on virtual hardware or real hardware on other architectures
Sponsored by: EMC / Isilon Storage Division
r276521:
Reset errno in :scalbnf_val and :scalbnl_val before running the tests so the tested errno isn't stale
This was needed in order for the test to pass on amd64 with stable/10 |
273827 |
29-Oct-2014 |
andrew |
MFC r269956:
From https://sourceware.org/ml/newlib/2014/msg00113.html By Richard Earnshaw at ARM > >GCC has for a number of years provides a set of pre-defined macros for >use with determining the ISA and features of the target during >pre-processing. However, the design was always somewhat cumbersome in >that each new architecture revision created a new define and then >removed the previous one. This meant that it was necessary to keep >updating the support code simply to recognise a new architecture being >added. > >The ACLE specification (ARM C Language Extentions) >(http://infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.set.swdev/index.html) >provides a much more suitable interface and GCC has supported this >since gcc-4.8. > >This patch makes use of the ACLE pre-defines to map to the internal >feature definitions. To support older versions of GCC a compatibility >header is provided that maps the traditional pre-defines onto the new >ACLE ones.
Stop using __FreeBSD_ARCH_armv6__ and switch to __ARM_ARCH >= 6 in the couple of places in tree. clang already implements ACLE. Add a define that says we implement version 1.1, even though the implementation isn't quite complete. |
271779 |
18-Sep-2014 |
tijl |
MFC r257770 r257818 r257823 r260066 r260067 r260089 r260145 r268587 r268588 r268589 r268590 r268593 r268597 r269758 r270845 r270847 r270893 r270932 r270947 r271147
Merge libm work by kargl, bde and das from the past few months. Besides optimisations and small bug fixes this includes new implementations for C99 functions expl, coshl, sinhl, tanhl, erfl and erfcl.
Approved by: re (kib) |
266314 |
17-May-2014 |
ian |
MFC 262989, 263210, 263230, 263231, 263239, 263242, 263243,
Export _libc_arm_fpu_present as a private symbol to be used by other system libraries, for example libm.
On armv6 access both the softfloat and, when available, the vfp to get and set the floating-point environment.
Build fenv-vfp.c with the softfp float abi. Without this gcc generates an incorrect assembly file that doesn't allow for vfp instructions.
Only build the vfp/softfp switching code on armv6 as we don't support vfp on anything earlier than this. This should fix the armeb and arm builds when using gcc.
Add an optimised version of the float and double helper functions. |
266133 |
15-May-2014 |
ian |
MFC r257207, r261161, r261163:
Update the hard-float version of the fenv functions to use the VFP unit. Any other floating-point unit is unsupported on ARM.
Use __fenv_static for all static inline functions.
Correctly shift the mask when masking/unmasking exceptions. |
263763 |
26-Mar-2014 |
dim |
MFC r262613:
Merge the projects/clang-sparc64 branch back to head. This brings in several updates from the llvm and clang trunks to make the sparc64 backend fully functional.
Apart from one patch to sys/sparc64/include/pcpu.h which is still under discussion, this makes it possible to let clang fully build world and kernel for sparc64.
Any assistance with testing this on actual sparc64 hardware is greatly appreciated, as there will unavoidably be bugs left.
Many thanks go to Roman Divacky for his upstream work on getting the sparc64 backend into shape.
MFC r262985:
Repair a few minor mismerges from r262261 in the clang-sparc64 project branch. This is also to minimize differences with upstream. |
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
|
255499 |
12-Sep-2013 |
dim |
After r255294, building lib/msun's symbol map (using clang as the preprocessor) gives the following error:
--- Version.map --- <stdin>:287:4: error: invalid preprocessing directive # Implemented as weak aliases for imprecise versions ^ 1 error generated.
Change the comment to a C-style one, to prevent this error.
Approved by: re (hrs)
|
255361 |
07-Sep-2013 |
andrew |
On ARM EABI double precision floating point values are stored in the endian the CPU is in, i.e. little-endian on most ARM cores.
This allows ARMv4 and ARMv5 boards to boot with the ARM EABI.
|
255294 |
06-Sep-2013 |
theraven |
Add stub implementations of the missing C++11 math functions. These are weak and so can be replaced by other versions in applications that choose to do so, and will give a linker warning when used so that applications that rely on the extra precision can avoid them.
Note that since the C/C++ specs only guarantee that long double has precision equal to double, code that actually relies on these functions having greater precision is unportable at best and broken at worst.
|
254994 |
28-Aug-2013 |
kargl |
* Whitespace.
|
254969 |
27-Aug-2013 |
kargl |
* s_erf.c: . Use integer literal constants instead of double literal constants.
* s_erff.c: . Use integer literal constants instead of casting double literal constants to float. . Update the threshold values from those carried over from erf() to values appropriate for float. . New sets of polynomial coefficients for the rational approximations. These coefficients have little, but positive, effect on the maximum error in ULP in the four intervals, but do improve the overall speed of execution. . Remove redundant GET_FLOAT_WORD(ix,x) as hx already contained the contents that is packed into ix. . Update the mask that is used to zero-out lower-order bits in x in the intervals [1.25, 2.857143] and [2.857143, 12]. In tests on amd64, this change improves the maximum error in ULP from 6.27739 and 63.8095 to 3.16774 and 2.92095 on these intervals for erffc().
Reviewed by: bde
|
253766 |
29-Jul-2013 |
theraven |
Restore the longer form of the _Generic. The short form does not work in C++.
|
253764 |
29-Jul-2013 |
theraven |
Reenable the isnan(double) / isinf(double) declarations when targeting C89 + SUSv2 mode.
|
253321 |
13-Jul-2013 |
theraven |
Cleaner support for type qualifiers.
Submitted by: Pasi Parviainen
|
253319 |
13-Jul-2013 |
theraven |
Ensure that the _Generic() macro in math.h works with qualified types.
tgmath.h contains the same bugs and so should be fixed in the same way.
|
253260 |
12-Jul-2013 |
theraven |
Fix the build with C++ where __builtin_types_compatible_p is not allowed.
|
253219 |
11-Jul-2013 |
theraven |
Fix some typoes in math.h cleanup.
|
253215 |
11-Jul-2013 |
theraven |
Cleanups to math.h that prevent namespace conflicts with C++.
Reviewed by: bde MFC after: 3 days
|
252170 |
24-Jun-2013 |
eadler |
Make the order of operations for lib/msun more clear.
Tested with md5 sum of object code
Reported by: swildner@DragonFlyBSD.org Submitted by: bde
|
251599 |
10-Jun-2013 |
das |
Add implementations of acoshl(), asinhl(), and atanhl(). This is a merge of the work done by bde and myself.
|
251404 |
05-Jun-2013 |
das |
Style fixes.
Submitted by: bde
|
251366 |
04-Jun-2013 |
das |
Add man links for expl(3) and expm1l(3).
|
251345 |
03-Jun-2013 |
kargl |
Change a comma to a semicolon.
Remove a blank line that crept into the declarations.
Fix a comment to show a sign on a NaN.
|
251343 |
03-Jun-2013 |
kargl |
ld80 and ld128 implementations of expm1l(). This code started life as a fairly faithful implementation of the algorithm found in
PTP Tang, "Table-driven implementation of the Expm1 function in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 18, 211-222 (1992).
Over the last 18-24 months, the code has under gone significant optimization and testing.
Reviewed by: bde Obtained from: bde (most of the optimizations)
|
251341 |
03-Jun-2013 |
kargl |
Fix two comments that got lost in the disentanglement of the larger diff.
|
251339 |
03-Jun-2013 |
kargl |
ld80/s_expl.c:
* Use integral numerical constants, and let the compiler do the conversion to long double.
ld128/s_expl.c:
* Use integral numerical constants, and let the compiler do the conversion to long double. * Use the ENTERI/RETURNI macros, which are no-ops on ld128. This however makes the ld80 and ld128 identical.
Reviewed by: bde (as part of larger diff)
|
251338 |
03-Jun-2013 |
kargl |
Micro-optimization: move the unary mius operator to operate on a literal constant.
Obtained from: bde
|
251337 |
03-Jun-2013 |
kargl |
Add a comment to note that bde supplied most, if not all, of the optimizations.
|
251335 |
03-Jun-2013 |
kargl |
ld80/s_expl.c:
* In the special case x = -Inf or -NaN, use a micro-optimization to eliminate the need to access u.xbits.man.
* Fix an off-by-one for small arguments |x| < 0x1p-65.
ld128/s_expl.c:
* In the special case x = -Inf or -NaN, use a micro-optimization to eliminate the need to access u.xbits.manh and u.xbits.manl.
* Fix an off-by-one for small arguments |x| < 0x1p-114.
Obtained from: bde
|
251334 |
03-Jun-2013 |
kargl |
ld80/s_expl.c:
* Update the evaluation of the polynomial. This allows the removal of the now unused variables t23 and t45.
ld128/s_expl.c:
* Update the evaluation of the polynomial and the intermediate result t. This update allows several numerical constants to be written as double rather than long double constants. Update the constants as appropriate.
Obtained from: bde
|
251333 |
03-Jun-2013 |
kargl |
Rename a few P2, P3, ... coefficients to A2, A3, ... missed in my previous commit.
|
251330 |
03-Jun-2013 |
kargl |
Update a comment to reflect that we are using an endpoint of an interval instead of a midpoint.
|
251328 |
03-Jun-2013 |
kargl |
Add a u suffix to the IEEEl2bits unions o_threshold and u_threshold, and use macros to access the e component of the unions. This allows the portions of the code in ld80 to be identical to the ld128 code.
Obtained from: bde
|
251327 |
03-Jun-2013 |
kargl |
Introduce the macro LOG2_INTERVAL, which is log2(number of intervals). Use the macroi as a micro-optimization to convert a subtraction and division to a shift.
Obtained from: bde
|
251325 |
03-Jun-2013 |
kargl |
Whitespace.
|
251321 |
03-Jun-2013 |
kargl |
* Rename the polynomial coefficients from P2, P3, ... to A2, A3, .... The names now coincide with the name used in PTP Tang's paper.
* Rename the variable from s to tbl to better reflect that this is a table, and to be consistent with the naming scheme in s_exp2l.c
Reviewed by: bde (as part of larger diff)
|
251316 |
03-Jun-2013 |
kargl |
* Style(9). Start non-Copyright fancy formatted comments with /**.
Reviewed by: bde (as part of larger diff)
|
251315 |
03-Jun-2013 |
kargl |
ld80/s_expl.c:
* Update Copyright years to include 2013.
ld128/s_expl.c:
* Correct and update Copyright years. This code originated from the ld80 version, so it should reflect the same time period.
Reviewed by: bde (as part of larger diff)
|
251292 |
03-Jun-2013 |
das |
Add logl, log2l, log10l, and log1pl.
Submitted by: bde
|
251121 |
30-May-2013 |
das |
I'm happy to finally commit stephen@'s implementations of cacos, cacosh, casin, casinh, catan, and catanh. Thanks to stephen@ and bde@ for working on these.
Submitted by: stephen@ Reviewed by: bde
|
251046 |
27-May-2013 |
kargl |
Style(9)
Approved by: das (implicit) Reported by: jh
|
251041 |
27-May-2013 |
kargl |
* Update polynomial coefficients. * Use ENTERI/RETURNI to allow the use of FP_PE on i386 target.
Reviewed by: das (and bde a long time ago) Approved by: das (mentor) Obtained from: bde (polynomial coefficients)
|
251024 |
27-May-2013 |
das |
Fix some regressions caused by the switch from gcc to clang. The fixes are workarounds for various symptoms of the problem described in clang bugs 3929, 8100, 8241, 10409, and 12958.
The regression tests did their job: they failed, someone brought it up on the mailing lists, and then the issue got ignored for 6 months. Oops. There may still be some regressions for functions we don't have test coverage for yet.
|
249724 |
21-Apr-2013 |
kib |
Merge the 386 and amd64 versions of the fenv.h, to make cc -m32 compilations which use fenv.h work.
Reviewed by: tjil Sponsored by: The FreeBSD Foundation
|
242879 |
10-Nov-2012 |
dim |
Only define isnan, isnanf, __isnan and __isnanf in libc.so, not in libc.a and libc_p.a. In addition, define isnan in libm.a and libm_p.a, but not in libm.so.
This makes it possible to statically link executables using both isnan and isnanf with libc and libm.
Tested by: kargl MFC after: 1 week
|
242692 |
07-Nov-2012 |
kevlo |
Fix typo; s/ouput/output
|
241887 |
22-Oct-2012 |
imp |
Revert r241756
|
241886 |
22-Oct-2012 |
imp |
Revert r241755
|
241756 |
19-Oct-2012 |
imp |
Document the method used to compute expf. Taken from exp, with changes to reflect differences in computation between the two.
|
241755 |
19-Oct-2012 |
imp |
Document the methods used to compute logf. Taken and edited from the double version, with adaptations for the differences between it and the float version.
|
241545 |
14-Oct-2012 |
joel |
mdoc: avoid nested displays.
|
241516 |
13-Oct-2012 |
kargl |
* Update the comment that explains the choice of values in the table and the requirement on trailing zero bits.
* Remove the __aligned() compiler directives as these were found to have a negative effect on the produced code.
Submitted by: bde Approved by: das (mentor)
|
241051 |
29-Sep-2012 |
kargl |
* src/math_private.h: . Change the API for the LD80C by removing the explicit passing of the sign bit. The sign can be determined from the last parameter of the macro. . On i386, load long double by bit manipulations to work around at least a gcc compiler issue. On non-i386 ld80 architectures, use a simple assignment.
* ld80/s_expl.c: . Update the only consumer of LD80C.
Submitted by: bde Approved by: das (mentor)
|
240866 |
23-Sep-2012 |
kargl |
* ld80/s_expl.c: . Fix the threshold for expl(x) where |x| is small. . Also update the previously incorrect comment to match the new threshold.
* ld128/s_expl.c: . Re-order logic in exceptional cases to match the logic used in other long double functions. . Fix the threshold for expl(x) where is |x| is small. . Also update the previously incorrect comment to match the new threshold.
Submitted by: bde Approved by: das (mentor)
|
240865 |
23-Sep-2012 |
kargl |
Fix whitespace issue.
Approved by: das (mentor, implicit)
|
240864 |
23-Sep-2012 |
kargl |
* ld80/s_expl.c: . Guard a comment from reformatting by indent(1). . Re-order variables in declarations to alphabetical order. . Remove a banal comment.
* ld128/s_expl.c: . Add a comment to point to ld80/s_expl.c for implementation details. . Move the #define of INTERVAL to reduce the diff with ld80/s_expl.c. . twom10000 does not need to be volatile, so move its declaration. . Re-order variables in declarations to alphabetical order. . Add a comment that describes the argument reduction. . Remove the same banal comment found in ld80/s_expl.c.
Reviewed by: bde Approved by: das (mentor)
|
240861 |
23-Sep-2012 |
kargl |
* Update the lookup table to use 53-bit high and low values. Also, update the comment to describe the choice of using a high and low decomposition of 2^(i/INTERNVAL) for 0 <= i <= INTERVAL in preparation for an implementation of expm1l.
* Move the #define of INTERVAL above the comment, because the comment refers to INTERVAL.
Reviewed by: bde Approved by: das (mentor)
|
240828 |
22-Sep-2012 |
kargl |
* Use ENTERI() and RETURNI() to toggle the rounding precision if necessary, so that cosl(), sinl() and tanl() work on i386 even for naive callers.
Suggested by: bde Reviewed by: bde Approved by: das (mentor)
|
240827 |
22-Sep-2012 |
kargl |
* Make STRICT_ALIGN() work for doubles as well as for floats. This only affects i386. The double case was intentionally left broken as an optimization, but we are getting closer to supporting applications and/or kernels that change the (FreeBSD i386) default rounding precision from FP_PD to FP_PE and never change it back, and this requires the STRICT_ALIGN()s that were added to support FP_PE to actually work in all precisions.
* Remove an extraneous semicolon at the end of a macro that was supposed to be function-like.
Submitted by: bde Approved by: das (mentor)
|
240566 |
16-Sep-2012 |
tijl |
Optimise i387 trigonometric functions. Replace "andw 0x400,%ax \ jnz" with "sahf \ jp", "fprem1" with "fprem" and "fstsw %ax" with "fnstsw %ax".
|
240368 |
11-Sep-2012 |
joel |
Minor mdoc fix.
|
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
|
238924 |
30-Jul-2012 |
kargl |
ieeefp.h is only needed on i386 class hardware.
Submitted by: bde Approved by: das (pre-approved)
|
238923 |
30-Jul-2012 |
kargl |
Whitespace.
Submitted by: bde Approved by: das (pre-approved)
|
238784 |
26-Jul-2012 |
kargl |
Replace the macro name NUM with INTERVALS. This change provides compatibility with the INTERVALS macro used in the soon-to-be-commmitted expm1l() and someday-to-be-committed log*l() functions.
Add a comment into ld128/s_expl.c noting at gcc issue that was deleted when rewriting ld80/e_expl.c as ld128/s_expl.c.
Requested by: bde Approved by: das (mentor)
|
238783 |
26-Jul-2012 |
kargl |
* ld80/expl.c: . Remove a few #ifdefs that should have been removed in the initial commit. . Sort fpmath.h to its rightful place.
* ld128/s_expl.c: . Replace EXPMASK with its actual value. . Sort fpmath.h to its rightful place.
Requested by: bde Approved by: das (mentor)
|
238782 |
26-Jul-2012 |
kargl |
Replace code that toggles between 53 and 64 bits on i386 class hardware with the ENTERI and RETURNI macros, which are now available in math_private.h.
Suggested by: bde Approved by: das (mentor)
|
238724 |
23-Jul-2012 |
kargl |
Hook ld80/s_expl.c or ld128/s_expl.c into the building of libm.
PR: standards/152415 Approved by: das (mentor)
|
238722 |
23-Jul-2012 |
kargl |
Compute the exponential of x for Intel 80-bit format and IEEE 128-bit format. These implementations are based on
PTP Tang, "Table-driven implementation of the exponential function in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 144-157 (1989).
PR: standards/152415 Submitted by: kargl Reviewed by: bde, das Approved by: das (mentor)
|
236148 |
27-May-2012 |
theraven |
Allow inclusion of libc++ <cmath> to work after including math.h
Submitted by: Yamaya Takashi Reviewed by: das MFC after: 1 week
|
235286 |
11-May-2012 |
gjb |
General mdoc(7) and typo fixes.
PR: 167734 Submitted by: Nobuyuki Koganemaru (kogane!jp.freebsd.org) MFC after: 3 days
|
234685 |
25-Apr-2012 |
des |
I stopped using my middle name years ago.
|
233992 |
07-Apr-2012 |
joel |
mdoc: fix column names, indentation, column separation within each row, and quotation. Also make sure we have the same amount of columns in each row as the number of columns we specify in the head arguments.
Reviewed by: brueffer
|
233973 |
07-Apr-2012 |
das |
Fix a bug in remquo{,f,l}, in which the quotient didn't always have the correct sign when the remainder was 0.
Fix a separate bug in remquo alone, in which the remainder and quotient were both off by a bit in certain cases involving subnormal remainders.
The bugs affected all platforms except amd64 and i386, on which the routines are implemented in assembly.
PR: 166463 Submitted by: Ilya Burylov MFC after: 2 weeks
|
233462 |
25-Mar-2012 |
joel |
Remove superfluous paragraph macro.
|
230371 |
20-Jan-2012 |
das |
Fix a small nit noted by bde: exp_x should be of type float, not double.
|
230370 |
20-Jan-2012 |
das |
Add a change I missed in r230367 (don't inline arm's fenv.h functions).
|
230367 |
20-Jan-2012 |
das |
Don't inline fenv.h functions on arm for now. Inlining makes sense: the function bodies require only 2 to 10 instructions. However, it leads to application binaries that refer to a private ABI, namely, the softfloat innards in libc. This could complicate future changes in the implementation of the floating-point emulation layer, so it seems best to have programs refer to the official fe* entry points in libm.
|
230192 |
16-Jan-2012 |
das |
Add an implementation of fenv.h intended for platforms that lack an FPU and use softfloat.
Thanks to Ian Lepore for testing and debugging this patch. The fenv regression tests pass (at least for Ian's arm chip) with this change.
|
229784 |
07-Jan-2012 |
uqs |
Convert files to UTF-8 and add some copyright markers where missing.
|
227472 |
12-Nov-2011 |
theraven |
Expose the unimplemented libm functions in the math.h header. This allows C++'s <cmath> to work without the compiler complaining that the C++ versions are calling implicitly-declared functions. You will still get a linker error when they are called. OpenBSD 5.0 claims to fully implement the C99 <math.h> stuff, so might be worth investigating...
Reviewed by: das Approved by: dim (mentor)
|
226618 |
21-Oct-2011 |
das |
Minor corrections and clarifications regarding exceptions.
|
226601 |
21-Oct-2011 |
das |
Fix a regression introduced in r226371: When the high part of x*y exactly cancels with z, return the low part of x*y instead of discarding it.
|
226600 |
21-Oct-2011 |
das |
Fix a corner case: tan(large + Inf i) == NaN + NaN i.
|
226599 |
21-Oct-2011 |
das |
Improved handling of large x in ccosh{,f}(): - Handle cases where exp(x) would overflow, but ccosh(x) ~= exp(x) / 2 shouldn't. - Use the ccosh(x) ~= exp(x) / 2 approximation to simplify the calculation when x is large.
Similarly for csinh(). Also fixed the return value of csinh(-Inf +- 0i).
|
226598 |
21-Oct-2011 |
das |
Use __ldexp_exp() to simplify things and improve accuracy for x near the overflow threshold.
|
226597 |
21-Oct-2011 |
das |
The cexp() and {,c}{cos,sin}h functions all need to be able to compute exp(x) scaled down by some factor, and the challenge is doing this accurately when exp(x) would overflow. This change replaces all of the tricks we've been using with common __ldexp_exp() and __ldexp_cexp() routines that handle all the scaling.
bde plans to improve on this further by moving the guts of exp() into k_exp.c and handling the scaling in a more direct manner. But the current approach is simple and adequate for now.
|
226596 |
21-Oct-2011 |
das |
Use STRICT_ASSIGN() to ensure that the compiler doesn't screw things up by storing x in a wider type than it's supposed to.
Submitted by: bde
|
226595 |
21-Oct-2011 |
das |
Per IEEE754r, pow(1, y) is 1 even if y is NaN, and pow(-1, +-Inf) is 1.
MFC after: 2 weeks
|
226594 |
21-Oct-2011 |
das |
Bugfix: feenableexcept() and fedisableexcept() should just return the old exception mask, not mask | ~FE_ALL_EXCEPT.
MFC after: 2 weeks
|
226460 |
17-Oct-2011 |
das |
It's no longer accurate to say that math.h "constitute[s] the C math library," since complex.h, tgmath.h, and fenv.h are also part of the math library. Replace the outdated sentence with some references to the other parts.
|
226458 |
17-Oct-2011 |
das |
Add c{cos,sin,tan}{,h}{,f} functions. This is joint work with bde and kargl.
|
226436 |
16-Oct-2011 |
eadler |
- change "is is" to "is" or "it is" - change "the the" to "the"
Approved by: lstewart Approved by: sahil (mentor) MFC after: 3 days
|
226415 |
16-Oct-2011 |
das |
Use #include "fenv.h" instead of #include <fenv.h>. This makes it more convenient to compile the math library by itself.
Requested by: bde
|
226414 |
16-Oct-2011 |
das |
Fix some non-standard variable declarations.
|
226413 |
16-Oct-2011 |
das |
Optimize the case of pure imaginary arguments. Calls like this are common, e.g., in DFT implementations.
Discussed with: bde, kargl
|
226412 |
16-Oct-2011 |
das |
Move the macros GET_LDBL_EXPSIGN() and SET_LDBL_EXPSIGN() into math_private.h, so they can be used elsewhere in the math library.
|
226411 |
16-Oct-2011 |
das |
Remove an unused variable.
|
226380 |
15-Oct-2011 |
das |
Remove some unnecessary initializations.
Obtained from: DragonFlyBSD
|
226376 |
15-Oct-2011 |
das |
Various changes to improve the accuracy and speed of log{2,10}{,f}. - Rename __kernel_log() to k_log1p(). - Move some of the work that was previously done in the kernel log into the callers. This enables further refactoring to improve accuracy or speed, although I don't recall the details. - Use extra precision when adding the final scaling term, which improves accuracy. - Describe and work around compiler problems that break some of the multiprecision calculations.
A fix for a small bug is also included: - Add a special case for log*(1). This is needed to ensure that log*(1) == +0 instead of -0, even when the rounding mode is FE_DOWNWARD.
Submitted by: bde
|
226375 |
15-Oct-2011 |
das |
Style fixes and updates to comments.
Submitted by: bde
|
226374 |
15-Oct-2011 |
das |
Don't define FP_FAST_FMA on sparc64; with the recent fixes, fma() is no longer "fast" on sparc64. (It really wasn't to begin with, since the old implementation was using long doubles, and long doubles are emulated in software on sparc64.)
|
226373 |
15-Oct-2011 |
das |
Add INSERT_WORD64 and EXTRACT_WORD64 macros for use in s_fma.c.
|
226372 |
15-Oct-2011 |
das |
Replace two lines accidentally removed in r226218. Thanks to bde for noticing this.
|
226371 |
15-Oct-2011 |
das |
Fix a double-rounding bug in fma{,f,l}. The bug would occur in round-to-nearest mode when the result, rounded to twice machine precision, was exactly halfway between two machine-precision values. The essence of the fix is to simulate a "sticky bit" in the pathological cases, which is how hardware implementations break the ties.
MFC after: 1 month
|
226245 |
11-Oct-2011 |
das |
Refactor this code by introducing separate functions to handle the extra-precision add and multiply operations. This simplifies future work but shouldn't result in any functional change.
|
226218 |
10-Oct-2011 |
das |
Provide external definitions of all of the standardized functions in fenv.h that are currently inlined.
The definitions are provided in fenv.c via 'extern inline' declaractions. This assumes the compiler handles 'extern inline' as specified in C99, which has been true under FreeBSD since 8.0.
The goal is to eventually remove the 'static' keyword from the inline definitions in fenv.h, so that non-inlined references all wind up pointing to the same external definition like they're supposed to. I am deferring the second step to provide a window where newly-compiled apps will still link against old math libraries. (This isn't supported, but there's no need to cause undue breakage.)
Reviewed by: stefanf, bde
|
223302 |
19-Jun-2011 |
kargl |
In the libm access macros for the double type, z can sometimes be used uninitialized. This can lead to spurious exceptions and bit clobbering.
Submitted by: bde Approved by: das (mentor)
|
223262 |
18-Jun-2011 |
benl |
Fix clang warnings.
Approved by: philip (mentor)
|
222508 |
30-May-2011 |
kargl |
Clean up the unneeded cpp macro INLINE_REM_PIO2L.
Reviewed by: das Approved by: das (mentor)
|
221234 |
29-Apr-2011 |
kargl |
Improve the accuracy from a max ULP of ~2000 to max ULP < 0.79 on i386-class hardware for sinl and cosl. The hand-rolled argument reduction have been replaced by e_rem_pio2l() implementations. To preserve history the following commands have been executed:
svn cp src/e_rem_pio2.c ld80/e_rem_pio2l.h mv ${HOME}/bde/ld80/e_rem_pio2l.c ld80/e_rem_pio2l.h
svn cp src/e_rem_pio2.c ld128/e_rem_pio2l.h mv ${HOME}/bde/ld128/e_rem_pio2l.c ld128/e_rem_pio2l.h
The ld80 version has been tested by bde, das, and kargl over the last few years (bde, das) and few months (kargl). An older ld128 version was tested by das. The committed version has only been compiled tested via 'make universe'.
Approved by: das (mentor) Obtained from: bde
|
219576 |
12-Mar-2011 |
kargl |
Take two. Add the missing file that should have been committed with r219571 and re-enable building of cbrtl.
Implement the long double version for the cube root function, cbrtl. The algorithm uses Newton's iterations with a crude estimate of the cube root to converge to a result.
Reviewed by: bde Approved by: das
|
219572 |
12-Mar-2011 |
kargl |
Temporary disable the building of cbrtl until I can determine why svn will not allow one to commit a new file.
Approved by: das (implicit)
|
219571 |
12-Mar-2011 |
kargl |
Implement the long double version for the cube root function, cbrtl. The algorithm uses Newton's iterations with a crude estimate of the cube root to converge to a result.
Reviewed by: bde Approved by: das
|
219366 |
07-Mar-2011 |
das |
Add cexp() to the complex(3) manpage. Thanks to bde for pointing out that I missed this.
|
219365 |
07-Mar-2011 |
das |
Remove part of an uncommitted change that snuck into the last commit.
|
219361 |
07-Mar-2011 |
das |
Convert log10f() to use __kernel_log(), which is more accurate and simpler.
|
219360 |
07-Mar-2011 |
das |
Convert log10() to use __kernel_log(), which is more accurate and simpler.
|
219359 |
07-Mar-2011 |
das |
Add cexp() and cexpf().
Reviewed by: bde (earlier version)
|
218909 |
21-Feb-2011 |
brucec |
Fix typos - remove duplicate "the".
PR: bin/154928 Submitted by: Eitan Adler <lists at eitanadler.com> MFC after: 3 days
|
218877 |
20-Feb-2011 |
murray |
Add complex(3) manual page documenting our partial support for C99 complex arithmetic in libm.
Reviewed by: David Schultz <das@FreeBSD.org> MFC after: 2 weeks
|
218511 |
10-Feb-2011 |
das |
Fix a bug where the wrong argument was passed to SET_FLOAT_WORD(). This bug results in a type mismatch that happens to be harmless because of the way SET_FLOAT_WORD() works.
Submitted by: bde
|
218510 |
10-Feb-2011 |
das |
Fix a bug where the wrong argument was passed to INSERT_WORDS(). This bug results in a type mismatch that happens to be harmless because of the way INSERT_WORDS() works.
Submitted by: bde
|
218509 |
10-Feb-2011 |
das |
For small arguments, these functions use simple approximations, e.g. cos(small) = 1, sin(small) = small. This commit tightens the thresholds at which the simple approximations are used.
Reviewed by: bde
|
218508 |
10-Feb-2011 |
das |
Fix a bogus threshold that was copied from the double precision version. This commit should have no effect on correctness; it merely changes the threshold at which a simpler approximation can be used.
Reviewed by: bde
|
218305 |
04-Feb-2011 |
kib |
Remove duplicate .note.GNU-stack section declaration.
Reported by: arundel
|
217108 |
07-Jan-2011 |
kib |
Add section .note.GNU-stack for assembly files used by 386 and amd64.
|
216248 |
07-Dec-2010 |
das |
Another minor nit: Make sure the constant here is a float so the compiler doesn't promote the entire expression to double.
|
216247 |
07-Dec-2010 |
das |
Fix various nits in style and comments that were pointed out by bde. Code changes verified with md5.
|
216211 |
05-Dec-2010 |
das |
Add log2() and log2f().
|
216210 |
05-Dec-2010 |
das |
Add a "kernel" log function, based on e_log.c, which is useful for 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).
|
216137 |
03-Dec-2010 |
das |
Disable gcc's built-in rint() function when compiling s_nearbyint.c. It results in incorrect optimizations that break nearbyint().
PR: 143358 Reviewed by: bde
|
215237 |
13-Nov-2010 |
uqs |
Fix bug in jn(3) and jnf(3) that led to -inf results
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
|
213337 |
01-Oct-2010 |
dim |
Use __FBSDID() instead of RCSID() in most .S files under lib/msun/i386, and one under lib/msun/amd64. This avoids adding the identifiers to the .text section, and moves them to the .comment section instead.
Suggested by: bde Approved by: rpaulo (mentor)
|
212531 |
13-Sep-2010 |
imp |
This is exactly the same as the .else, so remove it.
|
212518 |
13-Sep-2010 |
imp |
MFtbemd: Move to using MACHINE_CPUARCH, now that it is safe.
|
211934 |
28-Aug-2010 |
nwhitehorn |
Repair some build breakage introduced in r211725 and garbage collect some code made obsolete in the same commit.
|
211725 |
23-Aug-2010 |
imp |
MFtbemd:
Prefer MACHNE_CPUARCH to MACHINE_ARCH in most contexts where you want to test of all the CPUs of a given family conform.
|
211397 |
16-Aug-2010 |
joel |
Fix typos, spelling, formatting and mdoc mistakes found by Nobuyuki while translating these manual pages. Minor corrections by me.
Submitted by: Nobuyuki Koganemaru <n-kogane@syd.odn.ne.jp>
|
209877 |
10-Jul-2010 |
nwhitehorn |
powerpc64 floating-point is identical to powerpc, so use the same code on both architectures.
|
209110 |
12-Jun-2010 |
das |
Introduce __isnanf() as an alias for isnanf(), and make the isnan() 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>
|
208935 |
09-Jun-2010 |
uqs |
mdoc: spell out theta, the Unicode glyph is hard to read for terminal fonts
It is referred to as "theta" later in the document anyway, so stop being fancy.
|
208734 |
02-Jun-2010 |
uqs |
mdoc: spell macros correctly, there's no need for the backslash escape
|
208594 |
27-May-2010 |
uqs |
mdoc: Garbage collect unused/unneeded macros
|
208028 |
13-May-2010 |
uqs |
mdoc: move remaining sections into consistent order
This pertains mostly to FILES, HISTORY, EXIT STATUS and AUTHORS sections.
Found by: mdocml lint run Reviewed by: ru
|
208027 |
13-May-2010 |
uqs |
mdoc: move CAVEATS, BUGS and SECURITY CONSIDERATIONS sections to the bottom of the manpages and order them consistently.
GNU groff doesn't care about the ordering, and doesn't even mention CAVEATS and SECURITY CONSIDERATIONS as common sections and where to put them.
Found by: mdocml lint run Reviewed by: ru
|
205076 |
12-Mar-2010 |
uqs |
Fix several typos in macros or macro misusage.
Found by: make manlint Reviewed by: ru Approved by: philip (mentor)
|
203441 |
03-Feb-2010 |
kib |
Placate new binutils, by using 16-bit %ax instead of 32-bit %eax as an argument for fnstsw. Explicitely specify sizes for the XMM control and status word and X87 control and status words.
Reviewed by: das Tested by: avg MFC after: 2 weeks
|
194000 |
11-Jun-2009 |
ed |
Use the documented machine constraint for SSE registers.
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/
|
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.
|
192760 |
25-May-2009 |
attilio |
Use, in uncovered part, the END() macro in order to improve debugging. In this specific case, Valgrind won't get confused when analyzing such functions.
Sponsored by: Sandvine Incorporated Tested by: emaste MFC: 3 days
|
189805 |
14-Mar-2009 |
das |
Namespace: scalb() is withdrawn from POSIX.
|
189803 |
14-Mar-2009 |
das |
Eliminate __real__ and __imag__ gccisms.
|
188272 |
07-Feb-2009 |
das |
C99 TC2 now wants FP_FAST_FMA* to be defined to 1, if the macros are defined at all. See also: defect report #223.
|
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
|
186886 |
08-Jan-2009 |
das |
Fix the types of INFINITY and NAN, which were broken in r131851. They should both be floats, not doubles.
PR: 127795 Submitted by: Christoph Mallon MFC after: 2 weeks
|
186461 |
23-Dec-2008 |
marcel |
Add support for the FPA floating-point format on ARM. The 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
|
183714 |
09-Oct-2008 |
peter |
Clean out some empty mergeinfo records, presumably by people doing local cp/mv operations. The full repo-relative URL should be specified for the source in these cases.
|
181405 |
08-Aug-2008 |
das |
Remove some unused variables.
Reported by: Intel C Compiler
|
181402 |
08-Aug-2008 |
das |
In the line #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.
|
181377 |
07-Aug-2008 |
das |
Implement cproj{,f,l}().
|
181374 |
07-Aug-2008 |
das |
Use cpack() and the gcc extension __imag__ to implement cimag() and 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.
|
181258 |
03-Aug-2008 |
das |
Fix some style bogosity from fdlibm.
|
181257 |
03-Aug-2008 |
das |
Minor improvements: - Improve the order of some tests. - Fix style.
Submitted by: bde
|
181204 |
02-Aug-2008 |
das |
A few minor corrections, including some from bde: - 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.
|
181152 |
02-Aug-2008 |
das |
On i386, gcc truncates long double constants to double precision 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
|
181100 |
01-Aug-2008 |
das |
Fix some problems with asinf(), acosf(), atanf(), and atan2f():
- 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.
|
181074 |
31-Jul-2008 |
das |
Add implementations of acosl(), asinl(), atanl(), atan2l(), and cargl().
Reviewed by: bde sparc64 testing resources from: remko
|
181064 |
31-Jul-2008 |
das |
Set WARNS=1.
I believe I've committed all the bits necessary to make this compile on all supported architectures. :crosses fingers:
|
181063 |
31-Jul-2008 |
das |
The high part of the mantissa is 64 bits on sparc64.
|
181062 |
31-Jul-2008 |
das |
As in other parts of libm, mark a few constants as volatile to prevent 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.
|
180581 |
18-Jul-2008 |
das |
Sort the .PATH entries to give a more reasonable order of precedence: 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).
|
180074 |
28-Jun-2008 |
das |
Fix a typo in the cosl() prototype.
|
179882 |
19-Jun-2008 |
das |
Implement fmodl. Document fmodl and fix some errors in the fmod manpage.
|
178749 |
03-May-2008 |
gonzo |
Symbol.map is handled by cpp, so use C-style comments
Approved by: cognet (mentor)
|
178582 |
26-Apr-2008 |
imp |
Add mips support to libm, from mips2-jnpr perforce branch.
|
177875 |
03-Apr-2008 |
das |
Fix some corner cases: - 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.
|
177793 |
31-Mar-2008 |
das |
Remove a (bogus) remnant of debugging this on sparc64.
|
177768 |
30-Mar-2008 |
das |
Add assembly versions of remquol() and remainderl().
|
177766 |
30-Mar-2008 |
das |
Hook remquol() and remainderl() up to the build.
|
177765 |
30-Mar-2008 |
das |
Implement remainderl() as a wrapper around remquol(). The extra work remquol() performs to compute the quotient is negligible.
|
177764 |
30-Mar-2008 |
das |
Implement remquol() based on remquo().
|
177761 |
30-Mar-2008 |
das |
Implement csqrtl().
|
177760 |
30-Mar-2008 |
das |
Hook hypotl() and cabsl() up to the build.
|
177759 |
30-Mar-2008 |
das |
Document hypotl().
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
|
177758 |
30-Mar-2008 |
das |
Alias hypotl() and cabsl() for platforms where long double is the same as double.
|
177757 |
30-Mar-2008 |
das |
Implement cabsl() in terms of hypotl().
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
|
177756 |
30-Mar-2008 |
das |
Implement hypotl(). This is bde's conversion of fdlibm hypot(), with minor fixes for ld128 by me.
|
177751 |
30-Mar-2008 |
bde |
Use fabs[f]() instead of bit fiddling for setting absolute values. 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.
|
177749 |
30-Mar-2008 |
bde |
Use the expression fabs(x+0.0)-fabs(y+0.0) instead of 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.)
|
177746 |
30-Mar-2008 |
bde |
Fix a missing mask in a hi+lo decomposition. Thus bug made the extra 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).
|
177712 |
29-Mar-2008 |
das |
Include math.h for the fmaf() prototype.
|
177711 |
29-Mar-2008 |
das |
Fix some rather obscene code that has ambiguous if...if...else... constructs in it.
|
176747 |
02-Mar-2008 |
das |
1 << 47 needs to be written 1ULL << 47.
|
176722 |
02-Mar-2008 |
das |
Hook up sqrtl() to the build.
|
176721 |
02-Mar-2008 |
das |
MD implementations of sqrtl().
|
176720 |
02-Mar-2008 |
das |
MI implementation of sqrtl(). This is very slow and should be overridden when hardware sqrt is available.
|
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.
|
176530 |
24-Feb-2008 |
raj |
Let PowerPC world optionally build with -msoft-float. For FPU-less PowerPC 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
|
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).
|
176466 |
22-Feb-2008 |
bde |
Remove the "quick check no cancellation" optimization for 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.
|
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.
|
176462 |
22-Feb-2008 |
bde |
Add an irint() function in inline asm for amd64 and i386. irint() is 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.
|
176461 |
22-Feb-2008 |
bde |
Optimize the conversion to bits a little (by about 11 cycles or 16% 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.
|
176458 |
22-Feb-2008 |
bde |
Optimize the fixup for +-0 by using better classification for this case 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.
|
176456 |
22-Feb-2008 |
bde |
Fix rintl() on signaling NaNs and unsupported formats.
|
176451 |
22-Feb-2008 |
das |
s/rcsid/__FBSDID/
|
176450 |
22-Feb-2008 |
das |
Remove an unused variable.
|
176449 |
22-Feb-2008 |
das |
Eliminate some warnings.
|
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.
|
176409 |
19-Feb-2008 |
bde |
Optimize for 3pi/4 <= |x| <= 9pi/4 in much the same way as for pi/4 <= |x| <= 3pi/4. Use the same branch ladder as for float precision. Remove the optimization for |x| near pi/2 and don't do it near the multiples of pi/2 in the newly optimized range, since it requires fairly large code to handle only relativley few cases. Ifdef out optimization for |x| <= pi/4 since this case can't occur because it is done in callers.
On amd64 (A64), for cos() and sin() with uniformly distributed args, no cache misses, some parallelism in the caller, and good but not great CC and CFLAGS, etc., this saves about 40 cycles or 38% in the newly optimized range, or about 27% on average across the range |x| <= 2pi (~65 cycles for most args, while the A64 hardware fcos and fsin take ~75 cycles for half the args and 125 cycles for the other half). The speedup for tan() is much smaller, especially relatively. The speedup on i386 (A64) is slightly smaller, especially relatively. i386 is still much slower than amd64 here (unlike in the float case where it is slightly faster).
|
176408 |
19-Feb-2008 |
bde |
Rearrange the polynomial evaluation for better parallelism. This saves an average of about 8 cycles or 5% on A64 (amd64 and i386 -- more in cycles but about the same percentage on i386, and more with old versions of gcc) with good CFLAGS and some parallelism in the caller. As usual, it takes a couple more multiplications so it will be slower on old machines.
Convert to __FBSDID().
|
176390 |
18-Feb-2008 |
das |
Document return values better.
|
176388 |
18-Feb-2008 |
das |
Add tgammaf() as a simple wrapper around tgamma().
|
176387 |
18-Feb-2008 |
bde |
2 long double constants were missing L suffixes. This helped break tanl() on !(amd64 || i386). It gave slightly worse than double precision in some cases. tanl() now passes tests of 2^24 values on ia64.
|
176386 |
18-Feb-2008 |
bde |
Fix a typo which broke k_tanl.c on !(amd64 || i386).
|
176385 |
18-Feb-2008 |
bde |
Inline __ieee754__rem_pio2(). With gcc4-2, this gives an average optimization of about 10% for cos(x), sin(x) and tan(x) on |x| < 2**19*pi/2. We didn't do this before because __ieee754__rem_pio2() is too large and complicated for gcc-3.3 to inline very well. We don't do this for float precision because it interferes with optimization of the usual (?) case (|x| < 9pi/4) which is manually inlined for float precision only.
This has some rough edges: - some static data is duplicated unnecessarily. There isn't much after the recent move of large tables to k_rem_pio2.c, and some static data is duplicated to good affect (all the data static const, so that the compiler can evaluate expressions like 2*pio2 at compile time and generate even more static data for the constant for this). - extern inline is used (for the same reason as in previous inlining of k_cosf.c etc.), but C99 apparently doesn't allow extern inline functions with static data, and gcc will eventually warn about this.
Convert to __FBSDID().
Indent __ieee754_rem_pio2()'s declaration consistently (its style was made inconsistent with fdlibm a while ago, so complete this).
Fix __ieee754_rem_pio2()'s return type to match its prototype. Someone changed too many ints to int32_t's when fixing the assumption that all ints are int32_t's.
|
176373 |
17-Feb-2008 |
das |
Use volatile hacks to make sure exp() generates an underflow exception when it's supposed to. Previously, gcc -O2 was optimizing away the statement that generated it.
|
176361 |
17-Feb-2008 |
das |
Hook up sinl(), cosl(), and tanl() to the build.
|
176360 |
17-Feb-2008 |
das |
Add implementations of sinl(), cosl(), and tanl().
Submitted by: Steve Kargl <sgk@apl.washington.edu>
|
176359 |
17-Feb-2008 |
das |
Documentation for sinl(), cosl(), and tanl().
|
176358 |
17-Feb-2008 |
das |
Add kernel functions for 128-bit long doubles. These could be improved a bit, but access to a freebsd/sparc64 machine is needed.
Submitted by: bde and Steve Kargl <sgk@apl.washington.edu> (earlier version)
|
176357 |
17-Feb-2008 |
das |
Add kernel functions for 80-bit long doubles. Many thanks to Steve and Bruce for putting lots of effort into these; getting them right isn't easy, and they went through many iterations.
Submitted by: Steve Kargl <sgk@apl.washington.edu> with revisions from bde
|
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.
|
176305 |
15-Feb-2008 |
bde |
Sigh, the weak reference for ceill(), floorl() and truncl() was in unreachable code due to a missing include. This kept arm and powerpc broken.
Reported by: sam, grehan
|
176280 |
14-Feb-2008 |
bde |
Oops, the weak reference for ceill(), floorl() and truncl() was in the wrong file. This broke arm and powerpc.
Reported by: grehan
|
176277 |
14-Feb-2008 |
bde |
Use the expression fabs(x+0.0)+fabs(y+0.0) instad of a+b (where a is |x| or |y| and b is |y| or |x|) when mixing NaN arg(s).
hypot*() had its own foot shooting for mixing NaNs -- it swaps the args so that |x| in bits is largest, but does this before quieting signaling NaNs, so on amd64 (where the result of adding NaNs depends on the order) it gets inconsistent results if setting the quiet bit makes a difference, just like a similar ia64 and i387 hardware comparison. The usual fix (see e_powf.c 1.13 for more details) of mixing using (a+0.0)+-(b+0.0) doesn't work on amd64 if the args are swapped (since the rder makes a difference with SSE). Fortunately, the original args are unchanged and don't need to be swapped when we let the hardware decide the mixing after quieting them, but we need to take their absolute value.
hypotf() doesn't seem to have any real bugs masked by this non-bug. On amd64, its maximum error in 2^32 trials on amd64 is now 0.8422 ulps, and on i386 the maximum error is unchanged and about the same, except with certain CFLAGS it magically drops to 0.5 (perfect rounding).
Convert to __FBSDID().
|
176268 |
14-Feb-2008 |
bde |
Fix the hi+lo decomposition for 2/(3ln2). The decomposition needs to be into 12+24 bits of precision for extra-precision multiplication, but was into 13+24 bits. On i386 with -O1 the bug was hidden by accidental extra precision, but on amd64, in 2^32 trials the bug caused about 200000 errors of more than 1 ulp, with a maximum error of about 80 ulps. Now the maximum error in 2^32 trials on amd64 is 0.8573 ulps. It is still 0.8316 ulps on i386 with -O1.
The nearby decomposition of 1/ln2 and the decomposition of 2/(3ln2) in the double precision version seem to be sub-optimal but not broken.
|
176266 |
14-Feb-2008 |
bde |
Use the expression (x+0.0)-(y+0.0) instead of x+y when mixing NaN arg(s). This uses 2 tricks to improve consistency so that more serious problems aren't hidden in simple regression tests by noise for the NaNs:
- for a signaling NaN, adding 0.0 generates the invalid exception and converts to a quiet NaN, and doesn't have too many effects for other types of args (it converts -0 to +0 in some rounding modes, but that hopefully doesn't change the result after adding the NaN arg). This avoids some inconsistencies on i386 and ia64. On these arches, the result of an operation on 2 NaNs is apparently the largest or the smallest of the NaNs as bits (consistently largest or smallest for each arch, but the opposite). I forget which way the comparison goes and if the sign bit affects it. The quiet bit is is handled poorly by not always setting it before the comparision or ignoring it. Thus if one of the args was originally a signaling NaN and the other was originally a quiet NaN, then the result depends too much on whether the signaling NaN has been quieted at this point, which in turn depends on optimizations and promotions. E.g., passing float signaling NaNs to double functions must quiet them on conversion; on i387, loading a signaling NaN of type float or double (but not long double) into a register involves a conversion, so it quiets signaling NaNs, so if the addition has 2 register operands than it only sees quiet NaNs, but if the addition has a memory operand then it sees a signaling NaN iff it is in the memory operand.
- subtraction instead of addition is used to avoid a dubious optimization in old versions of gcc. For SSE operations, mixing of NaNs apparently always gives the target operand. This is not as good as the i387 and ia64 behaviour. It doesn't mix NaNs at all, and makes addition not quite commutative. Old versions of gcc sometimes rewrite x+y to y+x and thus give different results (in bits) for NaNs. gcc-3.3.3 rewrites x+y to y+x for one of pow() and powf() but not the other, so starting from float NaN args x and y, powf(x, y) was almost always different from pow(x, y).
These tricks won't give consistency of 2-arg float and double functions with long double ones on amd64, since long double ones use the i387 which has different semantics from SSE.
Convert to __FBSDID().
|
176245 |
13-Feb-2008 |
bde |
s_ceill.c s_floorl.c s_truncl.c
|
176243 |
13-Feb-2008 |
bde |
On arches where long double is the same as double, alias ceil(), floor() and trunc() to the corresponding long double functions. This is not just an optimization for these arches. The full long double functions have a wrong value for `huge', and the arches without full long doubles depended on it being wrong.
|
176236 |
13-Feb-2008 |
bde |
Fix the C version of ceill(x) for -1 < x <= -0 in all rounding modes. The result should be -0, but was +0.
|
176231 |
13-Feb-2008 |
bde |
Fix exp2*(x) on signaling NaNs by returning x+x as usual.
This has the side effect of confusing gcc-4.2.1's optimizer into more often doing the right thing. When it does the wrong thing here, it seems to be mainly making too many copies of x with dependency chains. This effect is tiny on amd64, but in some cases on i386 it is enormous. E.g., on i386 (A64) with -O1, the current version of exp2() should take about 50 cycles, but took 83 cycles before this change and 66 cycles after this change. exp2f() with -O1 only speeded up from 51 to 47 cycles. (exp2f() should take about 40 cycles, on an Athlon in either i386 or amd64 mode, and now takes 42 on amd64). exp2l() with -O1 slowed down from 155 cycles to 123 for some args; this is unimportant since the i386 exp2l() is a fake; the wrong thing for it seems to involve branch misprediction.
|
176229 |
13-Feb-2008 |
bde |
Rearrange the polynomial evaluation for better parallelism. This is faster on all machines tested (old Celeron (P2), A64 (amd64 and i386) and ia64) except on ia64 when compiled with -O1. It takes 2 more multiplications, so it will be slower on old machines. The speedup is about 8 cycles = 17% on A64 (amd64 and i386) with best CFLAGS and some parallelism in the caller.
Move the evaluation of 2**k up a bit so that it doesn't compete too much with the new polynomial evaluation. Unlike the previous optimization, this rearrangement cannot change the result, so compilers and CPU schedulers can do it, but they don't do it quite right yet. This saves a whole 1 or 2 cycles on A64.
|
176227 |
13-Feb-2008 |
bde |
Use hardware remainder on amd64 since it is 5 to 10 times faster than software remainder and is already used for remquo().
|
176207 |
12-Feb-2008 |
bde |
Fix remainder() and remainderf() in round-towards-minus-infinity mode when the result is +-0. IEEE754 requires (in all rounding modes) that if the result is +-0 then its sign is the same as that of the first arg, but in round-towards-minus-infinity mode an uncorrected implementation detail always reversed the sign. (The detail is that x-x with x's sign positive gives -0 in this mode only, but the algorithm assumed that x-x always has positive sign for such x.)
remquo() and remquof() seem to need the same fix, but I cannot test them yet.
Use long doubles when mixing NaN args. This trick improves consistency of results on at least amd64, so that more serious problems like the above aren't hidden in simple regression tests by noise for the NaNs. On amd64, hardware remainder should be used since it is about 10 times faster than software remainder and is already used for remquo(), but it involves using the i387 even for floats and doubles, and the i387 does NaN mixing which is better than but inconsistent with SSE NaN mixing. Software remainder() would probably have been inconsistent with software remainderl() for the same reason if the latter existed.
Signaling NaNs cause further inconsistencies on at least ia64 and i386.
Use __FBSDID().
|
176168 |
11-Feb-2008 |
bde |
Use double precision for z and thus for the entire calculation of exp2(i/TBLSIZE) * p(z) instead of only for the final multiplication and addition. This fixes the code to match the comment that the maximum error is 0.5010 ulps (except on machines that evaluate float expressions in extra precision, e.g., i386's, where the evaluation was already in extra precision).
Fix and expand the comment about use of double precision.
The relative roundoff error from evaluating p(z) in non-extra precision was about 16 times larger than in exp2() because the interval length is 16 times smaller. Its maximum was at least P1 * (1.0 ulps) * max(|z|) ~= log(2) * 1.0 * 1/32 ~= 0.0217 ulps (1.0 ulps from the addition in (1 + P1*z) with a cancelation error when z ~= -1/32). The actual final maximum was 0.5313 ulps, of which 0.0303 ulps must have come from the additional roundoff error in p(z). I can't explain why the additional roundoff error was almost 3/2 times larger than the rough estimate.
|
176132 |
09-Feb-2008 |
bde |
As usual, use a minimax polynomial that is specialized for float precision. The new polynomial has degree 4 instead of 10, and a maximum error of 2**-30.04 ulps instead of 2**-33.15. This doesn't affect the final error significantly; the maximum error was and is about 0.5015 ulps on i386 -O1, and the number of cases with an error of > 0.5 ulps is increased from 13851 to 14407.
Note that the error is only this close to 0.5 ulps due to excessive extra precision caused by compiler bugs on i386. The extra precision could be obtained intentionally, and is useful for keeping the error of the hyperbolic float functions below 1 ulp, since these functions are implemented using expm1f. My recent change for scaling by 2**k had the unintentional side effect of retaining extra precision for longer, so callers of expm1f see errors of more like 0.0015 ulps than 0.5015 ulps, and for the hyperbolic functions this reduces the maximum error from nearly about 2 ulps to about 0.75 ulps.
This is about 10% faster on i386 (A64). expm1* is still very slow, but now the float version is actually significantly faster. The algorithm is very sophisticated but not very good except on machines with fast division.
|
176128 |
09-Feb-2008 |
bde |
Fix a comment about coefficients and expand a related one.
|
176102 |
08-Feb-2008 |
bde |
Fix truncl() when the result should be -0.0L. When the result is +-0.0L, it must have the same sign as the arg in all rounding modes, but it was always +0.0L.
|
176101 |
08-Feb-2008 |
bde |
Oops, fix the fix in rev.1.10. logb() and logbf() were broken on denormals, and logb() remained broken after 1.10 because the fix for logbf() was incompletely translated.
Convert to __FBSDID().
|
176082 |
07-Feb-2008 |
bde |
Use a better method of scaling by 2**k. Instead of adding to the exponent bits of the reduced result, construct 2**k (hopefully in parallel with the construction of the reduced result) and multiply by it. This tends to be much faster if the construction of 2**k is actually in parallel, and might be faster even with no parallelism since adjustment of the exponent requires a read-modify-wrtite at an unfortunate time for pipelines.
In some cases involving exp2* on amd64 (A64), this change saves about 40 cycles or 30%. I think it is inherently only about 12 cycles faster in these cases and the rest of the speedup is from partly-accidentally avoiding compiler pessimizations (the construction of 2**k is now manually scheduled for good results, and -O2 doesn't always mess this up). In most cases on amd64 (A64) and i386 (A64) the speedup is about 20 cycles. The worst case that I found is expf on ia64 where this change is a pessimization of about 10 cycles or 5%. The manual scheduling for plain exp[f] is harder and not as tuned.
Details specific to expm1*: - the saving is closer to 12 cycles than to 40 for expm1* on i386 (A64). For some reason it is much larger for negative args. - also convert to __FBSDID().
|
176074 |
07-Feb-2008 |
bde |
Use a better method of scaling by 2**k. Instead of adding to the exponent bits of the reduced result, construct 2**k (hopefully in parallel with the construction of the reduced result) and multiply by it. This tends to be much faster if the construction of 2**k is actually in parallel, and might be faster even with no parallelism since adjustment of the exponent requires a read-modify-wrtite at an unfortunate time for pipelines.
In some cases involving exp2* on amd64 (A64), this change saves about 40 cycles or 30%. I think it is inherently only about 12 cycles faster in these cases and the rest of the speedup is from partly-accidentally avoiding compiler pessimizations (the construction of 2**k is now manually scheduled for good results, and -O2 doesn't always mess this up). In most cases on amd64 (A64) and i386 (A64) the speedup is about 20 cycles. The worst case that I found is expf on ia64 where this change is a pessimization of about 10 cycles or 5%. The manual scheduling for plain exp[f] is harder and not as tuned.
This change ld128/s_exp2l.c has not been tested.
|
176032 |
06-Feb-2008 |
bde |
As for the float trig functions and logf, use a minimax polynomial that is specialized for float precision. The new polynomial has degree 5 instead of 11, and a maximum error of 2**-27.74 ulps instead of 2**-30.64. This doesn't affect the final error significantly; the maximum error was and is about 0.9101 ulps on amd64 -01 and the number of cases with an error of > 0.5 ulps is actually reduced by epsilon despite the larger error in the polynomial.
This is about 15% faster on amd64 (A64), i386 (A64) and ia64. The asm version is still used instead of this on i386 since it is faster and more accurate.
|
175731 |
28-Jan-2008 |
das |
Adjust the exponent before converting the result from double to float precision. This fixes some double rounding problems for subnormals and simplifies things a bit.
|
175665 |
25-Jan-2008 |
bde |
Fix a harmless type error in 1.9.
|
175534 |
21-Jan-2008 |
bde |
Fix cutoffs. This is just a cleanup and an optimization for unusual cases which are used mainly by regression tests.
As usual, the cutoff for tiny args was not correctly translated to float precision. It was 2**-54 but 2**-24 works. It must be about 2**-precision, since the error from approximating log(1+x) by x is about the same as |x|. Exhaustive testing shows that 2**-24 gives perfect rounding in round-to-nearest mode.
Similarly for the cutoff for being small, except this is not used by so many other functions. It was 2**-29 but 2**-15 works. It must be a bit smaller than sqrt(2**-precision), since the error from approximating log(1+x) by x-x*x/2 is about the same as x*x. Exhaustive testing shows that 2**-15 gives a maximum error of 0.5052 ulps in round-to-nearest-mode. The algorithm for the general case is only good for 0.8388 ulps, so this is sufficient (but it loses slightly on i386 -- then extra precision gives 0.5032 ulps for the general case).
While investigating this, I noticed that optimizing the usual case by falling into a middle case involving a simple polynomial evaluation (return x-x*x/2 instead of x here) is not such a good idea since it gives an enormous pessimization of tinier args on machines for which denormals are slow. Float x*x/2 is denormal when |x| ~< 2**-64 and x*x/2 is evaluated in float precision, so it can easily be denormal for normal x. This is even more interesting for general polynomial evaluations. Multiplying out large powers of x is normally a good optimization since it reduces dependencies, but it creates denormals starting with quite large x.
|
175507 |
20-Jan-2008 |
bde |
Oops, when merging from the float version to the double versions, don't forget to translate "float" to "double".
ucbtest didn't detect the bug, but exhaustive testing of the float case relative to the double case eventually did. The bug only affects args x with |x| ~> 2**19*(pi/2) on non-i386 (i386 is broken in a different way for large args).
|
175505 |
19-Jan-2008 |
bde |
Remove the float version of the kernel of arg reduction for pi/2, since it should never have existed and it has not been used for many years (floats are reduced faster using doubles). All relevant changes (just the workaround for broken assignment) have been merged to the double version.
|
175503 |
19-Jan-2008 |
bde |
Do an ordinary assignment in STRICT_ASSIGN() except for floats until there is a problem with non-floats (when i386 defaults to extra precision). This essentially restores yesterday's behaviour for doubles on i386 (since generic rint() isn't used and everywhere else assumed working assignment), but for arches that use the generic rint() it finishes restoring some of 1995's behaviour (don't waste time doing unnecessary store/load).
|
175501 |
19-Jan-2008 |
bde |
Use STRICT_ASSIGN() for exp2f() and exp2() instead of a volatile variable hack for exp2f() only.
The volatile variable had a surprisingly large cost for exp2f() -- 19 cycles or 15% on i386 in the worst case observed. This is only partly explained by there being several references to the variable, only one of which benefited from it being volatile. Arches that have working assignment are likely to benefit even more from not having any volatile variable.
exp2() now has a chance of working with extra precision on i386.
exp2() has even more references to the variable, so it would have been pessimized more by simply declaring the variable as volatile. Even the temporary volatile variable for STRICT_ASSIGN costs 5-10% on i386, (A64) so I will change STRICT_ASSIGN() to do an ordinary assignment until i386 defaults to extra precision.
|
175499 |
19-Jan-2008 |
bde |
Use STRICT_ASSIGN() for _kernel_rem_pio2f() and _kernel_rem_pio2f() instead of a volatile cast hack for the float version only. The cast hack broke with gcc-4, but this was harmless since the float version hasn't been used for a few years. Merge from the float version so that the double version has a chance of working on i386 with extra precision.
See k_rem_pio2f.c rev.1.8 for the original hack.
Convert to _FBSDID().
|
175494 |
19-Jan-2008 |
bde |
Use STRICT_ASSIGN() for log1pf() and log1p() instead of a volatile cast hack for log1pf() only. The cast hack broke with gcc-4, resulting in ~1 million errors of more than 1 ulp, with a maximum error of ~1.5 ulps. Now the maximum error for log1pf() on i386 is 0.5034 ulps again (this depends on extra precision), and log1p() has a chance of working with extra precision.
See s_log1pf.c 1.8 for the original hack. (It claims only 62343 large errors).
Convert to _FBSDID(). Another thing broken with gcc-4 is the static const hack used for rcsids.
|
175480 |
19-Jan-2008 |
bde |
Use STRICT_ASSIGN() instead of assorted direct volatile hacks to work around assignments not working for gcc on i386. Now volatile hacks for rint() and rintf() don't needlessly pessimize so many arches and the remaining pessimizations (for arm and powerpc) can be avoided centrally.
This cleans up after s_rint.c 1.3 and 1.13 and s_rintf.c 1.3 and 1.9: - s_rint.c 1.13 broke 1.3 by only using a volatile cast hack in 1 place when it was needed in 2 places, and the volatile cast hack stopped working with gcc-4. These bugs only affected correctness tests on i386 since i386 normally uses asm rint() and doesn't support the extra precision mode that would break assignments of doubles. - s_rintf.c 1.9 improved(?) on 1.3 by using a volatile variable hack instead of an extra-precision variable hack, but it declared 2 variables as volatile when only 1 variable needed to be volatile. This only affected speed tests on i386 since i386 uses asm rintf().
|
175468 |
18-Jan-2008 |
das |
Use volatile hacks to make sure these functions generate an underflow exception when they're supposed to. Previously, gcc -O2 was optimizing away the statement that generated it.
|
175462 |
18-Jan-2008 |
das |
Hook up exp2l() and related docs to the build.
|
175461 |
18-Jan-2008 |
das |
Introduce a new log(3) manpage and move the relevant functions there. Document exp2l() in exp(3), and remove the quaint discussion of topics such as what these functions were called on the HP-71B's variant of BASIC.
|
175460 |
18-Jan-2008 |
das |
Implement exp2l(). There is one version for machines with 80-bit long doubles (i386, amd64, ia64) and one for machines with 128-bit long doubles (sparc64). Other platforms use the double version. I've only done runtime testing on i386.
Thanks to bde@ for helpful discussions and bugfixes.
|
175403 |
17-Jan-2008 |
bde |
Add a macro STRICT_ASSIGN() to help avoid the compiler bug that assignments and casts don't clip extra precision, if any. The implementation is to assign to a temporary volatile variable and read the result back to assign to the original lvalue.
lib/msun currently 2 different hard-coded hacks to avoid the problem in just a few places and needs it in a few more places. One variant uses volatile for the original lvalue. This works but is slower than necessary. Another temporarily casts the lvalue to volatile. This broke with gcc-4.2.1 or earlier (gcc now stores to the lvalue but doesn't load from it).
|
175371 |
15-Jan-2008 |
das |
Optimize this a bit better.
Submitted by: bde (although these aren't all of his changes)
|
175309 |
14-Jan-2008 |
das |
Implement rintl(), nearbyintl(), lrintl(), and llrintl(). Thanks to bde@ for feedback and testing of rintl().
|
175225 |
11-Jan-2008 |
das |
- Correct the range check in the double version to catch negative values that would overflow. - Style fixes and improved handling of NaNs suggested by bde.
|
174804 |
20-Dec-2007 |
das |
Grumble. DO declare logbl(), DON'T declare logl() just yet. bde is going to commit logl() Real Soon Now. I'm just trying to slow him down with merge conflicts.
Noticed by: bde
|
174801 |
20-Dec-2007 |
das |
Remove the declaration of logl(). The relevant bits haven't been committed yet, but the declaration leaked in when I added nan() and friends.
Reported by: pav
|
174759 |
18-Dec-2007 |
das |
Since nan() is supposed to work the same as strtod("nan(...)", NULL), my original implementation made both use the same code. Unfortunately, this meant libm depended on a vendor header at compile time and previously- unexposed vendor bits in libc at runtime.
Hence, I just wrote my own version of the relevant vendor routine. As it turns out, mine has a factor of 8 fewer of lines of code, and is a bit more readable anyway. The strtod() and *scanf() routines still use vendor code.
Reviewed by: bde
|
174732 |
18-Dec-2007 |
das |
Remove z_abs(). The z_*() functions were in libf77, and for some reason someone thought it would be a good idea to copy z_abs() to libm in 1994. However, it's never been declared or documented anywhere, and I'm reasonably confident that nobody uses it.
Discussed with: bde, deischen, kan
|
174720 |
17-Dec-2007 |
bde |
Oops, the previous commit was not needed -- the file was committed but not checked out due to my checkout error.
|
174719 |
17-Dec-2007 |
bde |
Translate from the i386 so that this compiles and runs.
I hope that this and the i386 version of it will not be needed, but this is currently about 16 cycles or 36% faster than the C version, and the i386 version is about 8 cycles or 19% faster than the C version, due to poor optimization of the C version.
|
174715 |
17-Dec-2007 |
bde |
Don't try to build s_nanl.c before it is committed.
|
174698 |
17-Dec-2007 |
das |
Add logbl(3) to libm.
|
174694 |
17-Dec-2007 |
das |
Document the fact that we have nan(3) now, and make some minor clarifications in other places.
|
174684 |
16-Dec-2007 |
das |
Implement and document nan(), nanf(), and nanl(). This commit adds two new directories in msun: ld80 and ld128. These are for long double functions specific to the 80-bit long double format used on x86-derived architectures, and the 128-bit format used on sparc64, respectively.
|
174618 |
15-Dec-2007 |
das |
1. Add csqrt{,f}(3). 2. Put carg{,f}(3) under the FBSD_1.1 namespace where it belongs (requested by kan@)
|
174617 |
15-Dec-2007 |
das |
Implement and document csqrt(3) and csqrtf(3).
|
174586 |
14-Dec-2007 |
das |
Update the standards section, and make a minor clarification about the return value of sqrt.
|
174584 |
14-Dec-2007 |
das |
Typo in previous commit
|
174583 |
14-Dec-2007 |
das |
Symbol.map additions for carg and cargf. (They're in C99, so I didn't add a new version for them.)
|
174563 |
12-Dec-2007 |
das |
s/C90/C99/
|
174562 |
12-Dec-2007 |
das |
Add a "STANDARDS" section.
|
174561 |
12-Dec-2007 |
das |
Implement carg(3) and cargf(3).
Rotting in an old src tree since: March 2005
|
170707 |
14-Jun-2007 |
bde |
Oops, back out previous commit since it was backwards to a wrong branch.
|
170706 |
14-Jun-2007 |
bde |
MFC: 1.11: fix the threshold for (not) using the simple Taylor approximation.
|
170550 |
11-Jun-2007 |
bde |
Fix an aliasing bug which was finally detected by gcc-4.2. fdlibm has hundreds of similar aliasing bugs, but all except this one seem to have been fixed by Cygnus and/or NetBSD before the modified version of fdlibm was imported into FreeBSD in 1994.
PR: standards/113147 Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
|
170089 |
29-May-2007 |
bde |
Merge the relevant part of rev.1.14 of s_cbrt.c (a micro-optimization involving moving the check for x == 0). The savings in cycles are smaller for cbrtf() than for cbrt(), and positive in all measured cases with gcc-3.4.4, but still very machine/compiler-dependent.
|
169807 |
21-May-2007 |
deischen |
Bump library versions in preparation for 7.0.
Ok'd by: kan
|
169524 |
13-May-2007 |
deischen |
Enable symbol versioning by default. Use WITHOUT_SYMVER to disable it. Warning, after symbol versioning is enabled, going back is not easy (use WITHOUT_SYMVER at your own risk).
Change the default thread library to libthr.
There most likely still needs to be a version bump for at least the thread libraries. If necessary, this will happen later.
|
169220 |
02-May-2007 |
bde |
Don't assume that int is signed 32-bits in one place. Keep assuming that ints have >= 31 value bits elsewhere. s/int/int32_t/ seems to have been done too globally for all other files in msun/src before msun/ was imported into FreeBSD.
Minor fixes in comments.
e_lgamma_r.c: Describe special cases in more detail: - exception for lgamma(0) and lgamma(neg.integer) - lgamma(-Inf) = Inf. This is wrong but is required by C99 Annex F. I hope to change this.
|
169212 |
02-May-2007 |
bde |
Fix tgamma() on some special args: (1) tgamma(-Inf) returned +Inf and failed to raise any exception, but should always have raised an exception, and should behave like tgamma(negative integer). (2) tgamma(negative integer) returned +Inf and raised divide-by-zero, but should return NaN and raise "invalid" on any IEEEish system. (3) About half of the 2**52 negative intgers between -2**53 and -2**52 were misclassified as non-integers by using floor(x + 0.5) to round to nearest, so tgamma(x) was wrong (+-0 instead of +Inf and now NaN) on these args. The floor() expression is hard to use since rounding of (x + 0.5) may give x or x + 1, depending on |x| and the current rounding mode. The fixed version uses ceil(x) to classify x before operating on x and ends up being more efficient since ceil(x) is needed anyway. (4) On at least the problematic args in (3), tgamma() raised a spurious inexact. (5) tgamma(large positive) raised divide-by-zero but should raise overflow. (6) tgamma(+Inf) raised divide-by-zero but should not raise any exception. (7) Raise inexact for tiny |x| in a way that has some chance of not being optimized away.
The fix for (5) and (6), and probably for (2), also prevents -O optimizing away the exception.
PR: 112180 (2) Standards: Annex F in C99 (IEC 60559 binding) requires (1), (2) and (6).
|
169209 |
02-May-2007 |
bde |
Document (in a comment) the current (slightly broken) handling of special values in more detail, and change the style of this comment to be closer to fdlibm and C99: - tgamma(-Inf) was undocumented and is wrong (+Inf, should be NaN) - tgamma(negative integer) is as intended (+Inf) but not best for IEEE-754 (NaN) - tgamma(-0) was documented as being wrong (+Inf) but was correct (-Inf) - documentation of setting of exceptions (overflow, etc.) was more complete here than in most of libm, but was further from matching the actual setting than in most of libm, due to various bugs here (primarily, always evaluating +Inf one/zero and getting unwanted divide-by-zero exceptions from this). Now the actual behaviour with gcc -O0 is documented. Optimization still breaks setting of exceptions all over libm, so nothing can depend on this working. - tgamma(NaN)'s exception was documented as being wrong (invalid) but was correct (no exception with IEEEish NaNs).
Finish (?) rev.1.5. gamma was not renamed to tgamma in one place.
Finish (?) rev.1.6. errno.h was not completely removed.
|
169092 |
29-Apr-2007 |
deischen |
Use C comments since we now preprocess these files with CPP.
|
165906 |
09-Jan-2007 |
imp |
Remove California Regent's clause 3, per letter
|
165855 |
07-Jan-2007 |
das |
Implement modfl().
|
165841 |
06-Jan-2007 |
das |
Fix a problem relating to fesetenv() clobbering i387 register stack.
Details: As a side-effect of restoring a saved FP environment, fesetenv() overwrites the tag word, which indicates which i387 registers are in use. Normally this isn't a problem because the calling convention requires the register stack to be empty on function entry and exit. However, fesetenv() is inlined, so we need to tell gcc explicitly that the i387 registers get clobbered.
PR: 85101
|
165840 |
06-Jan-2007 |
das |
Fix a cut-and-paste-o.
|
165839 |
06-Jan-2007 |
das |
Correctly handle NaN.
|
165838 |
06-Jan-2007 |
das |
Correctly handle inf/nan. This routine is currently unused because we seem to have assembly versions for all architectures, but it can't hurt to fix it.
|
165837 |
06-Jan-2007 |
das |
Remove modf from libm's symbol map. It's actually in libc for hysterical raisins.
|
165795 |
05-Jan-2007 |
das |
Remove an unneeded fnstcw instruction.
Noticed by: bde
|
165794 |
05-Jan-2007 |
das |
Remove a note pertaining to the Alpha.
|
163358 |
14-Oct-2006 |
bde |
Moved __BEGIN_DECLS up a little so that it covers __test_sse() and C++ isn't broken,
PR: 104425
|
161526 |
22-Aug-2006 |
ru |
Remove alpha left-overs.
|
160151 |
07-Jul-2006 |
bde |
Fixed the threshold for using the simple Taylor approximation.
In e_log.c, there was just a off-by-1 (1 ulp) error in the comment about the threshold. The precision of the threshold is unimportant, but the magic numbers in the code are easier to understand when the threshold is described precisely.
In e_logf.c, mistranslation of the magic numbers gave an off-by-1 (1 * 16 ulps) error in the intended negative bound for the threshold and an off-by-7 (7 * 16 ulps) error in the intended positive bound for the threshold, and the intended bounds were not translated from the double precision bounds so they were unnecessarily small by a factor of about 2048.
The optimization of using the simple Taylor approximation for args near a power of 2 is dubious since it only applies to a relatively small proportion of args, but if it is done then doing it 2048 times as often _may_ be more efficient. (My benchmarks show unexplained dependencies on the data that increase with further optimizations in this area.)
|
160122 |
05-Jul-2006 |
bde |
Fixed tanh(-0.0) on ia64 and optimizeed tanh(x) for 2**-55 <= |x| < 2**-28 as a side effect, by merging with the float precision version of tanh() and the double precision version of sinh().
For tiny x, tanh(x) ~= x, and we used the expression x*(one+x) to return this value (x) and set the inexact flag iff x != 0. This doesn't work on ia64 since gcc -O does the dubious optimization x*(one+x) = x+x*x so as to use fma, so the sign of -0.0 was lost.
Instead, handle tiny x in the same as sinh(), although this is imperfect: - return x directly and set the inexact flag in a less efficient way. - increased the threshold for non-tinyness from 2**-55 to 2**-28 so that many more cases are optimized than are pessimized.
Updated some comments and fixed bugs in others (ranges for half-open intervals mostly had the open end backwards, and there were nearby style bugs).
|
160118 |
05-Jul-2006 |
bde |
Removed the optimized asm versions of scalb() and scalbf(). These functions are only for compatibility with obsolete standards. They shouldn't be used, so they shouldn't be optimized. Use the generic versions instead.
This fixes scalbf() as a side effect. The optimized asm version left garbage on the FP stack. I fixed the corresponding bug in the optimized asm scalb() and scalbn() in 1996. NetBSD fixed it in scalb(), scalbn() and scalbnf() in 1999 but missed fixing it in scalbf(). Then in 2005 the bug was reimplemented in FreeBSD by importing NetBSD's scalbf().
The generic versions have slightly different error handling: - the asm versions blindly round the second parameter to a (floating point) integer and proceed, while the generic versions return NaN if this rounding changes the value. POSIX permits both behaviours (these functions are XSI extensions and the behaviour for a bogus non-integral second parameter is unspecified). Apart from this and the bug in scalbf(), the behaviour of the generic versions seems to be identical. (I only exhusatively tested generic_scalbf(1.0F, anyfloat) == asm_scalb(1.0F, anyfloat). This covers many representative corner cases involving NaNs and Infs but doesn't test exception flags. The brokenness of scalbf() showed up as weird behaviour after testing just 7 integer cases sequentially.)
|
160102 |
05-Jul-2006 |
bde |
Backed out rev.1.10. It tried to implement ldexpf() as a weak reference to scalbf(), but ldexpf() cannot be implemented in that way since the types of the second parameter differ. ldexpf() can be implemented as a weak or strong reference to scalbnf() (*) but that was already done long before rev.1.10 was committed. The old implementation uses a reference, so rev.1.10 had no effect on applications. The C files for the scalb() family are not used for amd64 or i386, so rev.1.10 had even less effect for these arches.
(*) scalbnf() raises the radix to the given exponent, while ldexpf() raises 2 to the given exponent. Thus the functions are equivalent except possibly for their error handling iff the radix is 2. Standards more or less require identical error handling. Under FreeBSD, the functions are equivalent except for more details being missing in scalbnf()'s man page.
|
157196 |
27-Mar-2006 |
deischen |
Add symbol versioning to libm.
|
154051 |
05-Jan-2006 |
bde |
Oops, on amd64 (and probably on all non-i386 systems), the previous commit broke the 2**24 cases where |x| > DBL_MAX/2. There are exponent range problems not just for denormals (underflow) but for large values (overflow). Doubles have more than enough exponent range to avoid the problems, but I forgot to convert enough terms to double, so there was an x+x term which was sometimes evaluated in float precision.
Unfortunately, this is a pessimization with some combinations of systems and compilers (it makes no difference on Athlon XP's, but on Athlon64's it gives a 5% pessimization with gcc-3.4 but not with gcc-3.3).
Exlain the problem better in comments.
|
154049 |
05-Jan-2006 |
bde |
Use double precision internally to optimize cbrtf(), and change the algorithm for the second step significantly to also get a perfectly rounded result in round-to-nearest mode. The resulting optimization is about 25% on Athlon64's and 30% on Athlon XP's (about 25 cycles out of 100 on the former).
Using extra precision, we don't need to do anything special to avoid large rounding errors in the third step (Newton's method), so we can regroup terms to avoid a division, increase clarity, and increase opportunities for parallelism. Rearrangement for parallelism loses the increase in clarity. We end up with the same number of operations but with a division reduced to a multiplication.
Using specifically double precision, there is enough extra precision for the third step to give enough precision for perfect rounding to float precision provided the previous steps are accurate to 16 bits. (They were accurate to 12 bits, which was almost minimal for imperfect rounding in the old version but would be more than enough for imperfect rounding in this version (9 bits would be enough now).) I couldn't find any significant time optimizations from optimizing the previous steps, so I decided to optimize for accuracy instead. The second step needed a division although a previous commit optimized it to use a polynomial approximation for its main detail, and this division dominated the time for the second step. Use the same Newton's method for the second step as for the third step since this is insignificantly slower than the division plus the polynomial (now that Newton's method only needs 1 division), significantly more accurate, and simpler. Single precision would be precise enough for the second step, but doesn't have enough exponent range to handle denormals without the special grouping of terms (as in previous versions) that requires another division, so we use double precision for both the second and third steps.
|
153548 |
20-Dec-2005 |
bde |
Extract the high and low words together. With gcc-3.4 on uniformly distributed non-large args, this saves about 14 of 134 cycles for Athlon64s and about 5 of 199 cycles for AthlonXPs.
Moved the check for x == 0 inside the check for subnormals. With gcc-3.4 on uniformly distributed non-large args, this saves another 5 cycles on Athlon64s and loses 1 cycle on AthlonXPs.
Use INSERT_WORDS() and not SET_HIGH_WORD() when converting the first approximation from bits to double. With gcc-3.4 on uniformly distributed non-large args, this saves another 4 cycles on both Athlon64s and and AthlonXPs.
Accessing doubles as 2 words may be an optimization on old CPUs, but on current CPUs it tends to cause extra operations and pipeline stalls, especially for writes, even when only 1 of the words needs to be accessed.
Removed an unused variable.
|
153520 |
19-Dec-2005 |
bde |
Use a minimax polynomial approximation instead of a Pade rational function approximation for the second step. The polynomial has degree 2 for cbrtf() and 4 for cbrt(). These degrees are minimal for the final accuracy to be essentially the same as before (slightly smaller). Adjust the rounding between steps 2 and 3 to match. Unfortunately, for cbrt(), this breaks the claimed accuracy slightly although incorrect rounding doesn't. Claim less accuracy since its not worth pessimizing the polynomial or relying on exhaustive testing to get insignificantly more accuracy.
This saves about 30 cycles on Athlons (mainly by avoiding 2 divisions) so it gives an overall optimization in the 10-25% range (a larger percentage for float precision, especially in 32-bit mode, since other overheads are more dominant for double precision, surprisingly more in 32-bit mode).
|
153517 |
18-Dec-2005 |
bde |
Fixed code to match comments and the algorithm: - in preparing for the third approximation, actually make t larger in magnitude than cbrt(x). After chopping, t must be incremented by 2 ulps to make it larger, not 1 ulp since chopping can reduce it by almost 1 ulp and it might already be up to half a different-sized-ulp smaller than cbrt(x). I have not found any cases where this is essential, but the think-time error bound depends on it. The relative smallness of the different-sized-ulp limited the bug. If there are cases where this is essential, then the final error bound would be 5/6+epsilon instead of of 4/6+epsilon ulps (still < 1). - in preparing for the third approximation, round more carefully (but still sloppily to avoid branches) so that the claimed error bound of 0.667 ulps is satisfied in all cases tested for cbrt() and remains satisfied in all cases for cbrtf(). There isn't enough spare precision for very sloppy rounding to work: - in cbrt(), even with the inadequate increment, the actual error was 0.6685 in some cases, and correcting the increment increased this a little. The fix uses sloppy rounding to 25 bits instead of very sloppy rounding to 21 bits, and starts using uint64_t instead of 2 words for bit manipulation so that rounding more bits is not much costly. - in cbrtf(), the 0.667 bound was already satisfied even with the inadequate increment, but change the code to almost match cbrt() anyway. There is not enough spare precision in the Newton approximation to double the inadequate increment without exceeding the 0.667 bound, and no spare precision to avoid this problem as in cbrt(). The fix is to round using an increment of 2 smaller-ulps before chopping so that an increment of 1 ulp is enough. In cbrt(), we essentially do the same, but move the chop point so that the increment of 1 is not needed.
Fixed comments to match code: - in cbrt(), the second approximation is good to 25 bits, not quite 26 bits. - in cbrt(), don't claim that the second approximation may be implemented in single precision. Single precision cannot handle the full exponent range without minor but pessimal changes to renormalize, and although single precision is enough, 25 bit precision is now claimed and used.
Added comments about some of the magic for the error bound 4/6+epsilon. I still don't understand why it is 4/6+ and not 6/6+ ulps.
Indent comments at the right of code more consistently.
|
153447 |
15-Dec-2005 |
bde |
Added comments about the apparently-magic rational function used in the second step of approximating cbrt(x). It turns out to be neither very magic not nor very good. It is just the (2,2) Pade approximation to 1/cbrt(r) at r = 1, arranged in a strange way to use fewer operations at a cost of replacing 4 multiplications by 1 division, which is an especially bad tradeoff on machines where some of the multiplications can be done in parallel. A Remez rational approximation would give at least 2 more bits of accuracy, but the (2,2) Pade approximation already gives 6 more bits than needed. (Changed the comment which essentially says that it gives 3 more bits.)
Lower order Pade approximations are not quite accurate enough for double precision but are plenty for float precision. A lower order Remez rational approximation might be enough for double precision too. However, rational approximations inherently require an extra division, and polynomial approximations work well for 1/cbrt(r) at r = 1, so I plan to switch to using the latter. There are some technical complications that tend to cost a division in another way.
|
153386 |
13-Dec-2005 |
bde |
Optimize by not doing excessive conversions for handling the sign bit. This gives an optimization of between 9 and 22% on Athlons (largest for cbrt() on amd64 -- from 205 to 159 cycles).
We extracted the sign bit and worked with |x|, and restored the sign bit as the last step. We avoided branches to a fault by using accesses to FP values as bits to clear and restore the sign bit. Avoiding branches is usually good, but the bit access macros are not so good (especially for setting FP values), and here they always caused pipeline stalls on Athlons. Even using branches would be faster except on args that give perfect branch misprediction, since only mispredicted branches cause stalls, but it possible to avoid touching the sign bit in FP values at all (except to preserve it in conversions from bits to FP not related to the sign bit). Do this. The results are identical except in 2 of the 3 unsupported rounding modes, since all the approximations use odd rational functions so they work right on strictly negative values, and the special case of -0 doesn't use an approximation.
|
153382 |
13-Dec-2005 |
bde |
Fixed some especially horrible style bugs (indentation that is neither KNF nor fdlibmNF combined with multiple statements per line).
|
153306 |
11-Dec-2005 |
bde |
Added comments about the magic behind <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).
|
153305 |
11-Dec-2005 |
bde |
Fixed the unexpectedly large maximum error after the previous commit. 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).
|
153303 |
11-Dec-2005 |
bde |
Fixed all 502518670 errors of more than 1 ulp for cbrtf() on amd64. 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.)
|
153302 |
11-Dec-2005 |
bde |
Fixed some magic numbers.
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.
|
153178 |
06-Dec-2005 |
obrien |
"Create" ldexpf for non-i386 architectures.
Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu>
|
153088 |
04-Dec-2005 |
bde |
Fixed the approximation to pio4. pio4_hi must be pio2_hi/2 since it 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.
|
153086 |
04-Dec-2005 |
bde |
For log1pf(), fixed the approximations to sqrt(2), sqrt(2)-1 and 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).
|
153078 |
04-Dec-2005 |
bde |
Use the usual volatile hack to trick gcc into clipping any extra precision 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].
|
153049 |
03-Dec-2005 |
bde |
Fixed fdlibm[+cygnus] logbf() and logb() on denormals. Adjustment 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.
|
153046 |
03-Dec-2005 |
bde |
Restored removal of the special handling needed for a result of +-0. 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
|
153043 |
03-Dec-2005 |
bde |
Simplified the fix in rev.1.3. Instead of using long double for 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.
|
153042 |
03-Dec-2005 |
bde |
Fixed rint(x) in the following cases: (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).
|
153017 |
02-Dec-2005 |
bde |
Fixed roundf(). The following cases never worked in FreeBSD: - 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.
|
152951 |
30-Nov-2005 |
bde |
Rearranged the polynomial evaluation to reduce dependencies, as in 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).
|
152949 |
30-Nov-2005 |
bde |
Fixed cosf(x) when x is a "negative" NaNs. I broke this in rev.1.10. 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.
|
152947 |
30-Nov-2005 |
bde |
Fixed the hi+lo approximation to log(2). The normal 17+24 bit decomposition 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.
|
152881 |
28-Nov-2005 |
bde |
Rearranged the polynomial evaluation some more to reduce dependencies. 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).
|
152878 |
28-Nov-2005 |
bde |
Fixed about 50 million errors of infinity ulps and about 3 million errors 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.
|
152872 |
28-Nov-2005 |
bde |
Exploit skew-symmetry to avoid an operation: -sin(x-A) = sin(A-x). This 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).
|
152871 |
28-Nov-2005 |
bde |
Try to use the "proximity" (~) operator consistently in comments (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.
|
152870 |
28-Nov-2005 |
bde |
Changed spelling of the request-to-inline macro name to match the change 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.
|
152869 |
28-Nov-2005 |
bde |
Use only double precision for "kernel" cosf and sinf (except for 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].
|
152766 |
24-Nov-2005 |
bde |
Minor cleanups and optimizations:
- 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.
|
152755 |
24-Nov-2005 |
joel |
s/5.5/6.0/ in HISTORY section.
Discussed with: ru
|
152741 |
24-Nov-2005 |
bde |
Optimized by eliminating the special case for 0.67434 <= |x| < pi/4.
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.
|
152713 |
23-Nov-2005 |
bde |
Use only double precision for "kernel" tanf (except for returning float). 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.
|
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.
|
152647 |
21-Nov-2005 |
bde |
Mess up the "kernel" float trig function .c files with ifdefs so that 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.
|
152641 |
21-Nov-2005 |
bde |
Use double precision to simplify and optimize a long division.
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.
|
152640 |
20-Nov-2005 |
bde |
Restored a cleanup in rev.1.9 tthat was lost in rev.1.10.
|
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.
|
152566 |
18-Nov-2005 |
bde |
Removed an unused declaration which was so old that it wasn't a prototype and thus just broke building at any nonzero WARNS level.
Fixed nearby style bugs.
|
152551 |
17-Nov-2005 |
ru |
-mdoc sweep.
|
152540 |
17-Nov-2005 |
bde |
Minor cleanups:
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.
|
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.
|
152353 |
13-Nov-2005 |
bde |
Fixed some magic numbers.
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.
|
152352 |
13-Nov-2005 |
bde |
Fixed some magic numbers.
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.
|
152343 |
12-Nov-2005 |
bde |
Imoproved comments for the minimax polynomial.
Removed an unused variable.
Fixed some wrong comments and some nearby misformatting.
|
152341 |
12-Nov-2005 |
bde |
Tweaked the minimax polynomial and improved its comments.
|
152340 |
12-Nov-2005 |
bde |
Improved comments for the minimax polynomial.
|
152335 |
12-Nov-2005 |
bde |
As for the float trig functions, use a minimax polynomial that is 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.
|
152284 |
10-Nov-2005 |
bde |
As for __kernel_cosf() and __kernel_sinf(), use a fairly optimal minimax 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.
|
152133 |
06-Nov-2005 |
bde |
Detach k_rem_pio2f.c from the build since it is now unused. It is a libm internal so this shouldn't cause version problems.
|
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.
|
152117 |
06-Nov-2005 |
bde |
The logb() functions are not just ieee754 "test" functions, but are standard in C99 and POSIX.1-2001+. They are also not deprecated, since apart from being standard they can handle special args slightly better than the ilogb() functions.
Move their documentation to ilogb.3. Try to use consistent and improved wording for both sets of functions. All of ieee854, C99 and POSIX have better wording and more details for special args.
Add history for the logb() functions and ilogbl(). Fix history for ilogb().
|
151969 |
02-Nov-2005 |
bde |
Moved the optimization for tiny x from __kernel_tan[f](x) to tan[f](x) 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.
|
151966 |
02-Nov-2005 |
bde |
Updated the comment about the optimization for tiny x (the previous 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.)
|
151963 |
02-Nov-2005 |
bde |
Removed dead code for handling tan[f]() on odd multiples of pi/2. This 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 { --- > {
|
151961 |
02-Nov-2005 |
bde |
Fixed some of the silliness related to rev.1.8. In 1.8, "double" in 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.
|
151880 |
30-Oct-2005 |
bde |
Fixed spelling of remquof() in its prototype.
|
151879 |
30-Oct-2005 |
bde |
Fixed some comments added in rev.1.5.
The log message for 1.5 said that some small (one or two ulp) inaccuracies were fixed, and a comment implied that the critical change is to switch the rounding mode to to-nearest, with a switch of the precision to extended at no extra cost. Actually, the errors are very large (ucbtest finds ones of several hundred ulps), and it is the switch of the precision that is critical.
Another comment was wrong about NaNs being handled sloppily.
|
151865 |
29-Oct-2005 |
bde |
Implement inline functions to give the complex result x+I*y from float 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>
|
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.
|
151796 |
28-Oct-2005 |
bde |
Use fairly optimal minimax polynomials for __kernel_cosf() and __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.
|
151698 |
26-Oct-2005 |
bde |
Use a better algorithm for reducing the error in __kernel_cos[f](). This supersedes the fix for the old algorithm in rev.1.8 of k_cosf.c.
I want this change mainly because it is an optimization. It helps make software cos[f](x) and sin[f](x) faster than the i387 hardware versions for small x. It is also a simplification, and reduces the maximum relative error for cosf() and sinf() on machines like amd64 from about 0.87 ulps to about 0.80 ulps. It was validated for cosf() and sinf() by exhaustive testing. Exhaustive testing is not possible for cos() and sin(), but ucbtest reports a similar reduction for the worst case found by non-exhaustive testing. ucbtest's non-exhaustive testing seems to be good enough to find problems in algorithms but not maximum relative errors when there are spikes. E.g., short runs of it find only 3 ulp error where the i387 hardware cos() has an error of about 2**40 ulps near pi/2.
|
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.
|
151620 |
24-Oct-2005 |
bde |
Moved the optimization for tiny x from __kernel_{cos,sin}[f](x) to {cos_sin}[f](x) so that x doesn't need to be reclassified in the "kernel" functions to determine if it is tiny (it still needs to be reclassified in the cosine case for other reasons that will go away).
This optimization is quite large for exponentially distributed x, since x is tiny for almost half of the domain, but it is a pessimization for uniformally distributed x since it takes a little time for all cases but rarely applies. Arg reduction on exponentially distributed x rarely gives a tiny x unless the reduction is null, so it is best to only do the optimization if the initial x is tiny, which is what this commit arranges. The imediate result is an average optimization of 1.4% relative to the previous version in a case that doesn't favour the optimization (double cos(x) on all float x) and a large pessimization for the relatively unimportant cases of lgamma[f][_r](x) on tiny, negative, exponentially distributed x. The optimization should be recovered for lgamma*() as part of fixing lgamma*()'s low-quality arg reduction.
Fixed various wrong constants for the cutoff for "tiny". For cosine, the cutoff is when x**2/2! == {FLT or DBL}_EPSILON/2. We round down to an integral power of 2 (and for cos() reduce the power by another 1) because the exact cutoff doesn't matter and would take more work to determine. For sine, the exact cutoff is larger due to the ration of terms being x**2/3! instead of x**2/2!, but we use the same cutoff as for cosine. We now use a cutoff of 2**-27 for double precision and 2**-12 for single precision. 2**-27 was used in all cases but was misspelled 2**27 in comments. Wrong and sloppy cutoffs just cause missed optimizations (provided the rounding mode is to nearest -- other modes just aren't supported).
|
151230 |
11-Oct-2005 |
bde |
Fixed range reduction for large multiples of pi/2 on systems with broken assignment to floats (e.g., i386 with gcc -O, but not amd64 or ia64; i386 with gcc -O0 worked accidentally).
Use an unnamed volatile temporary variable to trick gcc -O into clipping extra precision on assignment. It's surprising that only 1 place needed to be changed.
For tanf() on i386 with gcc -O, the bug caused errors > 1 ulp with a density of 2.3% for args larger in magnitude than 128*pi/2, with a maximum error of 1.624 ulps.
After this fix, exhaustive testing shows that range reduction for floats works as intended assuming that it is in within a factor of about 2^16 of working as intended for doubles. It provides >= 8 extra bits of precision for all ranges. On i386:
range max error in double/single ulps extra precision ----- ------------------------------- --------------- 0 to 3*pi/4 0x000d3132 / 0.0016 9+ bits 3*pi/4 to 128*pi/2 0x00160445 / 0.0027 8+ 128*pi/2 to +Inf 0x00000030 / 0.00000009 23+ 128*pi/2 up, -O0 before fix 0x00000030 / 0.00000009 23+ 128*pi/2 up, -O1 before fix 0x10000000 / 0.5 1
The 23+ bits of extra precision for large multiples corresponds to almost perfect reduction to a pair of floats (24 extra would be perfect).
After this fix, the maximum relative error (relative to the corresponding fdlibm double precision function) is < 1 ulp for all basic trig functions on all 2^32 float args on all machines tested:
amd64 ia64 i386-O0 i386-O1 ------ ------ ------ ------ cosf: 0.8681 0.8681 0.7927 0.5650 sinf: 0.8733 0.8610 0.7849 0.5651 tanf: 0.9708 0.9329 0.9329 0.7035
|
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.
|
151182 |
09-Oct-2005 |
bde |
Fix numerous errors of >= 1 ulp for cosf(x) and sinf(x) (1 line) and add a comment about related magic (many lines)).
__kernel_cos[f]() needs a trick to reduce the error to below 1 ulp when |x| >= 0.3 for the range-reduced x. Modulo other bugs, naive code that doesn't use the trick would have an error of >= 1 ulp in about 0.00006% of cases when |x| >= 0.3 for the unreduced x, with a maximum relative error of about 1.03 ulps. Mistransation of the trick from the double precision case resulted in errors in about 0.2% of cases, with a maximum relative error of about 1.3 ulps.
The mistranslation involved not doing implicit masking of the 32-bit float word corresponding to to implicit masking of the lower 32-bit double word by clearing it.
sinf() uses __kernel_cosf() for half of all cases so its errors from this bug are similar. tanf() is not affected.
The error bounds in the above and in my other recent commit messages are for amd64. Extra precision for floats on i386's accidentally masks this bug, but only if k_cosf.c is compiled with -O. Although the extra precision helps here, this is accidental and depends on longstanding gcc precision bugs (not clipping extra precision on assignment...), and the gcc bugs are mostly avoided by compiling without -O. I now develop libm mainly on amd64 systems to simplify error detection and debugging.
|
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.
|
150318 |
19-Sep-2005 |
bde |
Fixed aliasing bugs in TRUNC() by using the fdlibm macros for access to doubles as bits. fdlibm-1.1 had similar aliasing bugs, but these were fixed by NetBSD or Cygnus before a modified version of fdlibm was imported in 1994. TRUNC() is only used by tgamma() and some implementation-detail functions. The aliasing bugs were detected by compiling with gcc -O2 but don't seem to have broken tgamma() on i386's or amd64's. They broke my modified version of tgamma().
Moved the definition of TRUNC() to mathimpl.h so that it can be fixed in one place, although the general version is even slower than necessary because it has to operate on pointers to volatiles to handle its arg sometimes being volatile. Inefficiency of the fdlibm macros slows down libm generally, and tgamma() is a relatively unimportant part of libm. The macros act as if on 32-bit words in memory, so they are hard to optimize to direct actions on 64-bit double registers for (non-i386) machines where this is possible. The optimization is too hard for gcc on amd64's, and declaring variables as volatile makes it impossible.
|
150067 |
12-Sep-2005 |
das |
Add a missing ldexpf() alias for amd64.
Noticed by: bz@, tjr@
|
148297 |
22-Jul-2005 |
kensmith |
Bump the shared library version number of all libraries that have not been bumped since RELENG_5.
Reviewed by: ru Approved by: re (not needed for commit check but in principle...)
|
147446 |
16-Jun-2005 |
ru |
Markup nit.
Approved by: re (blanket)
|
147445 |
16-Jun-2005 |
ru |
Fixed compile warning.
Approved by: re (blanket)
|
147402 |
15-Jun-2005 |
ru |
Assorted markup fixes.
Approved by: re
|
145969 |
06-May-2005 |
deischen |
Prevent these functions from using stack outside of their frame.
Reported by: Marc Olzheim <marcolz at stack dot nl> OK'd by: das
|
145637 |
28-Apr-2005 |
stefanf |
Revert the last change, the conversion from long double to double can raise unwanted underflow exceptions.
Pointed out by: das
|
145399 |
22-Apr-2005 |
stefanf |
Use double additions to raise the inexact exception to work around problems with long double addition on sparc64.
|
145394 |
22-Apr-2005 |
stefanf |
Fix raising the inexact exception (FE_INEXACT) if the result differs from the argument.
Noticed by: das
|
145208 |
17-Apr-2005 |
ache |
Fix truncl.3 MLINKS
|
145171 |
16-Apr-2005 |
das |
More optimized math functions.
|
145170 |
16-Apr-2005 |
das |
Implement truncl() based on floorl().
|
144772 |
08-Apr-2005 |
das |
Add roundl(), lroundl(), and llroundl().
|
144771 |
08-Apr-2005 |
das |
These files should include s_lround.c instead of s_lrint.c. This only matters for efficiency, not for correctness.
|
144770 |
08-Apr-2005 |
das |
Fix a (coincidentally harmless) bug.
|
144690 |
05-Apr-2005 |
das |
Fix a long-standing bug in k_rem_pio2(), which led to large errors when tanf() was called with big arguments close to multiples of pi/2.
Reported by: ucbtest via bde
|
144650 |
05-Apr-2005 |
das |
Build exp2(), exp2f(), and related documentation.
|
144649 |
05-Apr-2005 |
das |
Document exp2() and exp2f(), and make other minor tweaks and updates.
|
144648 |
05-Apr-2005 |
das |
Implement exp2() and exp2f().
|
144091 |
25-Mar-2005 |
das |
Implement and document remquo() and remquof().
|
143780 |
18-Mar-2005 |
das |
Fix the double rounding problem with subnormals, and remove the XXX comments, which no longer apply.
|
143777 |
18-Mar-2005 |
das |
Add missing prototypes for fma() and fmaf(), and remove an inaccurate comment.
|
143769 |
17-Mar-2005 |
das |
Make the fenv.h routines work for programs that use SSE for floating-point arithmetic on i386. Now I'm going to make excuses for why this code is kinda scary:
- To avoid breaking the ABI with 5.3-RELEASE, we can't change sizeof(fenv_t). I stuck the saved mxcsr in some discontiguous reserved bits in the existing structure.
- Attempting to access the mxcsr on older processors results in an illegal instruction exception, so support for SSE must be detected at runtime. (The extra baggage is optimized away if either the application or libm is compiled with -msse{,2}.)
I didn't run tests to ensure that this doesn't SIGILL on older 486's lacking the cpuid instruction or on other processors lacking SSE. Results from running the fenv regression test on these processors would be appreciated. (You'll need to compile the test with -DNO_STRICT_DFL_ENV.) If you have an 80386, or if your processor supports SSE but the kernel didn't enable it, then you're probably out of luck.
Also, I un-inlined some of the functions that grew larger as a result of this change, moving them from fenv.h to fenv.c.
|
143722 |
16-Mar-2005 |
das |
Spell 'fedisableexcept' correctly.
|
143709 |
16-Mar-2005 |
das |
Document feenableexcept(), fedisableexcept(), and fegetexcept().
|
143708 |
16-Mar-2005 |
das |
Replace fegetmask() and fesetmask() with feenableexcept(), fedisableexcept(), and fegetexcept(). These two sets of routines provide the same functionality. I implemented the former as an undocumented internal interface to make the regression test easier to write. However, fe(enable|disable|get)except() is already part of glibc, and I would like to avoid gratuitous differences. The only major flaw in the glibc API is that there's no good way to report errors on processors that don't support all the unmasked exceptions.
|
143264 |
07-Mar-2005 |
das |
Replace strong references with weak references. There's no particularly good reason to do this, except that __strong_reference does type checking, whereas __weak_reference does not. On Alpha, the compiler won't accept a 'long double' parameter in place of a 'double' parameter even thought the two types are identical.
|
143260 |
07-Mar-2005 |
stefanf |
Remove an obsolete sentence from a comment.
|
143230 |
07-Mar-2005 |
das |
- If z is 0, one of x or y is 0, and the other is infinite, raise an invalid exception and return an NaN. - If a long double has 113 bits of precision, implement fma in terms of simple long double arithmetic instead of complicated double arithmetic. - If a long double is the same as a double, alias fma as fmal.
|
143227 |
07-Mar-2005 |
das |
Document scalbnl and scalblnl.
|
143226 |
07-Mar-2005 |
das |
Document nextafterl and nexttoward{,f,l}.
|
143225 |
07-Mar-2005 |
das |
Add nexttoward to the list of implemented functions, and explicitly list the four that are still missing.
|
143224 |
07-Mar-2005 |
das |
Document fmal.
|
143223 |
07-Mar-2005 |
das |
Remove ldexp and ldexpf. The former is in libc, and the latter is identical to scalbnf, which is now aliased as ldexpf. Note that the old implementations made the mistake of setting errno and were the only libm routines to do so.
|
143222 |
07-Mar-2005 |
das |
- Remove s_ldexpf.c (now aliased to scalbn.) - Add nexttoward{,f,l} and nextafterl. On all platforms, nexttowardl is an alias for nextafterl. - Add fmal. - Add man pages for new routines: fmal, nextafterl, nexttoward{,f,l}, scalb{,l}nl.
Note that on platforms where long double is the same as double, we generally just alias the double versions of the routines, since doing so avoids extra work on the source code level and redundant code in the binary. In particular:
ldbl53 ldbl64/113 fmal s_fma.c s_fmal.c ldexpl s_scalbn.c s_scalbnl.c nextafterl s_nextafter.c s_nextafterl.c nexttoward s_nextafter.c s_nexttoward.c nexttowardf s_nexttowardf.c s_nexttowardf.c nexttowardl s_nextafter.c s_nextafterl.c scalbnl s_scalbn.c s_scalbnl.c
|
143221 |
07-Mar-2005 |
das |
- Define FP_FAST_FMA for sparc64, since fma() is now implemented using sparc64's 128-bit long doubles. - Define FP_FAST_FMAL for ia64. - Prototypes for fmal, frexpl, ldexpl, nextafterl, nexttoward{,f,l}, scalblnl, and scalbnl.
|
143220 |
07-Mar-2005 |
das |
Alias scalbn as ldexpl and scalbnl on platforms where long double is the same as double.
|
143219 |
07-Mar-2005 |
das |
- Implement scalblnl. - In scalbln and scalblnf, check the bounds of the second argument. This is probably unnecessary, but strictly speaking, we should report an error if someone tries to compute scalbln(x, INT_MAX + 1ll).
|
143218 |
07-Mar-2005 |
das |
Implement nexttowardf. This is used on both platforms with 11-bit exponents and platforms with 15-bit exponents for long doubles.
|
143217 |
07-Mar-2005 |
das |
Implement nexttoward and nextafterl; the latter is also known as nexttowardl. These are not needed on machines where long doubles look like IEEE-754 doubles, so the implementation only supports the usual long double formats with 15-bit exponents.
Anything bizarre, such as machines where floating-point and integer data have different endianness, will cause problems. This is the case with big endian ia64 according to libc/ia64/_fpmath.h. Please contact me if you managed to get a machine running this way.
|
143216 |
07-Mar-2005 |
das |
- Try harder to trick gcc into not optimizing away statements that are intended to raise underflow and inexact exceptions. - On systems where long double is the same as double, nextafter should be aliased as nexttoward, nexttowardl, and nextafterl.
|
143213 |
07-Mar-2005 |
das |
Implement frexpl.
|
143212 |
07-Mar-2005 |
das |
Alias frexp as frexpl on platforms where a long double is the same as a double.
|
143211 |
07-Mar-2005 |
das |
Implement fmal.
|
143210 |
07-Mar-2005 |
das |
- Define the LDBL_PREC to be the number of significant bits in a long double's mantissa. - Add an assembly version of fmal.
|
143209 |
07-Mar-2005 |
das |
- Define the LDBL_PREC to be the number of significant bits in a long double's mantissa. - Add an assembly version of scalbnl.
|
143208 |
07-Mar-2005 |
das |
Define the LDBL_PREC to be the number of significant bits in a long double's mantissa.
|
143207 |
07-Mar-2005 |
das |
Add an assembly version of fmal.
|
143206 |
07-Mar-2005 |
das |
Add scalbnl, also known as as ldexpl.
|
143205 |
07-Mar-2005 |
das |
Alias scalbnf as ldexpf. The two are identical in binary floating-point formats.
|
143181 |
06-Mar-2005 |
das |
Fix a mistake in the exponent range.
|
143165 |
05-Mar-2005 |
das |
Work around a gcc bug. This fixes feholdexcept() et al. at -O1. Symptoms of the problem included assembler warnings and nondeterministic runtime behavior when a fe*() call that affects the fpsr is closely followed by a float point op.
The bug (at least, I think it's a bug) is that gcc does not insert a break between a volatile asm and a dependent instruction if the volatile asm came from an inlined function. Volatile asms seem to be fine in other circumstances, even without -mvolatile-asm-stop, so perhaps the compiler adds the stop bits before inlining takes place. The problem does not occur at -O0 because inlining is disabled, and it doesn't happen at -O2 because -fschedule-insns2 knows better.
|
142558 |
26-Feb-2005 |
das |
Un-document the non-extant exp10() and exp10f() functions. exp10() was a casualty of the transition away from the VAX.
|
142369 |
24-Feb-2005 |
das |
Revert rev 1.8, which causes small (e.g. 2 ulp) errors for some inputs. The trouble with replacing two floats with a double is that the latter has 6 extra bits of precision, which actually hurts accuracy in many cases. All of the constants are optimal when float arithmetic is used, and would need to be recomputed to do this right.
Noticed by: bde (ucbtest)
|
142182 |
21-Feb-2005 |
das |
Use hardware instructions for sqrt() and sqrtf().
|
142181 |
21-Feb-2005 |
das |
Use double arithmetic instead of simulating it with two floats. This results in a performance gain on the order of 10% for amd64 (sledge), ia64 (pluto1), i386+SSE (Pentium 4), and sparc64 (panther), and a negligible improvement for i386 without SSE. (The i386 port still uses the hardware instruction, though.)
|
142176 |
21-Feb-2005 |
das |
Remove the i387 versions of atan(), atan2(), and atan2f(). They are slower than the MI routines on modern hardware, except for degenerate cases such as the Pentium 4.
PR: 67469
|
142150 |
20-Feb-2005 |
das |
Remove i387 versions of asin() and acos(). Although the hardware instruction was faster on the 486, it's slower than our MD version on modern processors.
Determined by: bde PR: 67469
|
142149 |
20-Feb-2005 |
das |
Remove the float versions of the i387 trig functions obtained from NetBSD. They're buggy, giving particularly for inputs larger in magnitude than 2**63.
Noticed by: bde PR: 67469
|
141302 |
04-Feb-2005 |
das |
Fix a small scripting snafu in the previous revision.
|
141297 |
04-Feb-2005 |
das |
Remove another vestige of support for a non-IEEE libm.
|
141296 |
04-Feb-2005 |
das |
Reduce diffs against vendor source (Sun fdlibm 5.3).
|
141281 |
04-Feb-2005 |
das |
Move machine-dependent crud to its own makefile.
|
141280 |
04-Feb-2005 |
das |
Remove wrappers and other cruft intended to support SVID, mistakes in C90, and other arcana. Most of these features were never fully supported or enabled by default.
Ok: bde, stefanf
|
140949 |
28-Jan-2005 |
ru |
Typo.
|
140948 |
28-Jan-2005 |
ru |
Properly terminate sentence.
|
140890 |
27-Jan-2005 |
das |
- Move the functions presently described in in ieee(3) to their own manpages. They are not very related, so separating them makes it easier to add meaningful cross-references and extend some of the descriptions. - Move the part of math(3) that discusses IEEE 754 to the ieee(3) manpage.
|
140689 |
24-Jan-2005 |
cognet |
Define FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD and _ROUND_MASK to unbreak the build for arm.
|
140685 |
23-Jan-2005 |
das |
Update comment to reflect the code change in the previous revision.
Noticed by: ceri
|
140681 |
23-Jan-2005 |
das |
Many changes, including the following major ones: - Rearrange the list of functions into categories. - Remove the ulps column. It was appropriate for only some of the functions in the list, and correct for even fewer of them. - Add some new paragraphs, and remove some old ones about NaNs that may do more harm than good. - Document precisions other than double-precision.
|
140666 |
23-Jan-2005 |
das |
If x == y, return y, not x. C99 (though not IEEE 754) requires that nextafter(+0.0, -0.0) returns -0.0 and nextafter(-0.0, +0.0) returns +0.0.
|
140609 |
22-Jan-2005 |
das |
Add fma() and fmaf(), which implement a fused multiply-add operation.
|
140505 |
20-Jan-2005 |
ru |
Sort sections.
|
140355 |
16-Jan-2005 |
ru |
Use the \*(If string provided by mdoc(7), to represent infinity.
|
140354 |
16-Jan-2005 |
ru |
Removed redundant .br call.
|
140275 |
15-Jan-2005 |
das |
amd64 assembly versions of sqrt(), lrint(), and llrint() using SSE2.
|
140274 |
15-Jan-2005 |
das |
Most libm routines depend on the rounding mode and/or set exception flags, so they are not pure. Remove the __pure2 annotation from them. I believe that the following routines and their float and long double counterparts are the only ones here that can be __pure2:
copysign is* fabs finite fmax fmin fpclassify ilogb nan signbit
When gcc supports FENV_ACCESS, perhaps there will be a new annotation that allows the other functions to be considered pure when FENV_ACCESS is off.
Discussed with: bde
|
140270 |
15-Jan-2005 |
das |
Braino. Revert rev 1.50.
Pointy hat to: das
|
140269 |
14-Jan-2005 |
das |
Remove numerous references to VAX floating-point and the setting of errno, replacing them with a discussion of IEEE exceptions where appropriate. Cross-reference fenv(3) whenever exceptions are mentioned.
|
140265 |
14-Jan-2005 |
das |
Set math_errhandling to MATH_ERREXCEPT. Now that we have fenv.h, we basically support this, subject to gcc's lack of FENV_ACCESS support. In any case, the previous setting of math_errhandling to 0 is not allowed by POSIX.
|
140263 |
14-Jan-2005 |
das |
Remove some #if 0'd code.
|
140225 |
14-Jan-2005 |
ru |
Tiny markup nits.
|
140219 |
14-Jan-2005 |
das |
Mark all inline asms that read the floating-point control or status registers as volatile. Instructions that *wrote* to FP state were already marked volatile, but apparently gcc has license to move non-volatile asms past volatile asms. This broke amd64's feupdateenv at -O2 due to a WAR conflict between fnstsw and fldenv there.
|
140200 |
13-Jan-2005 |
stefanf |
Fixed too many of "the", and enclose multi-word argument in double quotes.
Obtained from: ru
|
140195 |
13-Jan-2005 |
das |
Import the subset of J.T. Conklin's single-precision x86-optimized math routines that appear to be (a) correct and (b) faster than their MI counterparts on my Pentium 4.
Obtained from: NetBSD
|
140188 |
13-Jan-2005 |
das |
The isnormal() in rev 1.2 should have been isfinite() so subnormals round correctly.
Noticed by: stefanf
|
140187 |
13-Jan-2005 |
das |
Things that are broken, unneeded, and unused since 1997 belong in the attic.
|
140177 |
13-Jan-2005 |
ru |
Markup nits.
|
140174 |
13-Jan-2005 |
ru |
Fixed too many of "the", and enclose multi-word argument in double quotes.
|
140172 |
13-Jan-2005 |
stefanf |
Implement and document ceill().
|
140171 |
13-Jan-2005 |
stefanf |
Bump .Dd for the last commit.
|
140143 |
12-Jan-2005 |
stefanf |
Hook up and document floorl().
|
140142 |
12-Jan-2005 |
stefanf |
Implement floorl().
|
140141 |
12-Jan-2005 |
stefanf |
Whitespace nit.
|
140088 |
11-Jan-2005 |
das |
Add MI implementations of [l]lrint[f]() and [l]lround[f]().
Discussed with: bde
|
140087 |
11-Jan-2005 |
das |
Document [l]lrint[f]() and [l]lround[f]().
|
140086 |
11-Jan-2005 |
das |
Faster lrint() and llrint() implementations for x86.
|
140085 |
11-Jan-2005 |
das |
Mark inline stmxcsr instructions as volatile, since this appears to be the only way to convince gcc that they read the MXCSR. The volatile annotation may be needed elsewhere as well.
|
140081 |
11-Jan-2005 |
ru |
Scheduled mdoc(7) sweep.
|
140079 |
11-Jan-2005 |
ru |
Sanitize the markup, as prompted.
|
138925 |
16-Dec-2004 |
das |
GC unused declaration
|
138924 |
16-Dec-2004 |
das |
Cosmetic changes only: - style - remove unused variables - de-support VAX
Inspired by: bin/42388
|
136395 |
11-Oct-2004 |
das |
More updates for math(3): - Make some minor rearrangements in the introduction. - Mention the problem with argument reduction on i386. - Add recently-implemented functions to the table. - Un-document the error bounds that only apply to the old 4BSD math library, and fill in the correct values where I know them. No attempt has been made to document bounds lower than 1 ulp, although smaller bounds are usually achievable in round-to-nearest mode.
|
136385 |
11-Oct-2004 |
stefanf |
Add and document ilogbl(), a long double version of ilogb().
|
136332 |
09-Oct-2004 |
stefanf |
Use the FP_ILOG macros from <math.h> rather than hardcoded return values. Also be prepared for FP_ILOGBNAN != INT_MAX.
Reviewed by: md5
|
136027 |
01-Oct-2004 |
kensmith |
Bump the library version numbers for the following libraries:
/lib/{libm,libreadline} /usr/lib/{libhistory,libopie,libpcap}
in preparation for doing the same thing to RELENG_5. HUGE amounts of help for determining what to bump provided by kris.
Discussed on: freebsd-current Approved by: re (not required for commit but something like this should be)
|
135360 |
17-Sep-2004 |
das |
Further refine some #ifs: - Simplify the logic by using __GNUC_PREREQ__. Suggested by stefanf. - Make math.h compile with old (pre-8.0) versions of icc. Submitted by sf [sic].
|
133289 |
07-Aug-2004 |
stefanf |
Add man pages for the cimag(), conj() and creal() functions.
|
133174 |
05-Aug-2004 |
cognet |
Only use rfs and wfs if ARM_HARD_FLOAT is defined, and use stubs if it is not, in order to unbreak arm make world. The right way to do it with soft floats will be figured out later. Discussed with: das
|
133147 |
05-Aug-2004 |
das |
Replace s_isnan.c and s_isnanf.c with the more compact s_isnan.c from libc. The externally-visible effect of this is to add __isnanl() to libm, which means that libm.so.2 can once again link against libc.so.4 when LD_BIND_NOW is set. This was broken by the addition of fdiml(), which calls __isnanl().
|
133146 |
05-Aug-2004 |
das |
Use isnormal() instead of fpclassify() to avoid dependency on libc.so.5.
|
132760 |
28-Jul-2004 |
kan |
Work around known GCC 3.4.x problem and use ANSI prototype for dremf().
|
132382 |
19-Jul-2004 |
das |
Fix two bugs in the signbit() macro, which was implemented last year:
- It was added to libc instead of libm. Hopefully no programs rely on this mistake.
- It didn't work properly on large long doubles because its argument was converted to type double, resulting in undefined behavior.
|
132292 |
17-Jul-2004 |
stefanf |
Fix minor namespace pollution: The prototypes for f{dim,max,min}(), nearbyint(), round() and trunc() shouldn't be visible when compiling with -D_XOPEN_SOURCE=500.
|
132246 |
16-Jul-2004 |
das |
Tweak the conditions under which certain gcc builtins are used:
- Unlike the builtin relational operators, builtin floating-point constants were not available until gcc 3.3, so account for this.[1]
- Apparently some versions of the Intel C Compiler fallaciously define __GNUC__ without actually being compatible with the claimed gcc version. Account for this, too.[2]
[1] Noticed by: Christian Hiris <4711@chello.at> [2] Submitted by: Alexander Leidinger <Alexander@Leidinger.net>
|
131865 |
09-Jul-2004 |
das |
Remove the declaration of isnan() from this file. It is no longer needed as of math.h v1.40, and its prototype is incorrect here.
|
131852 |
09-Jul-2004 |
das |
Implement the classification macros isfinite(), isinf(), isnan(), and isnormal() the hard way, rather than relying on fpclassify(). This is a lose in the sense that we need a total of 12 functions, but it is necessary for binary compatibility because we have never bumped libm's major version number. In particular, isinf(), isnan(), and isnanf() were BSD libc functions before they were C99 macros, so we can't reimplement them in terms of fpclassify() without adding a dependency on libc.so.5. I have tried to arrange things so that programs that could be compiled in FreeBSD 4.X will generate the same external references when compiled in 5.X. At the same time, the new macros should remain C99-compliant.
The isinf() and isnan() functions remain in libc for historical reasons; however, I have moved the functions that implement the macros isfinite() and isnormal() to libm where they belong. Moreover, half a dozen MD versions of isinf() and isnan() have been replaced with MI versions that work equally well.
Prodded by: kris
|
131851 |
09-Jul-2004 |
das |
Define the following macros in terms of [gi]cc builtins when the builtins are available: HUGE_VAL, HUGE_VALF, HUGE_VALL, INFINITY, and NAN. These macros now expand to floating-point constant expressions rather than external references, as required by C99. Other compilers will retain the historical behavior. Note that it is not possible say, e.g. #define HUGE_VAL 1.0e9999 because the above may result in diagnostics at translation time and spurious exceptions at runtime. Hence the need for compiler support for these features.
Also use builtins to implement the macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), and isunordered() when such builtins are available. Although the old macros are correct, the builtin versions are much faster, and they avoid double-expansion problems.
|
131676 |
06-Jul-2004 |
das |
Add C99's nearbyint{,f}() functions as wrappers around rint(). These trivial implementations are about 25 times slower than rint{,f}() on x86 due to the FP environment save/restore. They should eventually be redone in terms of fegetround() and bit fiddling.
|
131539 |
03-Jul-2004 |
ru |
Eliminate double whitespace.
|
131504 |
02-Jul-2004 |
ru |
Mechanically kill hard sentence breaks.
|
131421 |
01-Jul-2004 |
ru |
Markup, grammar, punctuation.
|
131320 |
30-Jun-2004 |
das |
Implement and document fdim{,f,l}, fmax{,f,l}, and fmin{,f,l}.
|
131001 |
24-Jun-2004 |
marcel |
s/ARCH/ARCH_SUBDIR/g -- This reduces the chance of possible conflicts with the user's environment.
Wondered why his cross-builds kept failing: marcel
|
130775 |
20-Jun-2004 |
stefanf |
Completely remove s_ilogb.S as the assembler implementation gives very little speed improvement to none at all over the MI version.
Submitted by: bde
|
130774 |
20-Jun-2004 |
das |
Uncomment some functions that we now support.
|
130770 |
20-Jun-2004 |
das |
Cross-reference round(3) and trunc(3) as appropriate.
|
130769 |
20-Jun-2004 |
das |
Connect scalbln(), trunc(), and the associated documentation to the build.
|
130768 |
20-Jun-2004 |
das |
Declare scalbln(), scalblnf(), trunc(), and truncf().
|
130767 |
20-Jun-2004 |
das |
Implement trunc() and truncf().
|
130766 |
20-Jun-2004 |
das |
Add trivial implementations of scalbln() and scalblnf(). These routines are specified in C99 for the sake of architectures where an int isn't big enough to represent the full range of floating-point exponents. However, even the 128-bit long double format has an exponent smaller than 15 bits, so for all practical purposes, scalbln() and scalblnf() are aliases for scalbn() and scalbnf(), respectively.
|
130715 |
19-Jun-2004 |
stefanf |
Document ilogb()'s return values in terms of the FP_ILOGB* macros.
|
130714 |
19-Jun-2004 |
stefanf |
Return the same result as the MI version for 0.0, INFINITY and NaN.
Reviewed by: standards@
|
130713 |
19-Jun-2004 |
stefanf |
Our MI implementation of ilogb() returns -INT_MAX for the argument 0.0 rather than INT_MIN, so adjust FP_ILOGB0 to reflect this. Use <machine/_limits.h> for INT_MAX's value while there.
Reviewed by: standards@
|
130706 |
19-Jun-2004 |
das |
Memory's free, but all the world ain't a VAX anymore. Bring math.3 kicking and screaming into the 1980's. This change converts most of the markup from man(7) to mdoc(7) format, and I believe it removes or updates everything that was flat out wrong. However, much work is still needed to sanitize the markup, improve coverage, and reduce overlap with other manpages. Some of the sections would better belong in a philosophy_of_w_kahan.3 manpage, but they are informative and remain at least as reminders of topics to cover.
Reviewed by: doc@, trhodes@
|
130371 |
12-Jun-2004 |
das |
The references to scalbn and scalbnf should be scalb and scalbf. (The former are actually useful, and ieee_test(3) only documents functions that aren't.) Add a sentence describing the domain of scalb() and scalbf().
|
130329 |
11-Jun-2004 |
das |
Shift the FPSR contents by the correct amount so feupdateenv() raises the correct exceptions from the old environment.
|
130328 |
11-Jun-2004 |
das |
Insert a missing '~' in feholdexcept(), so that it correctly clears the exception flags in the mxcsr as well as the x87 FPU.
|
130285 |
09-Jun-2004 |
das |
Fix a bug where rintf() rounded the wrong way in round-to-nearest mode on all inputs of the form x.75, where x is an even integer and log2(x) = 21. A similar problem occurred when rounding upward. The bug involves the following snippet copied from rint():
i>>=1; if((i0&i)!=0) i0 = (i0&(~i))|((0x100000)>>j0);
The constant 0x100000 should be 0x200000. Apparently this case was never tested.
It turns out that the bit manipulation is completely superfluous anyway, so remove it. (It tries to simulate 90% of the rounding process that the FPU does anyway.) Also, the special case of +-0 is handled twice (in different ways), so remove the second instance.
Throw in some related simplifications from bde:
- Work around a bug where gcc fails to clip to float precision by declaring two float variables as volatile. Previously, we tricked gcc into generating correct code by declaring some float constants as doubles.
- Remove additional superfluous bit manipulation.
- Minor reorganization.
- Include <sys/types.h> explicitly.
Note that some of the equivalent lines in rint() also appear to be unnecessary, but I'll defer to the numerical analysts who wrote it, since I can't test all 2^64 cases.
Discussed with: bde
|
130264 |
09-Jun-2004 |
das |
Include <sys/cdefs.h> earlier to get the various visibility constants. Previously, we were relying on <sys/_types.h> to include it implicitly.
|
130179 |
07-Jun-2004 |
das |
Add round(3) and roundf(3) and the associated documentation.
PR: 59797 Submitted by: "Steven G. Kargl" <kargl@troutmask.apl.washington.edu> Reviewed by: bde (earlier version, last year)
|
130149 |
06-Jun-2004 |
das |
Add fenv.h, fenv.c, and the associated documentation to the libm build. To facilitate this, add ${.CURDIR}/${ARCH} to make's search path unconditionally.
Reviewed by: standards@
|
130148 |
06-Jun-2004 |
das |
Add documentation for: - fenv(3) - feclearexcept(3), fegetexceptflag(3), feraiseexcept(3), fesetexceptflag(3), fetestexcept(3) - fegetround(3), fesetround(3) - fegetenv(3), feholdexcept(3), fesetenv(3), feupdateenv(3)
Reviewed by: standards@
|
130147 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the sparc64 port.
Reviewed by: standards@
|
130146 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the powerpc port.
Reviewed by: standards@
|
130145 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the ia64 port.
Reviewed by: standards@
|
130144 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the i386 port.
Reviewed by: standards@
|
130143 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the arm port.
It does not appear to be possible to cross-build arm from i386 at the moment, and I have no ARM hardware anyway. Thus, I'm sure there are bugs. I will gladly fix these when the arm port is more mature.
Reviewed by: standards@
|
130142 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the amd64 port.
Reviewed by: standards@
|
130141 |
06-Jun-2004 |
das |
Add an fenv.h implementation for the alpha port. All of the standard features appear to work, subject to the caveat that you tell gcc you want standard rather than recklessly fast behavior (-mieee-with-inexact -mfp-rounding-mode=d).
The non-standard feature of delivering a SIGFPE when an application raises an unmasked exception does not work, presumably due to a kernel bug. This isn't so bad given that floating-point exceptions on the Alpha architecture are not precise, so making them useful in userland requires a significant amount of wizardry.
Reviewed by: standards@
|
130003 |
02-Jun-2004 |
bde |
Fixed lots of 1 ULP errors caused by a broken approximation for pi/2. We approximate pi with more than float precision using pi_hi+pi_lo in the usual way (pi_hi is actually spelled pi in the source code), and expect (float)0.5*pi_lo to give the low part of the corresponding approximation for pi/2. However, the high part for pi/2 (pi_o_2) is rounded to nearest, which happens to round up, while the high part for pi was rounded down. Thus pi_o_2+(float)0.5*pi (in infinite precision) was a very bad approximation for pi/2 -- the low term has the wrong sign and increases the error drom less than half an ULP to a full ULP.
This fix rounds up instead of down for pi_hi. Consistently rounding down instead of up should work, and is the method used in e_acosf.c and e_asinf.c. The reason for the difference is that we sometimes want to return precisely pi/2 in e_atan2f.c, so it is convenient to have a correctly rounded (to nearest) value for pi/2 in a variable. a_acosf.c and e_asinf.c also differ in directly approximating pi/2 instead pi; they multiply by 2.0 instead of dividing by 0.5 to convert the approximation.
These complications are not directly visible in the double precision versions because rounding to nearest happens to round down.
|
129981 |
02-Jun-2004 |
das |
Port a bugfix from FDLIBM 5.3. The bug really only applies to tan() and not tanf() because float type can't represent numbers large enough to trigger the problem. However, there seems to be a precedent that the float versions of the fdlibm routines should mirror their double counterparts.
Also update to the FDLIBM 5.3 license.
Obtained from: FDLIBM Reviewed by: exhaustive comparison
|
129980 |
02-Jun-2004 |
das |
Merge a bugfix from FDLIBM 5.3 to ensure that the error in tan() is always less than 1 ulp. Also update to the 5.3 license.
Obtained from: FDLIBM
|
129959 |
01-Jun-2004 |
bde |
Merged from double precision case (e_pow.c 1.10: sign fixes).
|
129956 |
01-Jun-2004 |
bde |
Fixed the sign of the result in some overflow and underflow cases (ones where the exponent is an odd integer and the base is negative).
Obtained from: fdlibm-5.3
Sun finally released a new version of fdlibm just a coupe of weeks ago. It only fixes 3 bugs (this one, another one in pow() that we already have (rev.1.9), and one in tan(). I've learned too much about powf() lately, so this fix was easy to merge. The patch is not verbatim, because our base version has many differences for portability and I didn't like global renaming of an unrelated variable to keep it separate from the sign variable. This patch uses a new variable named sn for the sign.
|
129955 |
01-Jun-2004 |
bde |
Fixed another precision bug in powf(). This one is in the computation [t=p_l+p_h High]. We multiply t by lg2_h, and want the result to be exact. For the bogus float case of the high-low decomposition trick, we normally discard the lowest 12 bits of the fraction for the high part, keeping 12 bits of precision. That was used for t here, but it doesnt't work because for some reason we only discard the lowest 9 bits in the fraction for lg2_h. Discard another 3 bits of the fraction for t to compensate.
This bug gave wrong results like:
powf(0.9999999, -2.9999995) = 1.0000002 (should be 1.0000001) hex values: 3F7FFFFF C03FFFFE 3F800002 3F800001
As explained in the log for the previous commit, the bug is normally masked by doing float calculations in extra precision on i386's, but is easily detected by ucbtest on systems that don't have accidental extra precision.
This completes fixing all the bugs in powf() that were routinely found by ucbtest.
|
129950 |
01-Jun-2004 |
bde |
Fixed 2 bugs in the computation /* t_h=ax+bp[k] High */. (1) The bit for the 1.0 part of bp[k] was right shifted by 4. This seems to have been caused by a typo in converting e_pow.c to e_powf.c. (2) The lower 12 bits of ax+bp[k] were not discarded, so t_h was actually plain ax+bp[k]. This seems to have been caused by a logic error in the conversion.
These bugs gave wrong results like:
powf(-1.1, 101.0) = -15158.703 (should be -15158.707) hex values: BF8CCCCD 42CA0000 C66CDAD0 C66CDAD4
Fixing (1) gives a result wrong in the opposite direction (hex C66CDAD8), and fixing (2) gives the correct result.
ucbtest has been reporting this particular wrong result on i386 systems with unpatched libraries for 9 years. I finally figured out the extent of the bugs. On i386's they are normally hidden by extra precision. We use the trick of representing floats as a sum of 2 floats (one much smaller) to get extra precision in intermediate calculations without explicitly using more than float precision. This trick is just a pessimization when extra precision is available naturally (as it always is when dealing with IEEE single precision, so the float precision part of the library is mostly misimplemented). (1) and (2) break the trick in different ways, except on i386's it turns out that the intermediate calculations are done in enough precision to mask both the bugs and the limited precision of the float variables (as far as ucbtest can check).
ucbtest detects the bugs because it forces float precision, but this is not a normal mode of operation so the bug normally has little effect on i386's.
On systems that do float arithmetic in float precision, e.g., amd64's, there is no accidental extra precision and the bugs just give wrong results.
|
129864 |
30-May-2004 |
stefanf |
Add implementations for cimag{,f,l}, creal{,f,l} and conj{,f,l}. They are needed for cases where GCC's builtin functions cannot be used and for compilers that don't know about them.
Approved by: das (mentor)
|
129312 |
17-May-2004 |
das |
Remove some kludges designed to ensure that the compiler didn't round constants the wrong way on the VAX. Instead, use C99 hexadecimal floating-point constants, which are guaranteed to be exact on binary IEEE machines. (The correct hexadecimal values were already provided in the source, but not used.) Also, convert the constants to lowercase to work around a gcc bug that wasn't fixed until gcc 3.4.0.
Prompted by: stefanf
|
129040 |
07-May-2004 |
stefanf |
Add an implementation of copysignl(), a long double version of copysign().
Approved by: das (mentor)
|
129038 |
07-May-2004 |
stefanf |
Add an MLINK for fabsl().
Approved by: das (mentor)
|
128999 |
06-May-2004 |
stefanf |
The prototypes for cabs() and cabsf() are in <complex.h>. Fix their arguments' types and describe them briefly.
Reviewed by: ru, bde Approved by: das (mentor)
|
128628 |
25-Apr-2004 |
das |
Make sure that symbols are declared in math.h iff the appropriate namespaces are visible. Previously, math.h failed to hide some C99-, XSI-, and BSD-specific symbols in certain compilation environments.
The referenced PR has a nice listing of the appropriate conditions for making symbols visible in math.h. The only non-stylistic difference between the patch in the PR and this commit is that I superfluously test for __BSD_VISIBLE in a few places to be more explicit about which symbols have historically been part of the FreeBSD environment.
PR: 65939 Submitted by: Stefan Farfeleder <stefan@fafoe.narf.at>
|
128627 |
25-Apr-2004 |
das |
Remove a stale comment referring to values.h, which has never been part of FreeBSD.
PR: 65939
|
126871 |
12-Mar-2004 |
bde |
Initial support for C99's (or is it POSIX.1-2001's?) MATH_ERRNO, MATH_ERREXCEPTION and math_errhandling, so that C99 applications at least have the possibility of determining that errno is not set for math functions. Set math_errhandling to the non-standard-conforming value of 0 for now to indicate that we don't support either method of reporting errors. We intentionally don't support MATH_ERRNO because errno is a mistake, and we are missing support for MATH_ERREXCEPTION (<fenv.h>, compiler support for <fenv.h>, and actually setting the exception flags correctly).
|
121590 |
27-Oct-2003 |
das |
Fix a problem where libm compiled under 5.X would depend on features that are only in libc.so.5. This broke some 4.X applications linked to libm and run under 5.X.
Background: In C99, isinf() and isnan() cannot be implemented as regular functions. We use macros that call libc functions in 5.X, but for libm-internal use, we need to use the old versions until the next time libm's major version number is bumped.
Submitted by: bde Reported by: imp, kris
|
121513 |
25-Oct-2003 |
des |
Better safe than clever.
Submitted by: das
|
121502 |
25-Oct-2003 |
des |
Document fabsl(3).
Submitted by: Stefan Farfeleder <stefan@fafoe.narf.at>
|
121497 |
25-Oct-2003 |
des |
- fabsl.c should be named s_fabsl.c for consistency with libmsun's documented naming scheme (unfortunately the documentation isn't in the tree as far as I can tell); no repocopy is required as there is no history to preserve.
- replace simple and almost-correct implementation with slightly hackish but definitely correct implementation (tested on i386, alpha, sparc64) which requires pulling in fpmath.h and the MD _fpmath.h from libc.
- try not to make a mess of the Makefile in the process.
- enterprising minds are encouraged to implement more C99 long double functions.
|
121419 |
23-Oct-2003 |
des |
Connect fabsl.c to the build.
|
121418 |
23-Oct-2003 |
des |
Add prototypes for all long double functions in C99. Leave them all #if 0'd out, except for fabsl(3) which I've implemented.
|
121417 |
23-Oct-2003 |
des |
Implement fabsl(3), allowing the world to build with -fno-builtin.
|
119017 |
17-Aug-2003 |
gordon |
Stage 3 of dynamic root support. Make all the libraries needed to run binaries in /bin and /sbin installed in /lib. Only the versioned files reside in /lib, the .so symlink continues to live /usr/lib so the toolchain doesn't need to be modified.
|
117917 |
23-Jul-2003 |
bde |
Fixed some style bugs (misplacement and misformatting of some commented-out code).
|
117912 |
23-Jul-2003 |
peter |
Only provide one copy of the math functions. If we provide a MD function, do not also provide a __generic_XXX version as well. This is how we used to runtime select the generic vs i387 versions on the i386 platform.
This saves a pile of #defines in the src/math_private.h file to undo the __generic_XXX renames in some of the *.c files.
|
117910 |
23-Jul-2003 |
peter |
No longer need the internal __get_hw_float() function.
|
117909 |
23-Jul-2003 |
peter |
Now that we do not need to do runtime detection for the broken default fp emulator, stop doing the runtime selection of hardware or emulated floating point operations on i386. Note that I have not suppressed the duplicate compiles yet.
While here, fix the alpha. It has provided specific copysign/copysignf functions since the beginning of time, but they have never been used.
|
115230 |
22-May-2003 |
mike |
Fix two misuses of __BSD_VISIBLE.
Submitted by: bde Approved by: re
|
114331 |
30-Apr-2003 |
peter |
AMD64 support (another IEEEFP platform)
|
113077 |
04-Apr-2003 |
das |
Fix braino in definition of isfinite().
Noticed by: marcus Pointy hat to: das
|
111762 |
02-Mar-2003 |
ru |
mdoc(7) police: Nits.
|
111546 |
26-Feb-2003 |
imp |
- gamma_r, lgamma_r, gammaf_r, and lgammaf_r were protected by _REENTRANT in math.h; the consensus here was that __BSD_VISIBLE was correct instead.
- gamma_r, lgamma_r, gammaf_r, and lgammaf_r had no documentation in the lgamma(3) manpage.
Reviewed by: standards@ Submitted by: Ben Mesander
|
110769 |
12-Feb-2003 |
mike |
o Implement C99 classification macros isfinite(), isinf(), isnan(), isnormal(). The current isinf() and isnan() are perserved for binary compatibility with 5.0, but new programs will use the macros. o Implement C99 comparison macros isgreater(), isgreaterequal(), isless(), islessequal(), islessgreater(), isunordered().
Submitted by: David Schultz <dschultz@uclink.Berkeley.EDU>
|
110734 |
11-Feb-2003 |
mike |
Implement C99's signbit() macro.
|
110566 |
08-Feb-2003 |
mike |
Implement fpclassify(): o Add a MD header private to libc called _fpmath.h; this header contains bitfield layouts of MD floating-point types. o Add a MI header private to libc called fpmath.h; this header contains bitfield layouts of MI floating-point types. o Add private libc variables to lib/libc/$arch/gen/infinity.c for storing NaN values. o Add __double_t and __float_t to <machine/_types.h>, and provide double_t and float_t typedefs in <math.h>. o Add some C99 manifest constants (FP_ILOGB0, FP_ILOGBNAN, HUGE_VALF, HUGE_VALL, INFINITY, NAN, and return values for fpclassify()) to <math.h> and others (FLT_EVAL_METHOD, DECIMAL_DIG) to <float.h> via <machine/float.h>. o Add C99 macro fpclassify() which calls __fpclassify{d,f,l}() based on the size of its argument. __fpclassifyl() is never called on alpha because (sizeof(long double) == sizeof(double)), which is good since __fpclassifyl() can't deal with such a small `long double'.
This was developed by David Schultz and myself with input from bde and fenner.
PR: 23103 Submitted by: David Schultz <dschultz@uclink.Berkeley.EDU> (significant portions) Reviewed by: bde, fenner (earlier versions)
|
108533 |
01-Jan-2003 |
schweikh |
Correct typos, mostly s/ a / an / where appropriate. Some whitespace cleanup, especially in troff files.
|
108317 |
27-Dec-2002 |
schweikh |
english(4) police.
|
106268 |
31-Oct-2002 |
archie |
Re-apply the previously backed-out commit that fixes the problem where HUGE_VAL is not properly aligned on some architectures. The previous fix now works because the two versions of 'math.h' (include/math.h and lib/msun/src/math.h) have since been merged into one.
PR: bin/43544
|
105809 |
23-Oct-2002 |
markm |
Remove duplicate declaration.
|
104281 |
01-Oct-2002 |
bde |
Fixed a last-minute editing error in previous commit. nfs and/or cvs replaced a 14-byte change in the middle of the file with 14 NULs at EOF despite or because of aborting the initial commit to pick up the change.
|
104280 |
01-Oct-2002 |
bde |
Merged all interesting difference between the old math.h and the current one into the latter and removed the former.
This works around the bug that some broken Makefiles add -I.../src/include to CFLAGS, resulting in the old math.h being preferred and differences between the headers possibly being fatal.
The merge mainly involves declaring some functions as __pure2 although they are not yet all strictly free of side effects.
PR: 43544
|
103686 |
20-Sep-2002 |
archie |
Revert previous commit to unbreak world until we figure out the right way to do it.
|
103653 |
19-Sep-2002 |
archie |
Fix a problem with the definition of HUGE_VAL causing the gcc warning "cast increases required alignment of target type" on some platforms.
Reviewed by: bde
|
98349 |
17-Jun-2002 |
bde |
e_pow.c: Fixed pow(x, y) when x is very close to -1.0 and y is a very large odd integer. E.g., pow(-1.0 - pow(2.0, -52.0), 1.0 + pow(2.0, 52.0)) was 0.0 instead of being very close to -exp(1.0).
PR: 39236 Submitted by: Stephen L Moshier <steve@moshier.net>
e_powf.c: Apply the same patch although it is just cosmetic because odd integers large enough to cause the problem are too large to be precisely represented as floats.
MFC after: 1 week
|
97413 |
28-May-2002 |
alfred |
Fix formatting, this is hard to explain, so I'll show one example.
- float ynf(int n, float x) /* wrapper ynf */ +float +ynf(int n, float x) /* wrapper ynf */
This is because the __STDC__ stuff was indented.
Reviewed by: md5
|
97409 |
28-May-2002 |
alfred |
Assume __STDC__, remove non-__STDC__ code.
Reviewed by: md5
|
97407 |
28-May-2002 |
alfred |
Assume __STDC__, remove non-__STDC__ code.
Submitted by: keramida
|
97045 |
21-May-2002 |
benno |
Spread the word of PowerPC.
|
96462 |
12-May-2002 |
ru |
Added new bsd.incs.mk which handles installing of header files via INCS. Implemented INCSLINKS (equivalent to SYMLINKS) to handle symlinking include files. Allow for multiple groups of include files to be installed, with the powerful INCSGROUPS knob. Documentation to follow.
Added standard `includes' and `incsinstall' targets, use them in Makefile.inc1. Headers from the following makefiles were not installed before (during `includes' in Makefile.inc1):
kerberos5/lib/libtelnet/Makefile lib/libbz2/Makefile lib/libdevinfo/Makefile lib/libform/Makefile lib/libisc/Makefile lib/libmenu/Makefile lib/libmilter/Makefile lib/libpanel/Makefile
Replaced all `beforeinstall' targets for installing includes with the INCS stuff.
Renamed INCDIR to INCSDIR, for consistency with FILES and SCRIPTS, and for compatibility with NetBSD. Similarly for INCOWN, INCGRP, and INCMODE.
Consistently use INCLUDEDIR instead of /usr/include.
gnu/lib/libstdc++/Makefile and gnu/lib/libsupc++/Makefile changes were only lightly tested due to the missing contrib/libstdc++-v3. I fully tested the pre-WIP_GCC31 version of this patch with the contrib/libstdc++.295 stuff.
These changes have been tested on i386 with the -DNO_WERROR "make world" and "make release".
|
93211 |
26-Mar-2002 |
bde |
Resurrect Lite1's gamma() as C99's tgamma(). Minimal changes.
|
93204 |
26-Mar-2002 |
bde |
Fixed some bugs in the description of plain gamma() (and gammaf()). Give a more detailed and correct history of when gamma() was actually the gamma function.
|
93197 |
26-Mar-2002 |
bde |
Fixed some minor style bugs.
|
92917 |
21-Mar-2002 |
obrien |
Remove __P() usage.
|
92887 |
21-Mar-2002 |
obrien |
Fix SCM ID's.
|
91514 |
01-Mar-2002 |
obrien |
We need an frexp() function.
|
88801 |
02-Jan-2002 |
jake |
Add ifdef sparc64.
|
87805 |
13-Dec-2001 |
phantom |
Fix style bugs (mostly remove 'extern' from function prototypes)
Inspired by: conversation with bde
|
87804 |
13-Dec-2001 |
phantom |
* remove reference to m68k-dependent sources * fix comment
|
86712 |
21-Nov-2001 |
ru |
Grammar nit.
|
86680 |
20-Nov-2001 |
ru |
mdoc(7) police: fixed bugs from rev. 1.15.
|
86077 |
05-Nov-2001 |
dwmalone |
gamma(x) actually returns \log(|\Gamma(x)|), so correct the man page and add an historical note explaining this. This patch is based on Stephen's.
We still need someone to implement tgamma.
PR: 28972, 31764 Submitted by: Stephen Montgomery-Smith <stephen@math.missouri.edu>
|
84980 |
15-Oct-2001 |
dd |
Match parenthesis and don't give names to return values.
PR: 31214
|
84885 |
13-Oct-2001 |
bde |
Fixed missing quoting of >= (in ceil.3) and <= (in floor.3) by reverting to describing these operators in English. This completes the fix in rev.1.3 (rev.1.2 got this wrong by describing wrong operators in English).
Fixed bitrot and improved English in the DESCRIPTION section.
|
84882 |
13-Oct-2001 |
bde |
Fixed missing quoting of [-1, +1].
Submitted by: phantom
|
84881 |
13-Oct-2001 |
bde |
Use ".Lb libm" where it will have an effect (not just in the zombie man pages in libm).
Submitted by: phantom
|
84662 |
08-Oct-2001 |
dfr |
Port to ia64. Actually, just do like the alpha.
|
84403 |
03-Oct-2001 |
bde |
Don't install manpage links for the nonexistent functions exp2(), exp2f(), log2() and log2f().
|
84402 |
03-Oct-2001 |
bde |
Removed .Nm's for the nonexistent functions exp2() and exp2f().
|
84306 |
01-Oct-2001 |
ru |
mdoc(7) police: Use the new .In macro for #include statements.
|
84210 |
30-Sep-2001 |
dillon |
Add __FBSDID()s to libm
|
81462 |
10-Aug-2001 |
ru |
mdoc(7) police: join split punctuation to macro calls.
|
81351 |
09-Aug-2001 |
yar |
Tiny markup fix: `to' isn't a variable
|
79754 |
15-Jul-2001 |
dd |
Remove whitespace at EOL.
|
79531 |
10-Jul-2001 |
ru |
mdoc(7) police: removed HISTORY info from the .Os call.
|
78172 |
13-Jun-2001 |
ru |
Added skeleton <complex.h> (aligned with the POSIX.1-200x), mostly to fix the "-nostdinc WARNS=X" breakage caused by broken prototypes for cabs() and cabsl() in <math.h>.
Reimplemented cabs() and cabsl() using new complex numbers types and moved prototypes from <math.h> to <complex.h>.
|
75670 |
18-Apr-2001 |
ru |
mdoc(7) police: normalize .Nd.
|
74870 |
27-Mar-2001 |
ru |
MAN[1-9] -> MAN.
|
74804 |
26-Mar-2001 |
ru |
Don't use MANDEPEND and MANSRC.
|
73088 |
26-Feb-2001 |
ru |
.St -ansiC -> .St -isoC
|
72126 |
07-Feb-2001 |
ru |
mdoc(7) police: Change -filled displays (which just happen to be the same as -ragged in the current implementation) to -ragged. With mdocNG, -filled displays produce the correct output, formatted and justified to both margins.
|
71895 |
01-Feb-2001 |
ru |
mdoc(7) police: split punctuation characters + misc fixes.
|
70481 |
29-Dec-2000 |
ru |
Prepare for mdoc(7)NG.
|
69857 |
11-Dec-2000 |
ru |
mdoc(7) police: use canonical form of .Dd macro.
|
69051 |
22-Nov-2000 |
ru |
mdoc(7) police: Er macro usage cleanup.
|
68946 |
20-Nov-2000 |
ru |
mdoc(7) police: Nm -> Fn where appropriate.
|
67166 |
15-Oct-2000 |
brian |
Fix #include order
Spotted by: imura
|
61335 |
06-Jun-2000 |
bde |
Removed bogus 'l' suffixes in FP register to register instructions.
|
58647 |
27-Mar-2000 |
obrien |
MFS: Add a "magic" comment to help fixincludes realize it doesn't need to modify this file when building GCC 2.96 [by hand or via the port].
Submitted by: Zack Weinberg <zack@wolery.cumb.org>
|
57695 |
02-Mar-2000 |
sheldonh |
Remove more single-space hard sentence breaks.
|
57686 |
02-Mar-2000 |
sheldonh |
Remove single-space hard sentence breaks. These degrade the quality of the typeset output, tend to make diffs harder to read and provide bad examples for new-comers to mdoc.
|
53036 |
09-Nov-1999 |
phantom |
style fix
PR: docs/14737 Submitted by: Norihiro Kumagai <kuma@nk.rim.or.jp>
|
50476 |
28-Aug-1999 |
peter |
$Id$ -> $FreeBSD$
|
42044 |
24-Dec-1998 |
dfr |
Disable building with alpha software completion options until we upgrade compilers.
|
42029 |
23-Dec-1998 |
dfr |
Implement fpsetmask() and other fp*() functions. Programs should use
#include <ieeefp.h>
to access these functions instead of the i386 specific
#include <machine/floatingpoint.h>
Submitted by: Hidetoshi Shimokawa <simokawa@sat.t.u-tokyo.ac.jp>
|
35925 |
10-May-1998 |
jb |
There is no alpha asm code like on i386, so all the functions that the i386 builds with a __generic prefix need to have that stripped.
|
35399 |
23-Apr-1998 |
pst |
Back out last change
|
35380 |
22-Apr-1998 |
pst |
Fix cabs and cabsf definitions to be prototypes.
|
33662 |
20-Feb-1998 |
jb |
Add alpha support. m68k crept in too. Oops. 8-)
|
33107 |
04-Feb-1998 |
jlemon |
Document the fpgetprec/fpsetprec functions in their man page. Add cross-references to the elusive fpsetmask() function to various other man pages. Reviewed by: bde
|
32534 |
15-Jan-1998 |
danny |
PR: 5489 Submitted by: Steve G. Kargl <kargl@troutmask.apl.washington.edu> Repair corrupted text.
|
32405 |
10-Jan-1998 |
jb |
This is the only alpha math source that NetBSD has.
|
28971 |
31-Aug-1997 |
bde |
Hide the declaration of `struct exception' from C++, since it conflicts with the standard C++ `class exception'. This makes matherr() difficult to use in C++. Small loss.
|
27371 |
13-Jul-1997 |
bde |
Fixed minor bugs related to the addition of gammaf.
The major bug, that gamma is documented as really being gamma, is still unfixed.
|
25322 |
30-Apr-1997 |
bde |
Fixed wrong mnemonic `setnel' that gas happened to generate correct object code for.
Obtained from: a slightly different fix in NetBSD
|
24964 |
15-Apr-1997 |
bde |
Added -D_ARCH_INDIRECT=i387_ to CFLAGS. _ARCH_INDIRECT will soon be used to control generation of indirections in ENTRY(). Only msun needs it.
Use ${ARCH} consistently.
|
24011 |
19-Mar-1997 |
bde |
Fixed synopsis. Some float functions claimed to have the same name as the double version.
|
23579 |
09-Mar-1997 |
bde |
Use __ieee754_sqrt() instead of sqrt() internally. Similarly for the float versions. Using sqrt() was inefficient.
Obtained from: NetBSD
|
23577 |
09-Mar-1997 |
bde |
Include <machine/asm.h> instead of kernel-only <machine/asmacros.h>.
|
23397 |
05-Mar-1997 |
bde |
Fixed wrong magic numbers in scaling. hypotf() was very broken for large and small values:
hypotf(2.3819765e+38, 2.0416943e+38) was NaN instead of 3.1372484e+38 hypotf(-3.4028235e+38, 3.3886450e+38) was NaN instead of Inf hypotf(-2.8025969e-45, -2.8025969e-45) was 0 instead of 4.2038954e-45
Found by: ucbtest
|
22993 |
22-Feb-1997 |
peter |
Revert $FreeBSD$ to $Id$
|
22948 |
20-Feb-1997 |
bde |
Split up the Bessel function wrapper files so that most wrapper functions are in their own file.
|
22947 |
20-Feb-1997 |
bde |
Removed misplaced duplicate of comment about implementation details.
|
22946 |
20-Feb-1997 |
bde |
Compute (1 - x^2) as ((1 - x) * (1 + x)) instead of as (1 - x * x) to avoid easily avoidable loss of precision when |x| is nearly 1.
Extended (64-bit) precision only moves the meaning of "nearly" here.
This probably could be done better by splitting up the range into |x| <= 0.5 and |x| > 0.5 like the C version. However, ucbtest does't report any errors in this version. Perhaps the C version should be used anyway. It's only 25% slower now on a P5, provided the C version of sqrt() isn't used, and the C version could be optimized better.
Errors checked by: ucbtest
|
22808 |
16-Feb-1997 |
bde |
Select between the generic math functions and the i387-specific ones at runtime.
etc/make.conf: Nuked HAVE_FPU option.
lib/msun/Makefile: Always build the i387 objects. Copy the i387 source files at build time so that the i387 objects have different names. This is simpler than renaming the files in the cvs repository or repeating half of bsd.lib.mk to add explicit rules.
lib/msun/src/*.c: Renamed all functions that have an i387-specific version by adding `__generic_' to their names.
lib/msun/src/get_hw_float.c: New file for getting machdep.hw_float from the kernel.
sys/i386/include/asmacros.h: Abuse the ENTRY() macro to generate jump vectors and associated code. This works much like PIC PLT dynamic initialization. The PIC case is messy. The old i387 entry points are renamed. Renaming is easier here because the names are given by macro expansions.
|
22802 |
16-Feb-1997 |
bde |
Fixed the i87 version of exp(). It returned NaN for args +-Inf. It had some small (one or two ULP) inaccuracies.
Found by: ucbtest
|
22748 |
15-Feb-1997 |
jkh |
Put back .endif clobbered by the previous commit, breaking the build.
|
22731 |
15-Feb-1997 |
bde |
Disabled the i387 version if log1p(). It just evaluates log(1 + x). This defeats the point of log1p(). ucbtest reports errors of +-5e+15 ULPs. A correct version would use the i387 fyl2xp1 instruction for small x and maybe scale to small x. The C version does the scaling reasonably efficiently, and fyl2px1 is slow (at least on P5s), so not much is lost by always using the C version (only 25% for small x even with the broken i387 version; 50% for large x).
|
21907 |
20-Jan-1997 |
wosch |
Sort cross references.
|
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.
|
21435 |
08-Jan-1997 |
wollman |
Delete -D_POSIX_MODE and -D_MULTI_LIBM from CFLAGS. They never had any effect because _IEEE_LIBM always takes priority, so the definition just served to confuse.
Reviewed by: bde
|
20888 |
23-Dec-1996 |
wosch |
comma typos
|
20648 |
18-Dec-1996 |
bde |
Removed references to nonexistent functions log2() and log2f().
|
20447 |
14-Dec-1996 |
bde |
Fixed fiddling with the control word. Use the stack space reserved for that purpose instead of space below the stack.
|
20445 |
14-Dec-1996 |
bde |
Clean up the FPU stack before returning. One stack slot was leaked on every call. The damage was sometimes limited by other routines using and freeing stack slots that should have been empty to being with.
|
17958 |
30-Aug-1996 |
peter |
consistancy fixup
Submitted by: "Philippe Charnier" <charnier@xp11.frmug.org>
|
17935 |
30-Aug-1996 |
peter |
cmp -s || install -c ==> install -C
|
17860 |
28-Aug-1996 |
bde |
Made rintf() actually work. It was completely broken (when s_rint.c was compiled with -O) by the precision bug in the i386 version of gcc (assignments and casts don't clip the precision). E.g., rintf(12.3456789) was 12.125.
Avoid the same bug in rint(). It was only broken for the unusual case when the i387 precision is 64 bits. FreeBSD defaults to 53 bit precision to avoid problems like this, but the standard math emulator always uses 64 bit precision.
|
17756 |
21-Aug-1996 |
mpp |
Fix up the NAME lines forthe ceil and floor man pages to be less confusing.
Reviewed by: bde Partially obtained from: NetBSD-bugs
|
17141 |
12-Jul-1996 |
jkh |
General -Wall warning cleanup, part I. Submitted-By: Kent Vander Velden <graphix@iastate.edu>
|
16054 |
01-Jun-1996 |
bde |
Clean up the FP stack before returning. The i387 exp() leaked an FP register on its first call. Subsequent calls reused the register so the leak didn't accumulate. Fixes PR 1275.
|
14251 |
25-Feb-1996 |
bde |
Don't trash %ebp. Obtained from: NetBSD
|
14042 |
12-Feb-1996 |
mpp |
Fixed some minor formatting problems to silence manck some more. Corrected some bogus cross references to man pages that we don't/won't have and either deleted them, or found a more appropriate man page that we do have. Various other minor changes to silence manck.
Manck is currently down to about 200 lines of errors, down from the 500 - 600+ when I started all this.
|
13988 |
09-Feb-1996 |
mpp |
Correct one small typo in previous commit.
|
13987 |
09-Feb-1996 |
mpp |
Added some missing MLINKS for section 3 man pages. Also corrected a few minor formatting errors, file location and cross references in some of the section 3 man pages.
This shuts up a lot of the output from "manck" for section 3.
|
11682 |
22-Oct-1995 |
bde |
Undo the the changes in the previous revision (MANSRC now works right again). Use ${INSTALL} instead of install.
|
11674 |
22-Oct-1995 |
bde |
Fixed use of too many args for `.Em'. Fixed description of domain of y*(). Fixed description of error domain. (This description is still half redundant and half wrong, as in many other math man pages. fdlibm doesn't support the VAX or Tahoe.) Fixed capitalization of `Bessel'.
|
11136 |
02-Oct-1995 |
wollman |
Compress manual pages (if desired) in the obj directory rather than in the installation destination. Should make release-building substantially faster. The msun Makefile changes simple adapt to the new scheme.
|
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").
|
7658 |
07-Apr-1995 |
bde |
Submitted by: J.T. Conklin <jtc@wimsey.com>
First part of update to fdlibm 5.2: fix jn(n, x) and jnf(n, x). jn(-1, x) was too large by a factor of 3.
|
6953 |
08-Mar-1995 |
bde |
Obtained from: NetBSD
Remove common sources from ${SRCS} when they are replaced by arch-specific sources.
|
6794 |
01-Mar-1995 |
jkh |
Additions from Thomas Graichen to mention each functions' floating point counterpart. Submitted by: Thomas Graichen <graichen@sirius.physik.fu-berlin.de>
|
4363 |
11-Nov-1994 |
ljo |
Add missing z_abs. In BSD tradition this is in libm.a.
|
2574 |
08-Sep-1994 |
bde |
Install math.h.
|
2139 |
19-Aug-1994 |
jkh |
Latest fix from jtc: The fyl2xp1 instruction has such a limited range: -(1 - (sqrt(2) / 2)) <= x <= sqrt(2) - 1 it's not worth trying to use it.
Also, I'm not sure fyl2xp1's extra precision will matter once the result is converted from extended real (80 bits) back to double real (64 bits).
Reviewed by: jkh Submitted by: jtc
|
2122 |
19-Aug-1994 |
jkh |
Make this puppy actually compile now. Submitted by: jkh
|
2121 |
19-Aug-1994 |
jkh |
Do all the includes: <machine/asm.h> -> <machine/asmacros.h> Reviewed by: Submitted by:
|
2120 |
19-Aug-1994 |
jkh |
Change includes to reference <machine/asmacros.h>. Submitted by: jkh
|
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.
|
1573 |
27-May-1994 |
rgrimes |
BSD 4.4 Lite Lib Sources
|