#
e75a1bbc |
|
23-Apr-2024 |
Karl Tomlinson <karlt+@karlt.net> |
pow,powf(3),__ieee754_rem_pio2(f): Avoid negative integer left shift UB A compiler clever enough to know that z is positive with a non-zero biased exponent could, for example, optimize away the scalbnf(z,n) in pow() because behavior for left shift of negative values is undefined. `n` is negative when y*log2(|x|) < -0.5. i.e. |x^y| < sqrt(0.5) The intended behavior for operator<< in this code is to shift the two's complement representation of the first operand. In the pow() functions, the result is added to the IEEE 754 exponent of z = 2^y'. n may be negative enough to underflow the biased IEEE 754 exponent below zero, which is manifested in the sign bit of j (which would correspond to the IEEE 754 sign bit). The conversion from uint32_t to int32_t for out-of-int32_t-range values is implementation defined. The assumed behavior of interpreting the uint32_t value as a two's complement representation of a signed value is already assumed in many parts of the code, such as uses of GET_FLOAT_WORD() with signed integers. This code passes all the current tests, and makes some out of tree fuzzing tests pass again rather than hit UB (detailed in the commentary of the pull request). Signed-off-by: Karl Tomlinson <karlt+@karlt.net> Reviewed by: imp, steve kargl, dim Pull Request: https://github.com/freebsd/freebsd-src/pull/1137
|
#
0dd5a560 |
|
28-Jan-2024 |
Steve Kargl <kargl@FreeBSD.org> |
lib/msun: Cleanup after $FreeBSD$ removal Remove no longer needed explicit inclusion of sys/cdefs.h. PR: 276669 MFC after: 1 week
|
#
1d386b48 |
|
16-Aug-2023 |
Warner Losh <imp@FreeBSD.org> |
Remove $FreeBSD$: one-line .c pattern Remove /^[\s*]*__FBSDID\("\$FreeBSD\$"\);?\s*\n/
|
#
99843eb8 |
|
03-Aug-2023 |
Steve Kargl <kargl@FreeBSD.org> |
Clean up libm use of the __ieee754_ prefix This removes the __ieee754_ prefix from a number of the math functions. msun/src/math_private.h contains the statement that /* * ieee style elementary functions * * We rename functions here to improve other sources' diffability * against fdlibm. */ #define __ieee754_sqrt sqrt ... Here, fdlibm refers to https://netlib.org/fdlibm. It is seen from https://netlib.org/fdlibm/readme that this prefix was used to differentiate between different standards: Wrapper functions will twist the result of the ieee754 function to comply to the standard specified by the value of _LIB_VERSION if _LIB_VERSION = _IEEE_, return the ieee754 result; if _LIB_VERSION = _SVID_, return SVID result; if _LIB_VERSION = _XOPEN_, return XOPEN result; if _LIB_VERSION = _POSIX_, return POSIX/ANSI result. (These are macros, see fdlibm.h for their definition.) AFAICT, FreeBSD has never supported these wrappers. In addition, as C99, principally the long double, functions were added to libm, this convention was not maintained. Given that only 148 of 324 files under lib/msun contain a "Copyright (C) 1993 by Sun Microsystems" statement, the removal of the __ieee754_ prefix provides consistency across all source files. The last time someone compared lib/msun to fdlibm appears to be commit 3f70824172feb82ea3dcdb3866b54fe0eb7cd890 Author: David Schultz <das@FreeBSD.org> Date: Fri Feb 4 18:26:06 2005 +0000 Reduce diffs against vendor source (Sun fdlibm 5.3). The most recent fdlibm RCS string that appears in a Sun Microsystem copyrighted file is date "95/01/18". With Oracle Corporation's acquisition of Sun Microsystems in 2009, it is unlikely that fdlibm will ever be updated. A search for fdlibm at https://opensource.oracle.com/ yields no hits. Finally, OpenBSD removed the use of this prefix over 21 years ago. pSee revision 1.6 of OpenBSD's math_private.h. Note: this does not drop the __ieee754_ prefix from the trigonometric argument reduction functions, e.g., __ieee754_rem_pio2. These functions are internal to the libm and exported through Symbol.map; and thus, reserved for the implementation. PR: 272783 MFC after: 1 week
|
#
292815ea |
|
06-Sep-2021 |
Mark Murray <markm@FreeBSD.org> |
Fix powf(). Summary: From Steve Kargl: Paul Zimmermann has identified a bug in Openlibm's powf(), which is identical to FreeBSD's libm. Both derived from fdlibm. https://github.com/JuliaMath/openlibm/issues/212. Consider % cat h.c int main(void) { float x, y, z; x = 0x1.ffffecp-1F; y = -0x1.000002p+27F; z = 0x1.557a86p115F; printf("%e %e %e <-- should be %e\n", x, y, powf(x,y), z); return 0; } % cc -o h -fno-builtin h.c -lm && ./h 9.999994e-01 -1.342177e+08 inf <-- should be 5.540807e+34 Reviewers: manu Subscribers: imp, andrew, emaste Differential Revision: https://reviews.freebsd.org/D31865
|
#
93fc6789 |
|
08-Feb-2021 |
Steve Kargl <sgk@troutmask.apl.washington.edu> |
Fix incorrect powf(3) result with x near 1 and |y| much larger than 1 This adjusts the check to trigger overflow/underflow to a slightly lower value. Before: powf(9.999995e-01, -1.342177e+08) -> inf After: powf(9.999995e-01, -1.342177e+08) -> 1.858724e+31 MFC after: 1 week
|
#
6f1b8a07 |
|
17-Jul-2018 |
Bruce Evans <bde@FreeBSD.org> |
Add a macro nan_mix() and use it to get NaN results that are (bitwise) independent of the precision in most cases. This is mainly to simplify checking for errors. r176266 did this for e_pow[f].c using a less refined expression that often didn't work. r176276 fixes an error in the log message for r176266. The main refinement is to always expand to long double precision. See old log messages (especially these 2) and the comment on the macro for more general details. Specific details: - using nan_mix() consistently for the new and old pow*() functions was the only thing needed to make my consistency test for powl() vs pow() pass on amd64. - catrig[fl].c already had all the refinements, but open-coded. - e_atan2[fl].c, e_fmod[fl].c and s_remquo[fl] only had primitive NaN mixing. - e_hypot[fl].c already had a different refined version of r176266. Refine this further. nan_mix() is not directly usable here since we want to clear the sign bit. - e_remainder[f].c already had an earlier version of r176266. - s_ccosh[f].c,/s_csinh[f].c already had a version equivalent to r176266. Refine this further. nan_mix() is not directly usable here since the expression has to handle some non-NaN cases. - s_csqrt.[fl]: the mixing was special and mostly wrong. Partially fix the special version. - s_ctanh[f].c already had a version of r176266.
|
#
7525d42f |
|
08-Jul-2018 |
Mark Johnston <markj@FreeBSD.org> |
Reduce diff between msun/src/e_pow.c and msun/src/e_powf.c. Remove unnecessary casts, use integer literal constants instead of floating point constants where possible, and introduce three const static variables to hold 0.5, 0.25, and 1/3. PR: 229420 Submitted by: Steve Kargl <sgk@troutmask.apl.washington.edu> MFC after: 1 week
|
#
cd24d798 |
|
21-Oct-2011 |
David Schultz <das@FreeBSD.org> |
Per IEEE754r, pow(1, y) is 1 even if y is NaN, and pow(-1, +-Inf) is 1. MFC after: 2 weeks
|
#
a7d5f7eb |
|
19-Oct-2010 |
Jamie Gritton <jamie@FreeBSD.org> |
A new jail(8) with a configuration file, to replace the work currently done by /etc/rc.d/jail.
|
#
fe0506d7 |
|
09-Mar-2010 |
Marcel Moolenaar <marcel@FreeBSD.org> |
Create the altix project branch. The altix project will add support for the SGI Altix 350 to FreeBSD/ia64. The hardware used for porting is a two-module system, consisting of a base compute module and a CPU expansion module. SGI's NUMAFlex architecture can be an excellent platform to test CPU affinity and NUMA-aware features in FreeBSD.
|
#
d7f03759 |
|
19-Oct-2008 |
Ulf Lilleengen <lulf@FreeBSD.org> |
- Import the HEAD csup code which is the basis for the cvsmode work.
|
#
b4437c3d |
|
14-Feb-2008 |
Bruce Evans <bde@FreeBSD.org> |
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.
|
#
011cbae1 |
|
14-Feb-2008 |
Bruce Evans <bde@FreeBSD.org> |
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().
|
#
f88a48cc |
|
01-Jun-2004 |
Bruce Evans <bde@FreeBSD.org> |
Merged from double precision case (e_pow.c 1.10: sign fixes).
|
#
5f20e5ce |
|
01-Jun-2004 |
Bruce Evans <bde@FreeBSD.org> |
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.
|
#
12be4e0d |
|
01-Jun-2004 |
Bruce Evans <bde@FreeBSD.org> |
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.
|
#
3e2ec6ea |
|
17-Jun-2002 |
Bruce Evans <bde@FreeBSD.org> |
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
|
#
59b19ff1 |
|
28-May-2002 |
Alfred Perlstein <alfred@FreeBSD.org> |
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
|
#
2dcc2286 |
|
28-May-2002 |
Alfred Perlstein <alfred@FreeBSD.org> |
Assume __STDC__, remove non-__STDC__ code. Reviewed by: md5
|
#
7f3dea24 |
|
27-Aug-1999 |
Peter Wemm <peter@FreeBSD.org> |
$Id$ -> $FreeBSD$
|
#
9faa8dc6 |
|
09-Mar-1997 |
Bruce Evans <bde@FreeBSD.org> |
Use __ieee754_sqrt() instead of sqrt() internally. Similarly for the float versions. Using sqrt() was inefficient. Obtained from: NetBSD
|
#
7e546392 |
|
22-Feb-1997 |
Peter Wemm <peter@FreeBSD.org> |
Revert $FreeBSD$ to $Id$
|
#
1130b656 |
|
14-Jan-1997 |
Jordan K. Hubbard <jkh@FreeBSD.org> |
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.
|
#
6c06b4e2 |
|
29-May-1995 |
Rodney W. Grimes <rgrimes@FreeBSD.org> |
Remove trailing whitespace.
|
#
3a8617a8 |
|
19-Aug-1994 |
Jordan K. Hubbard <jkh@FreeBSD.org> |
J.T. Conklin's latest version of the Sun math library. -- Begin comments from J.T. Conklin: The most significant improvement is the addition of "float" versions of the math functions that take float arguments, return floats, and do all operations in floating point. This doesn't help (performance) much on the i386, but they are still nice to have. The float versions were orginally done by Cygnus' Ian Taylor when fdlibm was integrated into the libm we support for embedded systems. I gave Ian a copy of my libm as a starting point since I had already fixed a lot of bugs & problems in Sun's original code. After he was done, I cleaned it up a bit and integrated the changes back into my libm. -- End comments Reviewed by: jkh Submitted by: jtc
|