diff options
| author | Jeff Johnston <jjohnstn@redhat.com> | 2009-03-25 22:13:24 +0300 |
|---|---|---|
| committer | Jeff Johnston <jjohnstn@redhat.com> | 2009-03-25 22:13:24 +0300 |
| commit | 139f923bb4564fb021dce2b7ae460a3d8e3bb87d (patch) | |
| tree | 96c1d8794c6828033a418bce8d8f1cc6eccba037 /newlib/libm | |
| parent | baea11cd6a53ccde63b5fefa4c7ae2b5426a668d (diff) | |
2009-03-25 Craig Howland <howland@LGSInnovations.com>
* libc/include/math.h: (llround, llroundf): Declare.
* libm/common/s_llround.c: New file, implementing llround().
* libm/common/sf_llround.c: New file, implementing llroundf().
* libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS
version of function.
* libm/common/sf_lrint.c: Ditto.
* libm/common/sf_logb.c: Corrected return for subnormal argument
by replacing existing function with a version created from sf_ilogb.c.
* libm/common/s_logb.c: Ditto, except starting point s_ilogb.c. Also
added documentation for logb() and logbf().
* libm/common/s_signbit.c: Add signbit() documentation.
* libm/common/s_log2.c: Update return values to match what w_log2.c has,
since log2 uses log(); add note about being derived instead of direct.
* libm/common/sf_fma.c: Add casts to attempt to get correct results,
as well as comments pointing out problems with the implementation.
* libm/common/s_fma.c: Add fma() and fmaf() documentation.
* libm/common/sf_remquo.c: Incorrect quotient returns for large values
corrected by discarding existing function and replacing with Sun
verion, with some enhancements.
* libm/common/s_remquo.c: Ditto. Add remquo() and remquof()
documentation.
* libm/common/s_fmax.c: Add fmax() and fmaxf() documentation.
* libm/common/s_fmin.c: Add fmin() and fminf() documentation.
* libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf()
documentation.
* libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg.
* libm/common/s_trunc.c: Add trunc() and truncf() documentation.
* libm/common/s_rint.c: Add rint() and rintf() documentation.
* libm/common/s_round.c: Add round() and roundf() documentation.
* libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation.
* libm/common/s_infinity.c: Add infinity() and infinityf()
documentation.
* libm/common/s_lround.c: Add lround(), lroundf(), llround(), and
llroundf() documentation.
* libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf()
documentation.
* libm/common/isgreater.c: New file for documenting math.h function-like
macros isgreater(), isgreaterequal(), isless(), islessequal(),
islessgreater(), and isunordered().
* libm/common/s_isnan.c: Add documentation for function-like macros
fpclassify(), isfinite(), isinf(), isnan(), and isnormal().
* libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf()
documentation.
* libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc);
s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def,
s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def,
s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and
isgreater.def (chobj);
re-name all existing chew files (chobj) to match source file base
names (put in underscores), delete all special targets for chew files
(leaving all to be generated by rule).
* libm/common/Makefile.in: regenerate.
* libm/math/w_exp2.c: Add "base 2" to documentation description (and
delete TRAD_SYNOPSIS).
* libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along
with some history behind the function names.
* libm/math/math.tex: Add includes for newly-added documentation (see
.def additions to common/Makefile.am and math/Makefile.am in this
ChangeLog list), adjusted existing .def file names to match source file
base names (added underscores); add mention of HUGE_VALF; rename
"Version of library" section to "Error Handling" and add some text
about floating-point exception; added section "Standards Compliance And
Portability".
* libm/math/Makefile.am: Add w_exp2.def (chobj);
re-name all existing chew files (chobj) to match source file base
names, delete all special targets for chew files (leaving all to be
generated by rule).
* libm/math/Makefile.in: regenerated
* doc/makedoc.c: Change silent ignoring of commands < 5 characters
to a failure when reading macro file for commands < 4 characters;
add -v (verbose) option for printing some debugging information;
get rid of spurious translation of "@*" to "*" (no source files used @*,
so no existing doc pages were affected); clean up some compiler
warnings.
* doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl
which has them); Remove ITEM command (redundant with makedoc built-in
"o", not used in any present source file so nothing is lost, anyway).
* HOWTO: New file to hold information for maintainers regarding how
to do things. Initial sections on documentation and ELIX levels.
Diffstat (limited to 'newlib/libm')
33 files changed, 1370 insertions, 655 deletions
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am index 1369c25bf..95387f393 100644 --- a/newlib/libm/common/Makefile.am +++ b/newlib/libm/common/Makefile.am @@ -10,7 +10,7 @@ src = s_finite.c s_copysign.c s_modf.c s_scalbn.c \ s_log1p.c s_nan.c s_nextafter.c s_pow10.c \ s_rint.c s_logb.c s_log2.c s_matherr.c s_lib_ver.c \ s_fdim.c s_fma.c s_fmax.c s_fmin.c s_fpclassify.c s_lrint.c \ - s_lround.c s_nearbyint.c s_remquo.c s_round.c s_scalbln.c \ + s_lround.c s_llround.c s_nearbyint.c s_remquo.c s_round.c s_scalbln.c \ s_signbit.c s_trunc.c fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ @@ -19,7 +19,7 @@ fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ sf_log1p.c sf_nan.c sf_nextafter.c sf_pow10.c \ sf_rint.c sf_logb.c sf_log2.c \ sf_fdim.c sf_fma.c sf_fmax.c sf_fmin.c sf_fpclassify.c sf_lrint.c \ - sf_lround.c sf_nearbyint.c sf_remquo.c sf_round.c \ + sf_lround.c sf_llround.c sf_nearbyint.c sf_remquo.c sf_round.c \ sf_scalbln.c sf_trunc.c libcommon_la_LDFLAGS = -Xcompiler -nostdlib @@ -37,9 +37,13 @@ endif # USE_LIBTOOL include $(srcdir)/../../Makefile.shared -chobj = scbrt.def scopysign.def sexp10.def sexpm1.def silogb.def \ - sinfinity.def sisnan.def slog1p.def smatherr.def smodf.def \ - snan.def snextafter.def spow10.def sscalbn.def +chobj = s_cbrt.def s_copysign.def s_exp10.def s_expm1.def s_ilogb.def \ + s_infinity.def s_isnan.def s_log1p.def s_matherr.def s_modf.def \ + s_nan.def s_nextafter.def s_pow10.def s_scalbn.def \ + s_fdim.def s_fma.def s_fmax.def s_fmin.def \ + s_logb.def s_log2.def s_lrint.def s_lround.def s_nearbyint.def \ + s_remquo.def s_rint.def s_round.def s_signbit.def s_trunc.def \ + isgreater.def SUFFIXES = .def @@ -55,73 +59,6 @@ doc: $(chobj) CLEANFILES = $(chobj) *.ref -# Texinfo does not appear to support underscores in file names, so we -# name the .def files without underscores. - -scopysign.def: s_copysign.c - $(CHEW) < $(srcdir)/s_copysign.c >$@ 2>/dev/null - touch stmp-def - -scbrt.def: s_cbrt.c - $(CHEW) < $(srcdir)/s_cbrt.c >$@ 2>/dev/null - touch stmp-def - -serf.def: s_erf.c - $(CHEW) < $(srcdir)/s_serf.c >$@ 2>/dev/null - touch stmp-def - -sexp10.def: s_exp10.c - $(CHEW) < $(srcdir)/s_exp10.c >$@ 2>/dev/null - touch stmp-def - -sexpm1.def: s_expm1.c - $(CHEW) < $(srcdir)/s_expm1.c >$@ 2>/dev/null - touch stmp-def - -silogb.def: s_ilogb.c - $(CHEW) < $(srcdir)/s_ilogb.c >$@ 2>/dev/null - touch stmp-def - -sinfinity.def: s_infinity.c - $(CHEW) < $(srcdir)/s_infinity.c >$@ 2>/dev/null - touch stmp-def - -sisnan.def: s_isnan.c - $(CHEW) < $(srcdir)/s_isnan.c >$@ 2>/dev/null - touch stmp-def - -slog1p.def: s_log1p.c - $(CHEW) < $(srcdir)/s_log1p.c >$@ 2>/dev/null - touch stmp-def - -slog2.def: s_log2.c - $(CHEW) < $(srcdir)/s_log2.c >$@ 2>/dev/null - touch stmp-def - -smodf.def: s_modf.c - $(CHEW) < $(srcdir)/s_modf.c >$@ 2>/dev/null - touch stmp-def - -smatherr.def: s_matherr.c - $(CHEW) < $(srcdir)/s_matherr.c >$@ 2>/dev/null - touch stmp-def - -snan.def: s_nan.c - $(CHEW) < $(srcdir)/s_nan.c >$@ 2>/dev/null - touch stmp-def - -snextafter.def: s_nextafter.c - $(CHEW) < $(srcdir)/s_nextafter.c >$@ 2>/dev/null - touch stmp-def - -spow10.def: s_pow10.c - $(CHEW) < $(srcdir)/s_pow10.c >$@ 2>/dev/null - touch stmp-def - -sscalbn.def: s_scalbn.c - $(CHEW) < $(srcdir)/s_scalbn.c >$@ 2>/dev/null - touch stmp-def - # A partial dependency list. $(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h fdlibm.h diff --git a/newlib/libm/common/Makefile.in b/newlib/libm/common/Makefile.in index f84166cfb..144f9a282 100644 --- a/newlib/libm/common/Makefile.in +++ b/newlib/libm/common/Makefile.in @@ -72,9 +72,10 @@ am__objects_1 = lib_a-s_finite.$(OBJEXT) lib_a-s_copysign.$(OBJEXT) \ lib_a-s_fma.$(OBJEXT) lib_a-s_fmax.$(OBJEXT) \ lib_a-s_fmin.$(OBJEXT) lib_a-s_fpclassify.$(OBJEXT) \ lib_a-s_lrint.$(OBJEXT) lib_a-s_lround.$(OBJEXT) \ - lib_a-s_nearbyint.$(OBJEXT) lib_a-s_remquo.$(OBJEXT) \ - lib_a-s_round.$(OBJEXT) lib_a-s_scalbln.$(OBJEXT) \ - lib_a-s_signbit.$(OBJEXT) lib_a-s_trunc.$(OBJEXT) + lib_a-s_llround.$(OBJEXT) lib_a-s_nearbyint.$(OBJEXT) \ + lib_a-s_remquo.$(OBJEXT) lib_a-s_round.$(OBJEXT) \ + lib_a-s_scalbln.$(OBJEXT) lib_a-s_signbit.$(OBJEXT) \ + lib_a-s_trunc.$(OBJEXT) am__objects_2 = lib_a-sf_finite.$(OBJEXT) lib_a-sf_copysign.$(OBJEXT) \ lib_a-sf_modf.$(OBJEXT) lib_a-sf_scalbn.$(OBJEXT) \ lib_a-sf_cbrt.$(OBJEXT) lib_a-sf_exp10.$(OBJEXT) \ @@ -88,9 +89,10 @@ am__objects_2 = lib_a-sf_finite.$(OBJEXT) lib_a-sf_copysign.$(OBJEXT) \ lib_a-sf_fdim.$(OBJEXT) lib_a-sf_fma.$(OBJEXT) \ lib_a-sf_fmax.$(OBJEXT) lib_a-sf_fmin.$(OBJEXT) \ lib_a-sf_fpclassify.$(OBJEXT) lib_a-sf_lrint.$(OBJEXT) \ - lib_a-sf_lround.$(OBJEXT) lib_a-sf_nearbyint.$(OBJEXT) \ - lib_a-sf_remquo.$(OBJEXT) lib_a-sf_round.$(OBJEXT) \ - lib_a-sf_scalbln.$(OBJEXT) lib_a-sf_trunc.$(OBJEXT) + lib_a-sf_lround.$(OBJEXT) lib_a-sf_llround.$(OBJEXT) \ + lib_a-sf_nearbyint.$(OBJEXT) lib_a-sf_remquo.$(OBJEXT) \ + lib_a-sf_round.$(OBJEXT) lib_a-sf_scalbln.$(OBJEXT) \ + lib_a-sf_trunc.$(OBJEXT) @USE_LIBTOOL_FALSE@am_lib_a_OBJECTS = $(am__objects_1) \ @USE_LIBTOOL_FALSE@ $(am__objects_2) lib_a_OBJECTS = $(am_lib_a_OBJECTS) @@ -102,15 +104,16 @@ am__objects_3 = s_finite.lo s_copysign.lo s_modf.lo s_scalbn.lo \ s_log1p.lo s_nan.lo s_nextafter.lo s_pow10.lo s_rint.lo \ s_logb.lo s_log2.lo s_matherr.lo s_lib_ver.lo s_fdim.lo \ s_fma.lo s_fmax.lo s_fmin.lo s_fpclassify.lo s_lrint.lo \ - s_lround.lo s_nearbyint.lo s_remquo.lo s_round.lo s_scalbln.lo \ - s_signbit.lo s_trunc.lo + s_lround.lo s_llround.lo s_nearbyint.lo s_remquo.lo s_round.lo \ + s_scalbln.lo s_signbit.lo s_trunc.lo am__objects_4 = sf_finite.lo sf_copysign.lo sf_modf.lo sf_scalbn.lo \ sf_cbrt.lo sf_exp10.lo sf_expm1.lo sf_ilogb.lo sf_infinity.lo \ sf_isinf.lo sf_isinff.lo sf_isnan.lo sf_isnanf.lo sf_log1p.lo \ sf_nan.lo sf_nextafter.lo sf_pow10.lo sf_rint.lo sf_logb.lo \ sf_log2.lo sf_fdim.lo sf_fma.lo sf_fmax.lo sf_fmin.lo \ - sf_fpclassify.lo sf_lrint.lo sf_lround.lo sf_nearbyint.lo \ - sf_remquo.lo sf_round.lo sf_scalbln.lo sf_trunc.lo + sf_fpclassify.lo sf_lrint.lo sf_lround.lo sf_llround.lo \ + sf_nearbyint.lo sf_remquo.lo sf_round.lo sf_scalbln.lo \ + sf_trunc.lo @USE_LIBTOOL_TRUE@am_libcommon_la_OBJECTS = $(am__objects_3) \ @USE_LIBTOOL_TRUE@ $(am__objects_4) libcommon_la_OBJECTS = $(am_libcommon_la_OBJECTS) @@ -280,7 +283,7 @@ src = s_finite.c s_copysign.c s_modf.c s_scalbn.c \ s_log1p.c s_nan.c s_nextafter.c s_pow10.c \ s_rint.c s_logb.c s_log2.c s_matherr.c s_lib_ver.c \ s_fdim.c s_fma.c s_fmax.c s_fmin.c s_fpclassify.c s_lrint.c \ - s_lround.c s_nearbyint.c s_remquo.c s_round.c s_scalbln.c \ + s_lround.c s_llround.c s_nearbyint.c s_remquo.c s_round.c s_scalbln.c \ s_signbit.c s_trunc.c fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ @@ -289,7 +292,7 @@ fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \ sf_log1p.c sf_nan.c sf_nextafter.c sf_pow10.c \ sf_rint.c sf_logb.c sf_log2.c \ sf_fdim.c sf_fma.c sf_fmax.c sf_fmin.c sf_fpclassify.c sf_lrint.c \ - sf_lround.c sf_nearbyint.c sf_remquo.c sf_round.c \ + sf_lround.c sf_llround.c sf_nearbyint.c sf_remquo.c sf_round.c \ sf_scalbln.c sf_trunc.c libcommon_la_LDFLAGS = -Xcompiler -nostdlib @@ -300,9 +303,13 @@ libcommon_la_LDFLAGS = -Xcompiler -nostdlib @USE_LIBTOOL_FALSE@noinst_LIBRARIES = lib.a @USE_LIBTOOL_FALSE@lib_a_SOURCES = $(src) $(fsrc) @USE_LIBTOOL_FALSE@lib_a_CFLAGS = $(AM_CFLAGS) -chobj = scbrt.def scopysign.def sexp10.def sexpm1.def silogb.def \ - sinfinity.def sisnan.def slog1p.def smatherr.def smodf.def \ - snan.def snextafter.def spow10.def sscalbn.def +chobj = s_cbrt.def s_copysign.def s_exp10.def s_expm1.def s_ilogb.def \ + s_infinity.def s_isnan.def s_log1p.def s_matherr.def s_modf.def \ + s_nan.def s_nextafter.def s_pow10.def s_scalbn.def \ + s_fdim.def s_fma.def s_fmax.def s_fmin.def \ + s_logb.def s_log2.def s_lrint.def s_lround.def s_nearbyint.def \ + s_remquo.def s_rint.def s_round.def s_signbit.def s_trunc.def \ + isgreater.def SUFFIXES = .def CHEW = ../../doc/makedoc -f $(srcdir)/../../doc/doc.str @@ -555,6 +562,12 @@ lib_a-s_lround.o: s_lround.c lib_a-s_lround.obj: s_lround.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-s_lround.obj `if test -f 's_lround.c'; then $(CYGPATH_W) 's_lround.c'; else $(CYGPATH_W) '$(srcdir)/s_lround.c'; fi` +lib_a-s_llround.o: s_llround.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-s_llround.o `test -f 's_llround.c' || echo '$(srcdir)/'`s_llround.c + +lib_a-s_llround.obj: s_llround.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-s_llround.obj `if test -f 's_llround.c'; then $(CYGPATH_W) 's_llround.c'; else $(CYGPATH_W) '$(srcdir)/s_llround.c'; fi` + lib_a-s_nearbyint.o: s_nearbyint.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-s_nearbyint.o `test -f 's_nearbyint.c' || echo '$(srcdir)/'`s_nearbyint.c @@ -753,6 +766,12 @@ lib_a-sf_lround.o: sf_lround.c lib_a-sf_lround.obj: sf_lround.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_lround.obj `if test -f 'sf_lround.c'; then $(CYGPATH_W) 'sf_lround.c'; else $(CYGPATH_W) '$(srcdir)/sf_lround.c'; fi` +lib_a-sf_llround.o: sf_llround.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_llround.o `test -f 'sf_llround.c' || echo '$(srcdir)/'`sf_llround.c + +lib_a-sf_llround.obj: sf_llround.c + $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_llround.obj `if test -f 'sf_llround.c'; then $(CYGPATH_W) 'sf_llround.c'; else $(CYGPATH_W) '$(srcdir)/sf_llround.c'; fi` + lib_a-sf_nearbyint.o: sf_nearbyint.c $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_nearbyint.o `test -f 'sf_nearbyint.c' || echo '$(srcdir)/'`sf_nearbyint.c @@ -943,73 +962,6 @@ objectlist.awk.in: $(noinst_LTLIBRARIES) doc: $(chobj) -# Texinfo does not appear to support underscores in file names, so we -# name the .def files without underscores. - -scopysign.def: s_copysign.c - $(CHEW) < $(srcdir)/s_copysign.c >$@ 2>/dev/null - touch stmp-def - -scbrt.def: s_cbrt.c - $(CHEW) < $(srcdir)/s_cbrt.c >$@ 2>/dev/null - touch stmp-def - -serf.def: s_erf.c - $(CHEW) < $(srcdir)/s_serf.c >$@ 2>/dev/null - touch stmp-def - -sexp10.def: s_exp10.c - $(CHEW) < $(srcdir)/s_exp10.c >$@ 2>/dev/null - touch stmp-def - -sexpm1.def: s_expm1.c - $(CHEW) < $(srcdir)/s_expm1.c >$@ 2>/dev/null - touch stmp-def - -silogb.def: s_ilogb.c - $(CHEW) < $(srcdir)/s_ilogb.c >$@ 2>/dev/null - touch stmp-def - -sinfinity.def: s_infinity.c - $(CHEW) < $(srcdir)/s_infinity.c >$@ 2>/dev/null - touch stmp-def - -sisnan.def: s_isnan.c - $(CHEW) < $(srcdir)/s_isnan.c >$@ 2>/dev/null - touch stmp-def - -slog1p.def: s_log1p.c - $(CHEW) < $(srcdir)/s_log1p.c >$@ 2>/dev/null - touch stmp-def - -slog2.def: s_log2.c - $(CHEW) < $(srcdir)/s_log2.c >$@ 2>/dev/null - touch stmp-def - -smodf.def: s_modf.c - $(CHEW) < $(srcdir)/s_modf.c >$@ 2>/dev/null - touch stmp-def - -smatherr.def: s_matherr.c - $(CHEW) < $(srcdir)/s_matherr.c >$@ 2>/dev/null - touch stmp-def - -snan.def: s_nan.c - $(CHEW) < $(srcdir)/s_nan.c >$@ 2>/dev/null - touch stmp-def - -snextafter.def: s_nextafter.c - $(CHEW) < $(srcdir)/s_nextafter.c >$@ 2>/dev/null - touch stmp-def - -spow10.def: s_pow10.c - $(CHEW) < $(srcdir)/s_pow10.c >$@ 2>/dev/null - touch stmp-def - -sscalbn.def: s_scalbn.c - $(CHEW) < $(srcdir)/s_scalbn.c >$@ 2>/dev/null - touch stmp-def - # A partial dependency list. $(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h fdlibm.h diff --git a/newlib/libm/common/isgreater.c b/newlib/libm/common/isgreater.c new file mode 100644 index 000000000..bd4e95e05 --- /dev/null +++ b/newlib/libm/common/isgreater.c @@ -0,0 +1,75 @@ +/* isgreater.c: This file contains no source code, but rather only the + * man-page comments. All of the documented "functions" are actually macros + * defined in math.h (q.v.). */ +/* +FUNCTION +<<isgreater>>, <<isgreaterequal>>, <<isless>>, <<islessequal>>, <<islessgreater>>, and <<isunordered>>--comparison macros +INDEX + isgreater +INDEX + isgreaterequal +INDEX + isless +INDEX + islessequal +INDEX + islessgreater +INDEX + isunordered + +ANSI_SYNOPSIS + #include <math.h> + int isgreater(real-floating <[x]>, real-floating <[y]>); + int isgreaterequal(real-floating <[x]>, real-floating <[y]>); + int isless(real-floating <[x]>, real-floating <[y]>); + int islessequal(real-floating <[x]>, real-floating <[y]>); + int islessgreater(real-floating <[x]>, real-floating <[y]>); + int isunordered(real-floating <[x]>, real-floating <[y]>); + +DESCRIPTION +<<isgreater>>, <<isgreaterequal>>, <<isless>>, <<islessequal>>, +<<islessgreater>>, and <<isunordered>> are macros defined for use in +comparing floating-point numbers without raising any floating-point +exceptions. + +The relational operators (i.e. <, >, <=, and >=) support the usual mathematical +relationships between numeric values. For any ordered pair of numeric +values exactly one of the relationships--less, greater, and equal--is +true. Relational operators may raise the "invalid" floating-point +exception when argument values are NaNs. For a NaN and a numeric value, or +for two NaNs, just the unordered relationship is true (i.e., if one or both +of the arguments a NaN, the relationship is called unordered). The specified +macros are quiet (non floating-point exception raising) versions of the +relational operators, and other comparison macros that facilitate writing +efficient code that accounts for NaNs without suffering the "invalid" +floating-point exception. In the synopses shown, "real-floating" indicates +that the argument is an expression of real floating type. + +Please note that saying that the macros do not raise floating-point +exceptions, it is referring to the function that they are performing. It +is certainly possible to give them an expression which causes an exception. +For example: +o+ +o NaN < 1.0 + causes an "invalid" exception, +o isless(NaN, 1.0) + does not, and +o isless(NaN*0., 1.0) + causes an exception due to the "NaN*0.", but not from the +resultant reduced comparison of isless(NaN, 1.0). +o- + +RETURNS +@comment Formatting note: "$@" forces a new line +No floating-point exceptions are raised for any of the macros.@* +The <<isgreater>> macro returns the value of (x) > (y).@* +The <<isgreaterequal>> macro returns the value of (x) >= (y).@* +The <<isless>> macro returns the value of (x) < (y).@* +The <<islessequal>> macro returns the value of (x) <= (y).@* +The <<islessgreater>> macro returns the value of (x) < (y) || (x) > (y).@* +The <<isunordered>> macro returns 1 if either of its arguments is NaN and 0 otherwise. + +PORTABILITY +C99, POSIX. + +*/ diff --git a/newlib/libm/common/s_fdim.c b/newlib/libm/common/s_fdim.c index 0010bf59d..8eb878209 100644 --- a/newlib/libm/common/s_fdim.c +++ b/newlib/libm/common/s_fdim.c @@ -3,6 +3,39 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<<fdim>>, <<fdimf>>--positive difference +INDEX + fdim +INDEX + fdimf + +ANSI_SYNOPSIS + #include <math.h> + double fdim(double <[x]>, double <[y]>); + float fdimf(float <[x]>, float <[y]>); + +DESCRIPTION +The <<fdim>> functions determine the positive difference between their +arguments, returning: +. <[x]> - <[y]> if <[x]> > <[y]>, or + @ifnottex +. +0 if <[x]> <= <[y]>, or + @end ifnottex + @tex +. +0 if <[x]> $\leq$ <[y]>, or + @end tex +. NAN if either argument is NAN. +A range error may occur. + +RETURNS +The <<fdim>> functions return the positive difference value. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" @@ -17,7 +50,9 @@ #endif { int c = __fpclassifyd(x); - if (c == FP_NAN || c == FP_INFINITE) + if (c == FP_NAN) return(x); + if (__fpclassifyd(y) == FP_NAN) return(y); + if (c == FP_INFINITE) return HUGE_VAL; return x > y ? x - y : 0.0; diff --git a/newlib/libm/common/s_fma.c b/newlib/libm/common/s_fma.c index 86958afcb..964a8284a 100644 --- a/newlib/libm/common/s_fma.c +++ b/newlib/libm/common/s_fma.c @@ -1,3 +1,41 @@ +/* +FUNCTION +<<fma>>, <<fmaf>>--floating multiply add +INDEX + fma +INDEX + fmaf + +ANSI_SYNOPSIS + #include <math.h> + double fma(double <[x]>, double <[y]>, double <[z]>); + float fmaf(float <[x]>, float <[y]>, float <[z]>); + +DESCRIPTION +The <<fma>> functions compute (<[x]> * <[y]>) + <[z]>, rounded as one ternary +operation: they compute the value (as if) to infinite precision and round once +to the result format, according to the rounding mode characterized by the value +of FLT_ROUNDS. That is, they are supposed to do this: see below. + +RETURNS +The <<fma>> functions return (<[x]> * <[y]>) + <[z]>, rounded as one ternary +operation. + +BUGS +This implementation does not provide the function that it should, purely +returning "(<[x]> * <[y]>) + <[z]>;" with no attempt at all to provide the +simulated infinite precision intermediates which are required. DO NOT USE THEM. + +If double has enough more precision than float, then <<fmaf>> should provide +the expected numeric results, as it does use double for the calculation. But +since this is not the case for all platforms, this manual cannot determine +if it is so for your case. + +PORTABILITY +ANSI C, POSIX. + +*/ + #include "fdlibm.h" #ifndef _DOUBLE_IS_32BITS diff --git a/newlib/libm/common/s_fmax.c b/newlib/libm/common/s_fmax.c index 53e1559e3..4003f9a34 100644 --- a/newlib/libm/common/s_fmax.c +++ b/newlib/libm/common/s_fmax.c @@ -3,6 +3,31 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<<fmax>>, <<fmaxf>>--maximum +INDEX + fmax +INDEX + fmaxf + +ANSI_SYNOPSIS + #include <math.h> + double fmax(double <[x]>, double <[y]>); + float fmaxf(float <[x]>, float <[y]>); + +DESCRIPTION +The <<fmax>> functions determine the maximum numeric value of their arguments. +NaN arguments are treated as missing data: if one argument is a NaN and the +other numeric, then the <<fmax>> functions choose the numeric value. + +RETURNS +The <<fmax>> functions return the maximum numeric value of their arguments. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/s_fmin.c b/newlib/libm/common/s_fmin.c index 65e0b1c6d..2a059209c 100644 --- a/newlib/libm/common/s_fmin.c +++ b/newlib/libm/common/s_fmin.c @@ -3,6 +3,31 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<<fmin>>, <<fminf>>--minimum +INDEX + fmin +INDEX + fminf + +ANSI_SYNOPSIS + #include <math.h> + double fmin(double <[x]>, double <[y]>); + float fminf(float <[x]>, float <[y]>); + +DESCRIPTION +The <<fmin>> functions determine the minimum numeric value of their arguments. +NaN arguments are treated as missing data: if one argument is a NaN and the +other numeric, then the <<fmin>> functions choose the numeric value. + +RETURNS +The <<fmin>> functions return the minimum numeric value of their arguments. + +PORTABILITY +ANSI C, POSIX. + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/s_infinity.c b/newlib/libm/common/s_infinity.c index 9f9cbebee..02a545edb 100644 --- a/newlib/libm/common/s_infinity.c +++ b/newlib/libm/common/s_infinity.c @@ -5,7 +5,7 @@ /* FUNCTION - <<infinity>>, <<infinityf>>---representation of infinity + <<infinity>>, <<infinityf>>--representation of infinity INDEX infinity @@ -17,17 +17,18 @@ ANSI_SYNOPSIS double infinity(void); float infinityf(void); -TRAD_SYNOPSIS - #include <math.h> - double infinity(); - float infinityf(); - - DESCRIPTION <<infinity>> and <<infinityf>> return the special number IEEE infinity in double- and single-precision arithmetic respectively. +PORTABILITY +<<infinity>> and <<infinityf>> are neither standard C nor POSIX. C and +POSIX require macros HUGE_VAL and HUGE_VALF to be defined in math.h, which +Newlib defines to be infinities corresponding to these archaic infinity() +and infinityf() functions in floating-point implementations which do have +infinities. + QUICKREF infinity - pure diff --git a/newlib/libm/common/s_isnan.c b/newlib/libm/common/s_isnan.c index b0c40369a..5ae6c9b81 100644 --- a/newlib/libm/common/s_isnan.c +++ b/newlib/libm/common/s_isnan.c @@ -13,8 +13,21 @@ /* FUNCTION - <<isnan>>, <<isnanf>>, <<isinf>>, <<isinff>>, <<finite>>, <<finitef>>---test for exceptional numbers +<<fpclassify>>, <<isfinite>>, <<isinf>>, <<isnan>>, and <<isnormal>>--floating-point classification macros; <<finite>>, <<finitef>>, <<isinf>>, <<isinff>>, <<isnan>>, <<isnanf>>--test for exceptional numbers +@c C99 (start +INDEX + fpclassify +INDEX + isfinite +INDEX + isinf +INDEX + isnan +INDEX + isnormal +@c C99 end) +@c SUSv2 (start INDEX isnan INDEX @@ -28,8 +41,18 @@ INDEX isinff INDEX finitef +@c SUSv2 end) ANSI_SYNOPSIS + [C99 standard macros:] + #include <math.h> + int fpclassify(real-floating <[x]>); + int isfinite(real-floating <[x]>); + int isinf(real-floating <[x]>); + int isnan(real-floating <[x]>); + int isnormal(real-floating <[x]>); + + [Archaic SUSv2 functions:] #include <ieeefp.h> int isnan(double <[arg]>); int isinf(double <[arg]>); @@ -38,30 +61,64 @@ ANSI_SYNOPSIS int isinff(float <[arg]>); int finitef(float <[arg]>); -TRAD_SYNOPSIS - #include <ieeefp.h> - int isnan(<[arg]>) - double <[arg]>; - int isinf(<[arg]>) - double <[arg]>; - int finite(<[arg]>); - double <[arg]>; - int isnanf(<[arg]>); - float <[arg]>; - int isinff(<[arg]>); - float <[arg]>; - int finitef(<[arg]>); - float <[arg]>; - - DESCRIPTION - These functions provide information on the floating-point +<<fpclassify>>, <<isfinite>>, <<isinf>>, <<isnan>>, and <<isnormal>> are macros +defined for use in classifying floating-point numbers. This is a help because +of special "values" like NaN and infinities. In the synopses shown, +"real-floating" indicates that the argument is an expression of real floating +type. These function-like macros are C99 and POSIX-compliant, and should be +used instead of the now-archaic SUSv2 functions. + +The <<fpclassify>> macro classifies its argument value as NaN, infinite, normal, +subnormal, zero, or into another implementation-defined category. First, an +argument represented in a format wider than its semantic type is converted to +its semantic type. Then classification is based on the type of the argument. +The <<fpclassify>> macro returns the value of the number classification macro +appropriate to the value of its argument: + +o+ +o FP_INFINITE + <[x]> is either plus or minus infinity; +o FP_NAN + <[x]> is "Not A Number" (plus or minus); +o FP_NORMAL + <[x]> is a "normal" number (i.e. is none of the other special forms); +o FP_SUBNORMAL + <[x]> is too small be stored as a regular normalized number (i.e. loss of precision is likely); or +o FP_ZERO + <[x]> is 0 (either plus or minus). +o- + +The "<<is>>" set of macros provide a useful set of shorthand ways for +classifying floating-point numbers, providing the following equivalent +relations: + +o+ +o <<isfinite>>(<[x]>) +returns non-zero if <[x]> is finite. (It is equivalent to +(<<fpclassify>>(<[x]>) != FP_INFINITE && <<fpclassify>>(<[x]>) != FP_NAN).) + +o <<isinf>>(<[x]>) +returns non-zero if <[x]> is infinite. (It is equivalent to +(<<fpclassify>>(<[x]>) == FP_INFINITE).) + +o <<isnan>>(<[x]>) +returns non-zero if <[x]> is NaN. (It is equivalent to +(<<fpclassify>>(<[x]>) == FP_NAN).) + +o <<isnormal>>(<[x]>) +returns non-zero if <[x]> is normal. (It is equivalent to +(<<fpclassify>>(<[x]>) == FP_NORMAL).) +o- + + The archaic SUSv2 functions provide information on the floating-point argument supplied. - There are five major number formats: + There are five major number formats ("exponent" referring to the + biased exponent in the binary-encoded number): o+ o zero - A number which contains all zero bits. + A number which contains all zero bits, excluding the sign bit. o subnormal A number with a zero exponent but a nonzero fraction. o normal @@ -85,7 +142,21 @@ DESCRIPTION and <<isinf>> are macros that operate on multiple types of floating-point. The SUSv2 standard declares <<isnan>> as a function taking double. Newlib has decided to declare - them both as macros in math.h and as functions in ieeefp.h. + them both as macros in math.h and as functions in ieeefp.h to + maintain backward compatibility. + +RETURNS +@comment Formatting note: "$@" forces a new line +The fpclassify macro returns the value corresponding to the appropriate FP_ macro.@* +The isfinite macro returns nonzero if <[x]> is finite, else 0.@* +The isinf macro returns nonzero if <[x]> is infinite, else 0.@* +The isnan macro returns nonzero if <[x]> is an NaN, else 0.@* +The isnormal macro returns nonzero if <[x]> has a normal value, else 0. + +PORTABILITY +math.h macros are C99, POSIX. + +ieeefp.h funtions are outdated and should be avoided. QUICKREF isnan - pure diff --git a/newlib/libm/common/s_llround.c b/newlib/libm/common/s_llround.c new file mode 100644 index 000000000..923f885ca --- /dev/null +++ b/newlib/libm/common/s_llround.c @@ -0,0 +1,68 @@ +/* lround adapted to be llround for Newlib, 2009 by Craig Howland. */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +long long int +llround(double x) +{ + __int32_t sign, exponent_less_1023; + /* Most significant word, least significant word. */ + __uint32_t msw, lsw; + long long int result; + + EXTRACT_WORDS(msw, lsw, x); + + /* Extract sign. */ + sign = ((msw & 0x80000000) ? -1 : 1); + /* Extract exponent field. */ + exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023; + msw &= 0x000fffff; + msw |= 0x00100000; + + if (exponent_less_1023 < 20) + { + if (exponent_less_1023 < 0) + { + if (exponent_less_1023 < -1) + return 0; + else + return sign; + } + else + { + msw += 0x80000 >> exponent_less_1023; + result = msw >> (20 - exponent_less_1023); + } + } + else if (exponent_less_1023 < (8 * sizeof (long long int)) - 1) + { + if (exponent_less_1023 >= 52) + result = ((long long int) msw << (exponent_less_1023 - 20)) | (lsw << (exponent_less_1023 - 52)); + else + { + unsigned int tmp = lsw + (0x80000000 >> (exponent_less_1023 - 20)); + if (tmp < lsw) + ++msw; + result = ((long long int) msw << (exponent_less_1023 - 20)) | (tmp >> (52 - exponent_less_1023)); + } + } + else + /* Result is too large to be represented by a long long int. */ + return (long long int)x; + + return sign * result; +} + +#endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_log2.c b/newlib/libm/common/s_log2.c index 4683960eb..7ad2cd88b 100644 --- a/newlib/libm/common/s_log2.c +++ b/newlib/libm/common/s_log2.c @@ -14,7 +14,7 @@ /* FUNCTION - <<log2>>, <<log2f>>---base 2 logarithm + <<log2>>, <<log2f>>--base 2 logarithm INDEX log2 INDEX @@ -25,28 +25,40 @@ ANSI_SYNOPSIS double log2(double <[x]>); float log2f(float <[x]>); -TRAD_SYNOPSIS - #include <math.h> - double log2(<[x]>); - double <[x]>; - - float log2f(<[x]>); - float <[x]>; - DESCRIPTION - <<log2>> returns the base 2 logarithm of <[x]>. +The <<log2>> functions compute the base-2 logarithm of <[x]>. A domain error +occurs if the argument is less than zero. A range error occurs if the +argument is zero. + +The Newlib implementations are not full, intrinisic calculations, but +rather are derivatives based on <<log>>. (Accuracy might be slightly off from +a direct calculation.) In addition to functions, they are also implemented as +macros defined in math.h: +. #define log2(x) (log (x) / _M_LOG2_E) +. #define log2f(x) (logf (x) / (float) _M_LOG2_E) +To use the functions instead, just undefine the macros first. - <<log2f>> is identical, save that it takes and returns <<float>> values. +You can use the (non-ANSI) function <<matherr>> to specify error +handling for these functions, indirectly through the respective <<log>> +function. RETURNS - On success, <<log2>> and <<log2f>> return the calculated value. - If the result underflows, the returned value is <<0>>. If the - result overflows, the returned value is <<HUGE_VAL>>. In - either case, <<errno>> is set to <<ERANGE>>. +The <<log2>> functions return +@ifnottex +<<log base-2(<[x]>)>> +@end ifnottex +@tex +$log_2(x)$ +@end tex +on success. +When <[x]> is zero, the +returned value is <<-HUGE_VAL>> and <<errno>> is set to <<ERANGE>>. +When <[x]> is negative, the returned value is NaN (not a number) and +<<errno>> is set to <<EDOM>>. You can control the error behavior via +<<matherr>>. PORTABILITY - <<log2>> and <<log2f>> are required by ISO/IEC 9899:1999 and the - System V Interface Definition (Issue 6). +C99, POSIX, System V Interface Definition (Issue 6). */ /* diff --git a/newlib/libm/common/s_logb.c b/newlib/libm/common/s_logb.c index dbddd1994..1be4cddb3 100644 --- a/newlib/libm/common/s_logb.c +++ b/newlib/libm/common/s_logb.c @@ -1,5 +1,5 @@ - -/* @(#)s_logb.c 5.1 93/09/24 */ +/* 2009 for Newlib: Sun's s_ilogb.c converted to be s_logb.c. */ +/* @(#)s_ilogb.c 5.1 93/09/24 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -10,33 +10,101 @@ * is preserved. * ==================================================== */ - /* - * double logb(x) - * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. - * Use ilogb instead. +FUNCTION + <<logb>>, <<logbf>>--get exponent of floating-point number +INDEX + logb +INDEX + logbf + +ANSI_SYNOPSIS + #include <math.h> + double logb(double <[x]>); + float logbf(float <[x]>); + +DESCRIPTION +The <<logb>> functions extract the exponent of <[x]>, as a signed integer value +in floating-point format. If <[x]> is subnormal it is treated as though it were +normalized; thus, for positive finite <[x]>, +@ifnottex +1 <= (<[x]> * FLT_RADIX to the power (-logb(<[x]>))) < FLT_RADIX. +@end ifnottex +@tex +$1 \leq ( x \cdot FLT\_RADIX ^ {-logb(x)} ) < FLT\_RADIX$. +@end tex +A domain error may occur if the argument is zero. +In this floating-point implementation, FLT_RADIX is 2. Which also means +that for finite <[x]>, <<logb>>(<[x]>) = <<floor>>(<<log2>>(<<fabs>>(<[x]>))). + +All nonzero, normal numbers can be described as +@ifnottex +<[m]> * 2**<[p]>, where 1.0 <= <[m]> < 2.0. +@end ifnottex +@tex +$m \cdot 2^p$, where $1.0 \leq m < 2.0$. +@end tex +The <<logb>> functions examine the argument <[x]>, and return <[p]>. +The <<frexp>> functions are similar to the <<logb>> functions, but +returning <[m]> adjusted to the interval [.5, 1) or 0, and <[p]>+1. + +RETURNS +@comment Formatting note: "$@" forces a new line +When <[x]> is:@* ++inf or -inf, +inf is returned;@* +NaN, NaN is returned;@* +0, -inf is returned, and the divide-by-zero exception is raised;@* +otherwise, the <<logb>> functions return the signed exponent of <[x]>. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +frexp, ilogb +*/ + +/* double logb(double x) + * return the binary exponent of non-zero x + * logb(0) = -inf, raise divide-by-zero floating point exception + * logb(+inf|-inf) = +inf (no signal is raised) + * logb(NaN) = NaN (no signal is raised) + * Per C99 recommendation, a NaN argument is returned unchanged. */ #include "fdlibm.h" #ifndef _DOUBLE_IS_32BITS +double #ifdef __STDC__ - double logb(double x) +logb(double x) #else - double logb(x) - double x; +logb(x) +double x; #endif { - __int32_t lx,ix; - EXTRACT_WORDS(ix,lx,x); - ix &= 0x7fffffff; /* high |x| */ - if((ix|lx)==0) return -1.0/fabs(x); - if(ix>=0x7ff00000) return x*x; - if((ix>>=20)==0) /* IEEE 754 logb */ - return -1022.0; - else - return (double) (ix-1023); + __int32_t hx,lx,ix; + + EXTRACT_WORDS(hx,lx,x); + hx &= 0x7fffffff; /* high |x| */ + if(hx<0x00100000) { /* 0 or subnormal */ + if((hx|lx)==0) { + double xx; + /* arg==0: return -inf and raise divide-by-zero exception */ + INSERT_WORDS(xx,hx,lx); /* +0.0 */ + return -1./xx; /* logb(0) = -inf */ + } + else /* subnormal x */ + if(hx==0) { + for (ix = -1043; lx>0; lx<<=1) ix -=1; + } else { + for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1; + } + return (double) ix; + } + else if (hx<0x7ff00000) return (hx>>20)-1023; /* normal # */ + else if (hx>0x7ff00000 || lx) return x; /* x==NaN */ + else return HUGE_VAL; /* x==inf (+ or -) */ } #endif /* _DOUBLE_IS_32BITS */ diff --git a/newlib/libm/common/s_lrint.c b/newlib/libm/common/s_lrint.c index 7d3b561b1..5c15dd5e7 100644 --- a/newlib/libm/common/s_lrint.c +++ b/newlib/libm/common/s_lrint.c @@ -10,6 +10,44 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<lrint>>, <<lrintf>>, <<llrint>>, <<llrintf>>--round to integer +INDEX + lrint +INDEX + lrintf +INDEX + llrint +INDEX + llrintf + +ANSI_SYNOPSIS + #include <math.h> + long int lrint(double <[x]>); + long int lrintf(float <[x]>); + long long int llrint(double <[x]>); + long long int llrintf(float <[x]>); + +DESCRIPTION +The <<lrint>> and <<llrint>> functions round their argument to the nearest +integer value, using the current rounding direction. If the rounded value is +outside the range of the return type, the numeric result is unspecified. A +range error may occur if the magnitude of <[x]> is too large. +The "inexact" floating-point exception is raised in implementations that +support it when the result differs in value from the argument (i.e., when +a fraction actually has been truncated). + +RETURNS +<[x]> rounded to an integral value, using the current rounding direction. + +SEEALSO +<<lround>> + +PORTABILITY +ANSI C, POSIX + +*/ /* * lrint(x) diff --git a/newlib/libm/common/s_lround.c b/newlib/libm/common/s_lround.c index 7bd249f9c..b892a743c 100644 --- a/newlib/libm/common/s_lround.c +++ b/newlib/libm/common/s_lround.c @@ -8,6 +8,44 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<lround>>, <<lroundf>>, <<llround>>, <<llroundf>>--round to integer, to nearest +INDEX + lround +INDEX + lroundf +INDEX + llround +INDEX + llroundf + +ANSI_SYNOPSIS + #include <math.h> + long int lround(double <[x]>); + long int lroundf(float <[x]>); + long long int llround(double <[x]>); + long long int llroundf(float <[x]>); + +DESCRIPTION + The <<lround>> and <<llround>> functions round their argument to the + nearest integer value, rounding halfway cases away from zero, regardless + of the current rounding direction. If the rounded value is outside the + range of the return type, the numeric result is unspecified (depending + upon the floating-point implementation, not the library). A range + error may occur if the magnitude of x is too large. + +RETURNS +<[x]> rounded to an integral value as an integer. + +SEEALSO +See the <<round>> functions for the return being the same floating-point type +as the argument. <<lrint>>, <<llrint>>. + +PORTABILITY +ANSI C, POSIX + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/s_nearbyint.c b/newlib/libm/common/s_nearbyint.c index 25e629b65..e7a022072 100644 --- a/newlib/libm/common/s_nearbyint.c +++ b/newlib/libm/common/s_nearbyint.c @@ -8,6 +8,43 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<nearbyint>>, <<nearbyintf>>--round to integer +INDEX + nearbyint +INDEX + nearbyintf + +ANSI_SYNOPSIS + #include <math.h> + double nearbyint(double <[x]>); + float nearbyintf(float <[x]>); + +DESCRIPTION +The <<nearbyint>> functions round their argument to an integer value in +floating-point format, using the current rounding direction and +(supposedly) without raising the "inexact" floating-point exception. +See the <<rint>> functions for the same function with the "inexact" +floating-point exception being raised when appropriate. + +BUGS +Newlib does not support the floating-point exception model, so that +the floating-point exception control is not present and thereby what may +be seen will be compiler and hardware dependent in this regard. +The Newlib <<nearbyint>> functions are identical to the <<rint>> +functions with respect to the floating-point exception behavior, and +will cause the "inexact" exception to be raised for most targets. + +RETURNS +<[x]> rounded to an integral value, using the current rounding direction. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +<<rint>>, <<round>> +*/ #include <math.h> #include "fdlibm.h" diff --git a/newlib/libm/common/s_remquo.c b/newlib/libm/common/s_remquo.c index bd52f39fc..9e18ab132 100644 --- a/newlib/libm/common/s_remquo.c +++ b/newlib/libm/common/s_remquo.c @@ -1,39 +1,208 @@ -/* Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. +/* Adapted for Newlib, 2009. (Allow for int < 32 bits; return *quo=0 during + * errors to make test scripts easier.) */ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * - * Permission to use, copy, modify, and distribute this software - * is freely granted, provided that this notice is preserved. + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== */ +/* +FUNCTION +<<remquo>>, <<remquof>>--remainder and part of quotient +INDEX + remquo +INDEX + remquof -#include "fdlibm.h" +ANSI_SYNOPSIS + #include <math.h> + double remquo(double <[x]>, double <[y]>, int *<[quo]>); + float remquof(float <[x]>, float <[y]>, int *<[quo]>); + +DESCRIPTION +The <<remquo>> functions compute the same remainder as the <<remainder>> +functions; this value is in the range -<[y]>/2 ... +<[y]>/2. In the object +pointed to by <<quo>> they store a value whose sign is the sign of <<x>>/<<y>> +and whose magnitude is congruent modulo 2**n to the magnitude of the integral +quotient of <<x>>/<<y>>. (That is, <<quo>> is given the n lsbs of the +quotient, not counting the sign.) This implementation uses n=31 if int is 32 +bits or more, otherwise, n is 1 less than the width of int. + +For example: +. remquo(-29.0, 3.0, &<[quo]>) +returns -1.0 and sets <[quo]>=10, and +. remquo(-98307.0, 3.0, &<[quo]>) +returns -0.0 and sets <[quo]>=-32769, although for 16-bit int, <[quo]>=-1. In +the latter case, the actual quotient of -(32769=0x8001) is reduced to -1 +because of the 15-bit limitation for the quotient. + +RETURNS +When either argument is NaN, NaN is returned. If <[y]> is 0 or <[x]> is +infinite (and neither is NaN), a domain error occurs (i.e. the "invalid" +floating point exception is raised or errno is set to EDOM), and NaN is +returned. +Otherwise, the <<remquo>> functions return <[x]> REM <[y]>. -#ifndef _DOUBLE_IS_32BITS +BUGS +IEEE754-2008 calls for <<remquo>>(subnormal, inf) to cause the "underflow" +floating-point exception. This implementation does not. -#ifdef __STDC__ - double remquo(double x, double y, int *quo) /* wrapper remquo */ -#else - double remquo(x,y,quo) /* wrapper remquo */ - double x,y; - int *quo; +PORTABILITY +C99, POSIX. + +*/ + +#include <limits.h> +#include <math.h> +#include "fdlibm.h" + +/* For quotient, return either all 31 bits that can from calculation (using + * int32_t), or as many as can fit into an int that is smaller than 32 bits. */ +#if INT_MAX > 0x7FFFFFFFL + #define QUO_MASK 0x7FFFFFFF +# else + #define QUO_MASK INT_MAX #endif + +static const double Zero[] = {0.0, -0.0,}; + +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31--if that many fit-- + * because we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +double +remquo(double x, double y, int *quo) { - int signx, signy, signres; - int mswx; - int mswy; - double x_over_y; + __int32_t n,hx,hy,hz,ix,iy,sx,i; + __uint32_t lx,ly,lz,q,sxy; - GET_HIGH_WORD(mswx, x); - GET_HIGH_WORD(mswy, y); + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hy,ly,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ - signx = (mswx & 0x80000000) >> 31; - signy = (mswy & 0x80000000) >> 31; + /* purge off exception values */ + if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ + ((hy|((ly|-ly)>>31))>0x7ff00000)) { /* or y is NaN */ + *quo = 0; /* Not necessary, but return consistent value */ + return (x*y)/(x*y); + } + if(hx<=hy) { + if((hx<hy)||(lx<ly)) { + q = 0; + goto fixup; /* |x|<|y| return x or x-y */ + } + if(lx==ly) { + *quo = 1; + return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/ + } + } - signres = (signx ^ signy) ? -1 : 1; + /* determine ix = ilogb(x) */ + if(hx<0x00100000) { /* subnormal x */ + if(hx==0) { + for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; + } else { + for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; + } + } else ix = (hx>>20)-1023; - x_over_y = fabs(x / y); + /* determine iy = ilogb(y) */ + if(hy<0x00100000) { /* subnormal y */ + if(hy==0) { + for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; + } else { + for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; + } + } else iy = (hy>>20)-1023; - *quo = signres * (lrint(x_over_y) & 0x7f); + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -1022) + hx = 0x00100000|(0x000fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -1022-ix; + if(n<=31) { + hx = (hx<<n)|(lx>>(32-n)); + lx <<= n; + } else { + hx = lx<<(n-32); + lx = 0; + } + } + if(iy >= -1022) + hy = 0x00100000|(0x000fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -1022-iy; + if(n<=31) { + hy = (hy<<n)|(ly>>(32-n)); + ly <<= n; + } else { + hy = ly<<(n-32); + ly = 0; + } + } - return remainder(x,y); -} + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} + else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} + q <<= 1; + } + hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; + if(hz>=0) {hx=hz;lx=lz;q++;} -#endif /* defined(_DOUBLE_IS_32BITS) */ + /* convert back to floating value and restore the sign */ + if((hx|lx)==0) { /* return sign(x)*0 */ + q &= QUO_MASK; + *quo = (sxy ? -q : q); + return Zero[(__uint32_t)sx>>31]; + } + while(hx<0x00100000) { /* normalize x */ + hx = hx+hx+(lx>>31); lx = lx+lx; + iy -= 1; + } + if(iy>= -1022) { /* normalize output */ + hx = ((hx-0x00100000)|((iy+1023)<<20)); + } else { /* subnormal output */ + n = -1022 - iy; + if(n<=20) { + lx = (lx>>n)|((__uint32_t)hx<<(32-n)); + hx >>= n; + } else if (n<=31) { + lx = (hx<<(32-n))|(lx>>n); hx = sx; + } else { + lx = hx>>(n-32); hx = sx; + } + } +fixup: + INSERT_WORDS(x,hx,lx); + y = fabs(y); + if (y < 0x1p-1021) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5*y || (x==0.5*y && (q & 1))) { + q++; + x-=y; + } + GET_HIGH_WORD(hx,x); + SET_HIGH_WORD(x,hx^sx); + q &= QUO_MASK; + *quo = (sxy ? -q : q); + return x; +} diff --git a/newlib/libm/common/s_rint.c b/newlib/libm/common/s_rint.c index ba33e4786..4fa5ebc6c 100644 --- a/newlib/libm/common/s_rint.c +++ b/newlib/libm/common/s_rint.c @@ -10,6 +10,40 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<rint>>, <<rintf>>--round to integer +INDEX + rint +INDEX + rintf + +ANSI_SYNOPSIS + #include <math.h> + double rint(double <[x]>); + float rintf(float <[x]>); + +DESCRIPTION + The <<rint>> functions round their argument to an integer value in + floating-point format, using the current rounding direction. They + raise the "inexact" floating-point exception if the result differs + in value from the argument. See the <<nearbyint>> functions for the + same function with the "inexact" floating-point exception never being + raised. Newlib does not directly support floating-point exceptions. + The <<rint>> functions are written so that the "inexact" exception is + raised in hardware implementations that support it, even though Newlib + does not provide access. + +RETURNS +<[x]> rounded to an integral value, using the current rounding direction. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +<<nearbyint>>, <<round>> + +*/ /* * rint(x) diff --git a/newlib/libm/common/s_round.c b/newlib/libm/common/s_round.c index 09f19de64..b2cc0f558 100644 --- a/newlib/libm/common/s_round.c +++ b/newlib/libm/common/s_round.c @@ -8,6 +8,38 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<round>>, <<roundf>>--round to integer, to nearest +INDEX + round +INDEX + roundf + +ANSI_SYNOPSIS + #include <math.h> + double round(double <[x]>); + float roundf(float <[x]>); + +DESCRIPTION + The <<round>> functions round their argument to the nearest integer + value in floating-point format, rounding halfway cases away from zero, + regardless of the current rounding direction. (While the "inexact" + floating-point exception behavior is unspecified by the C standard, the + <<round>> functions are written so that "inexact" is not raised if the + result does not equal the argument, which behavior is as recommended by + IEEE 754 for its related functions.) + +RETURNS +<[x]> rounded to an integral value. + +PORTABILITY +ANSI C, POSIX + +SEEALSO +<<nearbyint>>, <<rint>> + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/s_scalbn.c b/newlib/libm/common/s_scalbn.c index 5edf6b766..50f20efec 100644 --- a/newlib/libm/common/s_scalbn.c +++ b/newlib/libm/common/s_scalbn.c @@ -13,38 +13,44 @@ /* FUNCTION -<<scalbn>>, <<scalbnf>>---scale by power of two +<<scalbn>>, <<scalbnf>>, <<scalbln>>, <<scalblnf>>--scale by power of FLT_RADIX (=2) INDEX scalbn INDEX scalbnf +INDEX + scalbln +INDEX + scalblnf ANSI_SYNOPSIS #include <math.h> - double scalbn(double <[x]>, int <[y]>); - float scalbnf(float <[x]>, int <[y]>); - -TRAD_SYNOPSIS - #include <math.h> - double scalbn(<[x]>,<[y]>) - double <[x]>; - int <[y]>; - float scalbnf(<[x]>,<[y]>) - float <[x]>; - int <[y]>; + double scalbn(double <[x]>, int <[n]>); + float scalbnf(float <[x]>, int <[n]>); + double scalbln(double <[x]>, long int <[n]>); + float scalblnf(float <[x]>, long int <[n]>); DESCRIPTION -<<scalbn>> and <<scalbnf>> scale <[x]> by <[n]>, returning <[x]> times -2 to the power <[n]>. The result is computed by manipulating the -exponent, rather than by actually performing an exponentiation or -multiplication. +The <<scalbn>> and <<scalbln>> functions compute + @ifnottex + <[x]> times FLT_RADIX to the power <[n]>. + @end ifnottex + @tex + $x \cdot FLT\_RADIX^n$. + @end tex +efficiently. The result is computed by manipulating the exponent, rather than +by actually performing an exponentiation or multiplication. In this +floating-point implementation FLT_RADIX=2, which makes the <<scalbn>> +functions equivalent to the <<ldexp>> functions. RETURNS -<[x]> times 2 to the power <[n]>. +<[x]> times 2 to the power <[n]>. A range error may occur. PORTABILITY -Neither <<scalbn>> nor <<scalbnf>> is required by ANSI C or by the System V -Interface Definition (Issue 2). +ANSI C, POSIX + +SEEALSO +<<ldexp>> */ diff --git a/newlib/libm/common/s_signbit.c b/newlib/libm/common/s_signbit.c index d7d49cbb1..746ab4662 100644 --- a/newlib/libm/common/s_signbit.c +++ b/newlib/libm/common/s_signbit.c @@ -3,6 +3,35 @@ * Permission to use, copy, modify, and distribute this software * is freely granted, provided that this notice is preserved. */ +/* +FUNCTION +<<signbit>>--Does floating-point number have negative sign? + +INDEX + signbit + +ANSI_SYNOPSIS + #include <math.h> + int signbit(real-floating <[x]>); + +DESCRIPTION +The <<signbit>> macro determines whether the sign of its argument value is +negative. The macro reports the sign of all values, including infinities, +zeros, and NaNs. If zero is unsigned, it is treated as positive. As shown in +the synopsis, the argument is "real-floating," meaning that any of the real +floating-point types (float, double, etc.) may be given to it. + +Note that because of the possibilities of signed 0 and NaNs, the expression +"<[x]> < 0.0" does not give the same result as <<signbit>> in all cases. + +RETURNS +The <<signbit>> macro returns a nonzero value if and only if the sign of its +argument value is negative. + +PORTABILITY +C99, POSIX. + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/s_trunc.c b/newlib/libm/common/s_trunc.c index c3ca154d9..daeb7910a 100644 --- a/newlib/libm/common/s_trunc.c +++ b/newlib/libm/common/s_trunc.c @@ -8,6 +8,35 @@ * is preserved. * ==================================================== */ +/* +FUNCTION +<<trunc>>, <<truncf>>--round to integer, towards zero +INDEX + trunc +INDEX + truncf + +ANSI_SYNOPSIS + #include <math.h> + double trunc(double <[x]>); + float truncf(float <[x]>); + +DESCRIPTION + The <<trunc>> functions round their argument to the integer value, in + floating format, nearest to but no larger in magnitude than the + argument, regardless of the current rounding direction. (While the + "inexact" floating-point exception behavior is unspecified by the C + standard, the <<trunc>> functions are written so that "inexact" is not + raised if the result does not equal the argument, which behavior is as + recommended by IEEE 754 for its related functions.) + +RETURNS +<[x]> truncated to an integral value. + +PORTABILITY +ANSI C, POSIX + +*/ #include "fdlibm.h" diff --git a/newlib/libm/common/sf_fdim.c b/newlib/libm/common/sf_fdim.c index 80221ace1..fe349098b 100644 --- a/newlib/libm/common/sf_fdim.c +++ b/newlib/libm/common/sf_fdim.c @@ -15,8 +15,10 @@ #endif { int c = __fpclassifyf(x); - if (c == FP_NAN || c == FP_INFINITE) - return HUGE_VAL; + if (c == FP_NAN) return(x); + if (__fpclassifyf(y) == FP_NAN) return(y); + if (c == FP_INFINITE) + return HUGE_VALF; return x > y ? x - y : 0.0; } diff --git a/newlib/libm/common/sf_fma.c b/newlib/libm/common/sf_fma.c index 3b4bcc1ce..4360f400b 100644 --- a/newlib/libm/common/sf_fma.c +++ b/newlib/libm/common/sf_fma.c @@ -15,8 +15,14 @@ float z; #endif { - /* Let the implementation handle this. */ - return (x * y) + z; + /* NOTE: The floating-point exception behavior of this is not as + * required. But since the basic function is not really done properly, + * it is not worth bothering to get the exceptions right, either. */ + /* Let the implementation handle this. */ /* <= NONSENSE! */ + /* In floating-point implementations in which double is larger than float, + * computing as double should provide the desired function. Otherwise, + * the behavior will not be as specified in the standards. */ + return (float) (((double) x * (double) y) + (double) z); } #ifdef _DOUBLE_IS_32BITS diff --git a/newlib/libm/common/sf_llround.c b/newlib/libm/common/sf_llround.c new file mode 100644 index 000000000..fe2b4bd23 --- /dev/null +++ b/newlib/libm/common/sf_llround.c @@ -0,0 +1,55 @@ +/* lroundf adapted to be llroundf for Newlib, 2009 by Craig Howland. */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +long long int +llroundf(float x) +{ + __int32_t exponent_less_127; + __uint32_t w; + long long int result; + __int32_t sign; + + GET_FLOAT_WORD (w, x); + exponent_less_127 = ((w & 0x7f800000) >> 23) - 127; + sign = (w & 0x80000000) != 0 ? -1 : 1; + w &= 0x7fffff; + w |= 0x800000; + + if (exponent_less_127 < (int)((8 * sizeof (long long int)) - 1)) + { + if (exponent_less_127 < 0) + return exponent_less_127 < -1 ? 0 : sign; + else if (exponent_less_127 >= 23) + result = (long long int) w << (exponent_less_127 - 23); + else + { + w += 0x400000 >> exponent_less_127; + result = w >> (23 - exponent_less_127); + } + } + else + return (long long int) x; + + return sign * result; +} + +#ifdef _DOUBLE_IS_32BITS + +long long int +llround(double x) +{ + return llroundf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_logb.c b/newlib/libm/common/sf_logb.c index f193f91f6..75336a1e0 100644 --- a/newlib/libm/common/sf_logb.c +++ b/newlib/libm/common/sf_logb.c @@ -1,7 +1,4 @@ -/* sf_logb.c -- float version of s_logb.c. - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ - +/* 2009 for Newlib: Sun's sf_ilogb.c converted to be sf_logb.c. */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. @@ -13,24 +10,41 @@ * ==================================================== */ +/* float logb(float x) + * return the binary exponent of non-zero x + * logbf(0) = -inf, raise divide-by-zero floating point exception + * logbf(+inf|-inf) = +inf (no signal is raised) + * logbf(NaN) = NaN (no signal is raised) + * Per C99 recommendation, a NaN argument is returned unchanged. + */ + #include "fdlibm.h" +float #ifdef __STDC__ - float logbf(float x) +logbf(float x) #else - float logbf(x) - float x; +logbf(x) +float x; #endif { - __int32_t ix; - GET_FLOAT_WORD(ix,x); - ix &= 0x7fffffff; /* high |x| */ - if(FLT_UWORD_IS_ZERO(ix)) return (float)-1.0/fabsf(x); - if(!FLT_UWORD_IS_FINITE(ix)) return x*x; - if((ix>>=23)==0) /* IEEE 754 logb */ - return -126.0; - else - return (float) (ix-127); + __int32_t hx,ix; + + GET_FLOAT_WORD(hx,x); + hx &= 0x7fffffff; + if(FLT_UWORD_IS_ZERO(hx)) { + float xx; + /* arg==0: return -inf and raise divide-by-zero exception */ + SET_FLOAT_WORD(xx,hx); /* +0.0 */ + return -1./xx; /* logbf(0) = -inf */ + } + if(FLT_UWORD_IS_SUBNORMAL(hx)) { + for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1; + return (float) ix; + } + else if (FLT_UWORD_IS_INFINITE(hx)) return HUGE_VALF; /* x==+|-inf */ + else if (FLT_UWORD_IS_NAN(hx)) return x; + else return (float) ((hx>>23)-127); } #ifdef _DOUBLE_IS_32BITS diff --git a/newlib/libm/common/sf_lrint.c b/newlib/libm/common/sf_lrint.c index 3165705d1..045a62f03 100644 --- a/newlib/libm/common/sf_lrint.c +++ b/newlib/libm/common/sf_lrint.c @@ -95,7 +95,7 @@ TWO23[2]={ double x; #endif { - return (double) lrintf((float) x); + return lrintf((float) x); } #endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_lround.c b/newlib/libm/common/sf_lround.c index eee420b93..e1f2fa100 100644 --- a/newlib/libm/common/sf_lround.c +++ b/newlib/libm/common/sf_lround.c @@ -56,7 +56,7 @@ double x; #endif { - return (double) lroundf((float) x); + return lroundf((float) x); } #endif /* defined(_DOUBLE_IS_32BITS) */ diff --git a/newlib/libm/common/sf_remquo.c b/newlib/libm/common/sf_remquo.c index 41c2038d4..7a756af8d 100644 --- a/newlib/libm/common/sf_remquo.c +++ b/newlib/libm/common/sf_remquo.c @@ -1,50 +1,130 @@ -/* Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. +/* Adapted for Newlib, 2009. (Allow for int < 32 bits; return *quo=0 during + * errors to make test scripts easier.) */ +/* @(#)e_fmod.c 1.3 95/01/18 */ +/*- + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * - * Permission to use, copy, modify, and distribute this software - * is freely granted, provided that this notice is preserved. + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== */ +#include <math.h> #include "fdlibm.h" -#ifdef __STDC__ - float remquof(float x, float y, int *quo) /* wrapper remquof */ -#else - float remquof(x,y,quo) /* wrapper remquof */ - float x,y; - int *quo; +/* For quotient, return either all 31 bits that can from calculation (using + * int32_t), or as many as can fit into an int that is smaller than 32 bits. */ +#if INT_MAX > 0x7FFFFFFFL + #define QUO_MASK 0x7FFFFFFF +# else + #define QUO_MASK INT_MAX #endif -{ - int signx, signy, signres; - int wx; - int wy; - float x_over_y; - GET_FLOAT_WORD(wx, x); - GET_FLOAT_WORD(wy, y); +static const float Zero[] = {0.0, -0.0,}; - signx = (wx & 0x80000000) >> 31; - signy = (wy & 0x80000000) >> 31; +/* + * Return the IEEE remainder and set *quo to the last n bits of the + * quotient, rounded to the nearest integer. We choose n=31--if that many fit-- + * we wind up computing all the integer bits of the quotient anyway as + * a side-effect of computing the remainder by the shift and subtract + * method. In practice, this is far more bits than are needed to use + * remquo in reduction algorithms. + */ +float +remquof(float x, float y, int *quo) +{ + __int32_t n,hx,hy,hz,ix,iy,sx,i; + __uint32_t q,sxy; - signres = (signx ^ signy) ? -1 : 1; + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hy,y); + sxy = (hx ^ hy) & 0x80000000; + sx = hx&0x80000000; /* sign of x */ + hx ^=sx; /* |x| */ + hy &= 0x7fffffff; /* |y| */ - x_over_y = fabsf(x / y); + /* purge off exception values */ + if(hy==0||hx>=0x7f800000||hy>0x7f800000) { /* y=0,NaN;or x not finite */ + *quo = 0; /* Not necessary, but return consistent value */ + return (x*y)/(x*y); + } + if(hx<hy) { + q = 0; + goto fixup; /* |x|<|y| return x or x-y */ + } else if(hx==hy) { + *quo = 1; + return Zero[(__uint32_t)sx>>31]; /* |x|=|y| return x*0*/ + } - *quo = signres * (lrintf(x_over_y) & 0x7f); + /* determine ix = ilogb(x) */ + if(hx<0x00800000) { /* subnormal x */ + for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; + } else ix = (hx>>23)-127; - return remainderf(x,y); -} + /* determine iy = ilogb(y) */ + if(hy<0x00800000) { /* subnormal y */ + for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1; + } else iy = (hy>>23)-127; -#ifdef _DOUBLE_IS_32BITS + /* set up {hx,lx}, {hy,ly} and align y to x */ + if(ix >= -126) + hx = 0x00800000|(0x007fffff&hx); + else { /* subnormal x, shift x to normal */ + n = -126-ix; + hx <<= n; + } + if(iy >= -126) + hy = 0x00800000|(0x007fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -126-iy; + hy <<= n; + } -#ifdef __STDC__ - double remquo(double x, double y, int *quo) /* wrapper remquof */ -#else - double remquo(x,y,quo) /* wrapper remquof */ - double x,y; - int *quo; -#endif -{ - return (double) remquof((float) x, (float) y, quo); -} + /* fix point fmod */ + n = ix - iy; + q = 0; + while(n--) { + hz=hx-hy; + if(hz<0) hx = hx << 1; + else {hx = hz << 1; q++;} + q <<= 1; + } + hz=hx-hy; + if(hz>=0) {hx=hz;q++;} -#endif /* defined(_DOUBLE_IS_32BITS) */ + /* convert back to floating value and restore the sign */ + if(hx==0) { /* return sign(x)*0 */ + *quo = (sxy ? -q : q); + return Zero[(__uint32_t)sx>>31]; + } + while(hx<0x00800000) { /* normalize x */ + hx <<= 1; + iy -= 1; + } + if(iy>= -126) { /* normalize output */ + hx = ((hx-0x00800000)|((iy+127)<<23)); + } else { /* subnormal output */ + n = -126 - iy; + hx >>= n; + } +fixup: + SET_FLOAT_WORD(x,hx); + y = fabsf(y); + if (y < 0x1p-125f) { + if (x+x>y || (x+x==y && (q & 1))) { + q++; + x-=y; + } + } else if (x>0.5f*y || (x==0.5f*y && (q & 1))) { + q++; + x-=y; + } + GET_FLOAT_WORD(hx,x); + SET_FLOAT_WORD(x,hx^sx); + q &= 0x7fffffff; + *quo = (sxy ? -q : q); + return x; +} diff --git a/newlib/libm/math/Makefile.am b/newlib/libm/math/Makefile.am index 668254145..335817920 100644 --- a/newlib/libm/math/Makefile.am +++ b/newlib/libm/math/Makefile.am @@ -65,14 +65,14 @@ endif # USE_LIBTOOL include $(srcdir)/../../Makefile.shared -chobj = wacos.def wacosh.def wasin.def sasinh.def \ - satan.def watan2.def watanh.def wj0.def \ - wcosh.def serf.def wexp.def \ - sfabs.def sfloor.def wfmod.def sfrexp.def \ - wgamma.def whypot.def sldexp.def wlog.def \ - wlog10.def \ - wpow.def wremainder.def ssin.def wsinh.def \ - wsqrt.def stan.def stanh.def +chobj = w_acos.def w_acosh.def w_asin.def s_asinh.def \ + s_atan.def w_atan2.def w_atanh.def w_j0.def \ + w_cosh.def s_erf.def w_exp.def w_exp2.def \ + s_fabs.def s_floor.def w_fmod.def s_frexp.def \ + w_gamma.def w_hypot.def s_ldexp.def w_log.def \ + w_log10.def \ + w_pow.def w_remainder.def s_sin.def w_sinh.def \ + w_sqrt.def s_tan.def s_tanh.def SUFFIXES = .def @@ -89,94 +89,6 @@ doc: $(chobj) CLEANFILES = $(chobj) *.ref -# Texinfo does not appear to support underscores in file names, so we -# name the .def files without underscores. - -wacos.def: w_acos.c - $(CHEW) < $(srcdir)/w_acos.c >$@ 2>/dev/null - touch stmp-def -wacosh.def: w_acosh.c - $(CHEW) < $(srcdir)/w_acosh.c >$@ 2>/dev/null - touch stmp-def -wasin.def: w_asin.c - $(CHEW) < $(srcdir)/w_asin.c >$@ 2>/dev/null - touch stmp-def -sasinh.def: s_asinh.c - $(CHEW) < $(srcdir)/s_asinh.c >$@ 2>/dev/null - touch stmp-def -satan.def: s_atan.c - $(CHEW) < $(srcdir)/s_atan.c >$@ 2>/dev/null - touch stmp-def -watan2.def: w_atan2.c - $(CHEW) < $(srcdir)/w_atan2.c >$@ 2>/dev/null - touch stmp-def -watanh.def: w_atanh.c - $(CHEW) < $(srcdir)/w_atanh.c >$@ 2>/dev/null - touch stmp-def -wj0.def: w_j0.c - $(CHEW) < $(srcdir)/w_j0.c >$@ 2>/dev/null - touch stmp-def -scopysign.def: s_copysign.c - $(CHEW) < $(srcdir)/../common/s_copysign.c >$@ 2>/dev/null - touch stmp-def -wcosh.def: w_cosh.c - $(CHEW) < $(srcdir)/w_cosh.c >$@ 2>/dev/null - touch stmp-def -serf.def: s_erf.c - $(CHEW) < $(srcdir)/s_erf.c >$@ 2>/dev/null - touch stmp-def -wexp.def: w_exp.c - $(CHEW) < $(srcdir)/w_exp.c >$@ 2>/dev/null - touch stmp-def -sfabs.def: s_fabs.c - $(CHEW) < $(srcdir)/s_fabs.c >$@ 2>/dev/null - touch stmp-def -sfloor.def: s_floor.c - $(CHEW) < $(srcdir)/s_floor.c >$@ 2>/dev/null - touch stmp-def -wfmod.def: w_fmod.c - $(CHEW) < $(srcdir)/w_fmod.c >$@ 2>/dev/null - touch stmp-def -sfrexp.def: s_frexp.c - $(CHEW) < $(srcdir)/s_frexp.c >$@ 2>/dev/null - touch stmp-def -wgamma.def: w_gamma.c - $(CHEW) < $(srcdir)/w_gamma.c >$@ 2>/dev/null - touch stmp-def -whypot.def: w_hypot.c - $(CHEW) < $(srcdir)/w_hypot.c >$@ 2>/dev/null - touch stmp-def -sldexp.def: s_ldexp.c - $(CHEW) < $(srcdir)/s_ldexp.c >$@ 2>/dev/null - touch stmp-def -wlog.def: w_log.c - $(CHEW) < $(srcdir)/w_log.c >$@ 2>/dev/null - touch stmp-def -wlog10.def: w_log10.c - $(CHEW) < $(srcdir)/w_log10.c >$@ 2>/dev/null - touch stmp-def -wpow.def: w_pow.c - $(CHEW) < $(srcdir)/w_pow.c >$@ 2>/dev/null - touch stmp-def -wremainder.def: w_remainder.c - $(CHEW) < $(srcdir)/w_remainder.c >$@ 2>/dev/null - touch stmp-def -ssin.def: s_sin.c - $(CHEW) < $(srcdir)/s_sin.c >$@ 2>/dev/null - touch stmp-def -wsinh.def: w_sinh.c - $(CHEW) < $(srcdir)/w_sinh.c >$@ 2>/dev/null - touch stmp-def -wsqrt.def: w_sqrt.c - $(CHEW) < $(srcdir)/w_sqrt.c >$@ 2>/dev/null - touch stmp-def -stan.def: s_tan.c - $(CHEW) < $(srcdir)/s_tan.c >$@ 2>/dev/null - touch stmp-def -stanh.def: s_tanh.c - $(CHEW) < $(srcdir)/s_tanh.c >$@ 2>/dev/null - touch stmp-def - # A partial dependency list. $(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h $(srcdir)/../common/fdlibm.h diff --git a/newlib/libm/math/Makefile.in b/newlib/libm/math/Makefile.in index 599623892..c0c3167f8 100644 --- a/newlib/libm/math/Makefile.in +++ b/newlib/libm/math/Makefile.in @@ -373,14 +373,14 @@ libmath_la_LDFLAGS = -Xcompiler -nostdlib @USE_LIBTOOL_FALSE@noinst_LIBRARIES = lib.a @USE_LIBTOOL_FALSE@lib_a_SOURCES = $(src) $(fsrc) @USE_LIBTOOL_FALSE@lib_a_CFLAGS = $(AM_CFLAGS) -chobj = wacos.def wacosh.def wasin.def sasinh.def \ - satan.def watan2.def watanh.def wj0.def \ - wcosh.def serf.def wexp.def \ - sfabs.def sfloor.def wfmod.def sfrexp.def \ - wgamma.def whypot.def sldexp.def wlog.def \ - wlog10.def \ - wpow.def wremainder.def ssin.def wsinh.def \ - wsqrt.def stan.def stanh.def +chobj = w_acos.def w_acosh.def w_asin.def s_asinh.def \ + s_atan.def w_atan2.def w_atanh.def w_j0.def \ + w_cosh.def s_erf.def w_exp.def w_exp2.def \ + s_fabs.def s_floor.def w_fmod.def s_frexp.def \ + w_gamma.def w_hypot.def s_ldexp.def w_log.def \ + w_log10.def \ + w_pow.def w_remainder.def s_sin.def w_sinh.def \ + w_sqrt.def s_tan.def s_tanh.def SUFFIXES = .def CHEW = ../../doc/makedoc -f $(srcdir)/../../doc/doc.str @@ -1424,94 +1424,6 @@ objectlist.awk.in: $(noinst_LTLIBRARIES) doc: $(chobj) cat $(srcdir)/math.tex >> $(TARGETDOC) -# Texinfo does not appear to support underscores in file names, so we -# name the .def files without underscores. - -wacos.def: w_acos.c - $(CHEW) < $(srcdir)/w_acos.c >$@ 2>/dev/null - touch stmp-def -wacosh.def: w_acosh.c - $(CHEW) < $(srcdir)/w_acosh.c >$@ 2>/dev/null - touch stmp-def -wasin.def: w_asin.c - $(CHEW) < $(srcdir)/w_asin.c >$@ 2>/dev/null - touch stmp-def -sasinh.def: s_asinh.c - $(CHEW) < $(srcdir)/s_asinh.c >$@ 2>/dev/null - touch stmp-def -satan.def: s_atan.c - $(CHEW) < $(srcdir)/s_atan.c >$@ 2>/dev/null - touch stmp-def -watan2.def: w_atan2.c - $(CHEW) < $(srcdir)/w_atan2.c >$@ 2>/dev/null - touch stmp-def -watanh.def: w_atanh.c - $(CHEW) < $(srcdir)/w_atanh.c >$@ 2>/dev/null - touch stmp-def -wj0.def: w_j0.c - $(CHEW) < $(srcdir)/w_j0.c >$@ 2>/dev/null - touch stmp-def -scopysign.def: s_copysign.c - $(CHEW) < $(srcdir)/../common/s_copysign.c >$@ 2>/dev/null - touch stmp-def -wcosh.def: w_cosh.c - $(CHEW) < $(srcdir)/w_cosh.c >$@ 2>/dev/null - touch stmp-def -serf.def: s_erf.c - $(CHEW) < $(srcdir)/s_erf.c >$@ 2>/dev/null - touch stmp-def -wexp.def: w_exp.c - $(CHEW) < $(srcdir)/w_exp.c >$@ 2>/dev/null - touch stmp-def -sfabs.def: s_fabs.c - $(CHEW) < $(srcdir)/s_fabs.c >$@ 2>/dev/null - touch stmp-def -sfloor.def: s_floor.c - $(CHEW) < $(srcdir)/s_floor.c >$@ 2>/dev/null - touch stmp-def -wfmod.def: w_fmod.c - $(CHEW) < $(srcdir)/w_fmod.c >$@ 2>/dev/null - touch stmp-def -sfrexp.def: s_frexp.c - $(CHEW) < $(srcdir)/s_frexp.c >$@ 2>/dev/null - touch stmp-def -wgamma.def: w_gamma.c - $(CHEW) < $(srcdir)/w_gamma.c >$@ 2>/dev/null - touch stmp-def -whypot.def: w_hypot.c - $(CHEW) < $(srcdir)/w_hypot.c >$@ 2>/dev/null - touch stmp-def -sldexp.def: s_ldexp.c - $(CHEW) < $(srcdir)/s_ldexp.c >$@ 2>/dev/null - touch stmp-def -wlog.def: w_log.c - $(CHEW) < $(srcdir)/w_log.c >$@ 2>/dev/null - touch stmp-def -wlog10.def: w_log10.c - $(CHEW) < $(srcdir)/w_log10.c >$@ 2>/dev/null - touch stmp-def -wpow.def: w_pow.c - $(CHEW) < $(srcdir)/w_pow.c >$@ 2>/dev/null - touch stmp-def -wremainder.def: w_remainder.c - $(CHEW) < $(srcdir)/w_remainder.c >$@ 2>/dev/null - touch stmp-def -ssin.def: s_sin.c - $(CHEW) < $(srcdir)/s_sin.c >$@ 2>/dev/null - touch stmp-def -wsinh.def: w_sinh.c - $(CHEW) < $(srcdir)/w_sinh.c >$@ 2>/dev/null - touch stmp-def -wsqrt.def: w_sqrt.c - $(CHEW) < $(srcdir)/w_sqrt.c >$@ 2>/dev/null - touch stmp-def -stan.def: s_tan.c - $(CHEW) < $(srcdir)/s_tan.c >$@ 2>/dev/null - touch stmp-def -stanh.def: s_tanh.c - $(CHEW) < $(srcdir)/s_tanh.c >$@ 2>/dev/null - touch stmp-def - # A partial dependency list. $(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h $(srcdir)/../common/fdlibm.h diff --git a/newlib/libm/math/math.tex b/newlib/libm/math/math.tex index 7bc56ce9a..608268ab9 100644 --- a/newlib/libm/math/math.tex +++ b/newlib/libm/math/math.tex @@ -9,6 +9,7 @@ Two definitions from @file{math.h} are of particular interest. @item The representation of infinity as a @code{double} is defined as @code{HUGE_VAL}; this number is returned on overflow by many functions. +The macro @code{HUGE_VALF} is a corresponding value for @code{float}. @item The structure @code{exception} is used when you write customized error @@ -44,43 +45,59 @@ machines---are available when you include @file{fastmath.h} instead of * atan:: Arctangent * atan2:: Arctangent of y/x * atanh:: Inverse hyperbolic tangent -* jN:: Bessel functions (jN, yN) +* jN:: Bessel functions (jN, yN) * cbrt:: Cube root * copysign:: Sign of Y, magnitude of X * cosh:: Hyperbolic cosine -* erf:: Error function (erf, erfc) -* exp:: Exponential -* expm1:: Exponential of x, - 1 +* erf:: Error function (erf, erfc) +* exp:: Exponential, base e +* exp2:: Exponential, base 2 +* expm1:: Exponential, base e, of x - 1 * fabs:: Absolute value (magnitude) -* floor:: Floor and ceiling (floor, ceil) +* fdim:: Positive difference +* floor:: Floor and ceiling (floor, ceil) +* fma:: Floating multiply add +* fmax:: Maximum +* fmin:: Minimum * fmod:: Floating-point remainder (modulo) +* fpclassify:: Floating-point classification macro * frexp:: Split floating-point number * gamma:: Logarithmic gamma function * hypot:: Distance from origin * ilogb:: Get exponent * infinity:: Floating infinity -* isnan:: Check type of number -* ldexp:: Load exponent +* isgreater:: Comparison macros +* ldexp:: Scale by a power of 2 * log:: Natural logarithms * log10:: Base 10 logarithms * log1p:: Log of 1 + X +* log2:: Base 2 logarithms +* logb:: Get exponent +* lrint:: Round to integer +* lround:: Round to integer, away from zero (lround, llround) * matherr:: Modifiable math error handler * modf:: Split fractional and integer parts * nan:: Floating Not a Number +* nearbyint:: Round to integer * nextafter:: Get next representable number * pow:: X to the power Y * remainder:: remainder of X divided by Y -* scalbn:: scalbn +* remquo:: Remainder and part of quotient +* rint:: Round to integer +* round:: Round to integer, away from zero +* scalbn:: Scale by a power of FLT_RADIX (2) +* signbit:: Does floating-point number have negative sign? * sin:: Sine or cosine (sin, cos) * sinh:: Hyperbolic sine * sqrt:: Positive square root * tan:: Tangent * tanh:: Hyperbolic tangent +* trunc:: Round to integer, towards zero @end menu @page @node version -@section Version of library +@section Error Handling There are four different versions of the math library routines: IEEE, POSIX, X/Open, or SVID. The version may be selected at runtime by @@ -114,119 +131,89 @@ log: DOMAIN error The library is set to X/Open mode by default. -@page -@include math/wacos.def - -@page -@include math/wacosh.def - -@page -@include math/wasin.def - -@page -@include math/sasinh.def - -@page -@include math/satan.def - -@page -@include math/watan2.def - -@page -@include math/watanh.def - -@page -@include math/wj0.def - -@page -@include common/scbrt.def - -@page -@include common/scopysign.def - -@page -@include math/wcosh.def - -@page -@include math/serf.def - -@page -@include math/wexp.def - -@page -@include common/sexpm1.def - -@page -@include math/sfabs.def - -@page -@include math/sfloor.def - -@page -@include math/wfmod.def - -@page -@include math/sfrexp.def - -@page -@include math/wgamma.def - -@page -@include math/whypot.def - -@page -@include common/silogb.def - -@page -@include common/sinfinity.def - -@page -@include common/sisnan.def - -@page -@include math/sldexp.def - -@page -@include math/wlog.def - -@page -@include math/wlog10.def - -@page -@include common/slog1p.def - -@page -@include common/smatherr.def - -@page -@include common/smodf.def - -@page -@include common/snan.def - -@page -@include common/snextafter.def - -@page -@include math/wpow.def - -@page -@include math/wremainder.def - -@page -@include common/sscalbn.def - -@page -@include math/wsqrt.def - -@page -@include math/ssin.def - -@page -@include math/wsinh.def - -@page -@include math/stan.def - -@page -@include math/stanh.def +The aforementioned error reporting is the supported Newlib libm error +handling method. However, the majority of the functions are written +so as to produce the floating-point exceptions (e.g. "invalid", +"divide-by-zero") as required by the C and POSIX standards, for +floating-point implementations that support them. Newlib does not provide +the floating-point exception access routines defined in the standards +for fenv.h, though, which is why they are considered unsupported. It is +mentioned in case you have separately-provided access routines so that +you are aware that they can be caused. + +@section Standards Compliance And Portability +Most of the individual function descriptions describe the standards to which +each function complies. However, these descriptions are mostly out of date, +having been written before C99 was released. One of these days we'll get +around to updating the rest of them. (If you'd like to help, please let us +know.) + +``C99'' refers to ISO/IEC 9899:1999, ``Programming languages--C''. +``POSIX'' refers to IEEE Standard 1003.1. POSIX@registeredsymbol{} is a +registered trademark of The IEEE. + +@c To sort the include list easily, keep the indentation right because want to +@c skip the s_|w_ at the start of most--but not all--of the file names. +@c (e.g., isgreater.def does not have a leading s nor w.) Then, sort +@c based on the column. For example: "sort -t@ -k3.17" +@c A few hand-edits might be appropriate after a sort, although not necessary +@c and are a nuisance as ought to be kept in sync with menu list above: +@c atan2 after atan, exp2 after exp, log first in log list, and w_j0 to place +@c to reflect function name of Bessel (as opposed to j; e.g. after atanh, +@c before cbrt). + +@page @include math/w_acos.def +@page @include math/w_acosh.def +@page @include math/w_asin.def +@page @include math/s_asinh.def +@page @include math/s_atan.def +@page @include math/w_atan2.def +@page @include math/w_atanh.def +@page @include math/w_j0.def +@page @include common/s_cbrt.def +@page @include common/s_copysign.def +@page @include math/w_cosh.def +@page @include math/s_erf.def +@page @include math/w_exp.def +@page @include math/w_exp2.def +@page @include common/s_expm1.def +@page @include math/s_fabs.def +@page @include common/s_fdim.def +@page @include math/s_floor.def +@page @include common/s_fma.def +@page @include common/s_fmax.def +@page @include common/s_fmin.def +@page @include math/w_fmod.def +@page @include math/s_frexp.def +@page @include math/w_gamma.def +@page @include math/w_hypot.def +@page @include common/s_ilogb.def +@page @include common/s_infinity.def +@page @include common/isgreater.def +@page @include common/s_isnan.def +@page @include math/s_ldexp.def +@page @include math/w_log.def +@page @include math/w_log10.def +@page @include common/s_log1p.def +@page @include common/s_log2.def +@page @include common/s_logb.def +@page @include common/s_lrint.def +@page @include common/s_lround.def +@page @include common/s_matherr.def +@page @include common/s_modf.def +@page @include common/s_nan.def +@page @include common/s_nearbyint.def +@page @include common/s_nextafter.def +@page @include math/w_pow.def +@page @include math/w_remainder.def +@page @include common/s_remquo.def +@page @include common/s_rint.def +@page @include common/s_round.def +@page @include common/s_scalbn.def +@page @include common/s_signbit.def +@page @include math/s_sin.def +@page @include math/w_sinh.def +@page @include math/w_sqrt.def +@page @include math/s_tan.def +@page @include math/s_tanh.def +@page @include common/s_trunc.def diff --git a/newlib/libm/math/w_exp2.c b/newlib/libm/math/w_exp2.c index ed0bc39e9..efb676df1 100644 --- a/newlib/libm/math/w_exp2.c +++ b/newlib/libm/math/w_exp2.c @@ -13,7 +13,7 @@ /* FUNCTION - <<exp2>>, <<exp2f>>---exponential + <<exp2>>, <<exp2f>>--exponential, base 2 INDEX exp2 INDEX @@ -24,14 +24,6 @@ ANSI_SYNOPSIS double exp2(double <[x]>); float exp2f(float <[x]>); -TRAD_SYNOPSIS - #include <math.h> - double exp2(<[x]>); - double <[x]>; - - float exp2f(<[x]>); - float <[x]>; - DESCRIPTION <<exp2>> and <<exp2f>> calculate 2 ^ <[x]>, that is, @ifnottex @@ -50,6 +42,9 @@ RETURNS result overflows, the returned value is <<HUGE_VAL>>. In either case, <<errno>> is set to <<ERANGE>>. +PORTABILITY + ANSI C, POSIX. + */ /* diff --git a/newlib/libm/math/w_gamma.c b/newlib/libm/math/w_gamma.c index fad40496d..3717f53e5 100644 --- a/newlib/libm/math/w_gamma.c +++ b/newlib/libm/math/w_gamma.c @@ -12,11 +12,20 @@ * */ +/* BUG: FIXME? + According to Linux man pages for tgamma, lgamma, and gamma, the gamma +function was originally defined in BSD as implemented here--the log of the gamma +function. BSD 4.3 changed the name to lgamma, apparently removing gamma. BSD +4.4 re-introduced the gamma name with the more intuitive, without logarithm, +plain gamma function. The C99 standard apparently wanted to avoid a problem +with the poorly-named earlier gamma and used tgamma when adding a plain +gamma function. + So the current gamma is matching an old, bad definition, and not +matching a newer, better definition. */ /* FUNCTION - <<gamma>>, <<gammaf>>, <<lgamma>>, <<lgammaf>>, <<gamma_r>>, - <<gammaf_r>>, <<lgamma_r>>, <<lgammaf_r>>---logarithmic gamma - function + <<gamma>>, <<gammaf>>, <<lgamma>>, <<lgammaf>>, <<gamma_r>>, <<gammaf_r>>, <<lgamma_r>>, <<lgammaf_r>>, <<tgamma>>, and <<tgammaf>>--logarithmic and plain gamma functions + INDEX gamma INDEX @@ -33,6 +42,10 @@ INDEX lgamma_r INDEX lgammaf_r +INDEX +tgamma +INDEX +tgammaf ANSI_SYNOPSIS #include <math.h> @@ -44,6 +57,8 @@ double gamma_r(double <[x]>, int *<[signgamp]>); float gammaf_r(float <[x]>, int *<[signgamp]>); double lgamma_r(double <[x]>, int *<[signgamp]>); float lgammaf_r(float <[x]>, int *<[signgamp]>); +double tgamma(double <[x]>); +float tgammaf(float <[x]>); TRAD_SYNOPSIS #include <math.h> @@ -67,15 +82,19 @@ int <[signgamp]>; float lgammaf_r(<[x]>, <[signgamp]>) float <[x]>; int <[signgamp]>; +double tgamma(<[x]>) +double <[x]>; +float tgammaf(<[x]>) +float <[x]>; DESCRIPTION -<<gamma>> calculates +<<gamma>> calculates @tex -$\mit ln\bigl(\Gamma(x)\bigr)$, +$\mit ln\bigl(\Gamma(x)\bigr)$, @end tex the natural logarithm of the gamma function of <[x]>. The gamma function (<<exp(gamma(<[x]>))>>) is a generalization of factorial, and retains -the property that +the property that @ifnottex <<exp(gamma(N))>> is equivalent to <<N*exp(gamma(N-1))>>. @end ifnottex @@ -83,13 +102,13 @@ the property that $\mit \Gamma(N)\equiv N\times\Gamma(N-1)$. @end tex Accordingly, the results of the gamma function itself grow very -quickly. <<gamma>> is defined as +quickly. <<gamma>> is defined as @tex $\mit ln\bigl(\Gamma(x)\bigr)$ rather than simply $\mit \Gamma(x)$ @end tex @ifnottex the natural log of the gamma function, rather than the gamma function -itself, +itself, @end ifnottex to extend the useful range of results representable. @@ -113,8 +132,17 @@ variable <<signgam>> is not used. These functions may be used for reentrant calls (but they will still set the global variable <<errno>> if an error occurs). +<<tgamma>> and <<tgammaf>> are the "true gamma" functions, returning +@tex +$\mit \Gamma(x)$, +@end tex +the gamma function of <[x]>--without a logarithm. +(They are apparently so named because of the prior existence of the old, +poorly-named <<gamma>> functions which returned the log of gamma up +through BSD 4.2.) + RETURNS -Normally, the computed result is returned. +Normally, the computed result is returned. When <[x]> is a nonpositive integer, <<gamma>> returns <<HUGE_VAL>> and <<errno>> is set to <<EDOM>>. If the result overflows, <<gamma>> @@ -123,7 +151,12 @@ returns <<HUGE_VAL>> and <<errno>> is set to <<ERANGE>>. You can modify this error treatment using <<matherr>>. PORTABILITY -Neither <<gamma>> nor <<gammaf>> is ANSI C. */ +Neither <<gamma>> nor <<gammaf>> is ANSI C. It is better not to use either +of these; use <<lgamma>> or <<tgamma>> instead.@* +<<lgamma>>, <<lgammaf>>, <<tgamma>>, and <<tgammaf>> are nominally C standard +in terms of the base return values, although the <<matherr>> error-handling +is not standard, nor is the <[signgam]> global for <<lgamma>>. +*/ /* double gamma(double x) * Return the logarithm of the Gamma function of x. |