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.
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
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.
[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.
(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.
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
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>
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).
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
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.
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.
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.
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.
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
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>
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)
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
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
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
- 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
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".
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>
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.
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>.
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.
#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>
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
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
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.
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).
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.
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.
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.
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.
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'.
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").
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
-- 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
