From 139f923bb4564fb021dce2b7ae460a3d8e3bb87d Mon Sep 17 00:00:00 2001 From: Jeff Johnston Date: Wed, 25 Mar 2009 19:13:24 +0000 Subject: 2009-03-25 Craig Howland * 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. --- newlib/libm/common/Makefile.am | 81 ++------------- newlib/libm/common/Makefile.in | 116 ++++++--------------- newlib/libm/common/isgreater.c | 75 ++++++++++++++ newlib/libm/common/s_fdim.c | 37 ++++++- newlib/libm/common/s_fma.c | 38 +++++++ newlib/libm/common/s_fmax.c | 25 +++++ newlib/libm/common/s_fmin.c | 25 +++++ newlib/libm/common/s_infinity.c | 15 +-- newlib/libm/common/s_isnan.c | 113 ++++++++++++++++---- newlib/libm/common/s_llround.c | 68 ++++++++++++ newlib/libm/common/s_log2.c | 46 +++++--- newlib/libm/common/s_logb.c | 104 +++++++++++++++---- newlib/libm/common/s_lrint.c | 38 +++++++ newlib/libm/common/s_lround.c | 38 +++++++ newlib/libm/common/s_nearbyint.c | 37 +++++++ newlib/libm/common/s_remquo.c | 219 ++++++++++++++++++++++++++++++++++----- newlib/libm/common/s_rint.c | 34 ++++++ newlib/libm/common/s_round.c | 32 ++++++ newlib/libm/common/s_scalbn.c | 44 ++++---- newlib/libm/common/s_signbit.c | 29 ++++++ newlib/libm/common/s_trunc.c | 29 ++++++ newlib/libm/common/sf_fdim.c | 6 +- newlib/libm/common/sf_fma.c | 10 +- newlib/libm/common/sf_llround.c | 55 ++++++++++ newlib/libm/common/sf_logb.c | 46 +++++--- newlib/libm/common/sf_lrint.c | 2 +- newlib/libm/common/sf_lround.c | 2 +- newlib/libm/common/sf_remquo.c | 150 ++++++++++++++++++++------- 28 files changed, 1195 insertions(+), 319 deletions(-) create mode 100644 newlib/libm/common/isgreater.c create mode 100644 newlib/libm/common/s_llround.c create mode 100644 newlib/libm/common/sf_llround.c (limited to 'newlib/libm/common') 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 +<>, <>, <>, <>, <>, and <>--comparison macros +INDEX + isgreater +INDEX + isgreaterequal +INDEX + isless +INDEX + islessequal +INDEX + islessgreater +INDEX + isunordered + +ANSI_SYNOPSIS + #include + 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 +<>, <>, <>, <>, +<>, and <> 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 <> macro returns the value of (x) > (y).@* +The <> macro returns the value of (x) >= (y).@* +The <> macro returns the value of (x) < (y).@* +The <> macro returns the value of (x) <= (y).@* +The <> macro returns the value of (x) < (y) || (x) > (y).@* +The <> 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 +<>, <>--positive difference +INDEX + fdim +INDEX + fdimf + +ANSI_SYNOPSIS + #include + double fdim(double <[x]>, double <[y]>); + float fdimf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> 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 <> 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 +<>, <>--floating multiply add +INDEX + fma +INDEX + fmaf + +ANSI_SYNOPSIS + #include + double fma(double <[x]>, double <[y]>, double <[z]>); + float fmaf(float <[x]>, float <[y]>, float <[z]>); + +DESCRIPTION +The <> 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 <> 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 <> 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 +<>, <>--maximum +INDEX + fmax +INDEX + fmaxf + +ANSI_SYNOPSIS + #include + double fmax(double <[x]>, double <[y]>); + float fmaxf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> 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 <> functions choose the numeric value. + +RETURNS +The <> 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 +<>, <>--minimum +INDEX + fmin +INDEX + fminf + +ANSI_SYNOPSIS + #include + double fmin(double <[x]>, double <[y]>); + float fminf(float <[x]>, float <[y]>); + +DESCRIPTION +The <> 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 <> functions choose the numeric value. + +RETURNS +The <> 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 - <>, <>---representation of infinity + <>, <>--representation of infinity INDEX infinity @@ -17,17 +17,18 @@ ANSI_SYNOPSIS double infinity(void); float infinityf(void); -TRAD_SYNOPSIS - #include - double infinity(); - float infinityf(); - - DESCRIPTION <> and <> return the special number IEEE infinity in double- and single-precision arithmetic respectively. +PORTABILITY +<> and <> 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 - <>, <>, <>, <>, <>, <>---test for exceptional numbers +<>, <>, <>, <>, and <>--floating-point classification macros; <>, <>, <>, <>, <>, <>--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 + 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 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 - 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 +<>, <>, <>, <>, and <> 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 <> 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 <> 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 "<>" set of macros provide a useful set of shorthand ways for +classifying floating-point numbers, providing the following equivalent +relations: + +o+ +o <>(<[x]>) +returns non-zero if <[x]> is finite. (It is equivalent to +(<>(<[x]>) != FP_INFINITE && <>(<[x]>) != FP_NAN).) + +o <>(<[x]>) +returns non-zero if <[x]> is infinite. (It is equivalent to +(<>(<[x]>) == FP_INFINITE).) + +o <>(<[x]>) +returns non-zero if <[x]> is NaN. (It is equivalent to +(<>(<[x]>) == FP_NAN).) + +o <>(<[x]>) +returns non-zero if <[x]> is normal. (It is equivalent to +(<>(<[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 <> are macros that operate on multiple types of floating-point. The SUSv2 standard declares <> 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 - <>, <>---base 2 logarithm + <>, <>--base 2 logarithm INDEX log2 INDEX @@ -25,28 +25,40 @@ ANSI_SYNOPSIS double log2(double <[x]>); float log2f(float <[x]>); -TRAD_SYNOPSIS - #include - double log2(<[x]>); - double <[x]>; - - float log2f(<[x]>); - float <[x]>; - DESCRIPTION - <> returns the base 2 logarithm of <[x]>. +The <> 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 <>. (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. - <> is identical, save that it takes and returns <> values. +You can use the (non-ANSI) function <> to specify error +handling for these functions, indirectly through the respective <> +function. RETURNS - On success, <> and <> return the calculated value. - If the result underflows, the returned value is <<0>>. If the - result overflows, the returned value is <>. In - either case, <> is set to <>. +The <> functions return +@ifnottex +<)>> +@end ifnottex +@tex +$log_2(x)$ +@end tex +on success. +When <[x]> is zero, the +returned value is <<-HUGE_VAL>> and <> is set to <>. +When <[x]> is negative, the returned value is NaN (not a number) and +<> is set to <>. You can control the error behavior via +<>. PORTABILITY - <> and <> 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 + <>, <>--get exponent of floating-point number +INDEX + logb +INDEX + logbf + +ANSI_SYNOPSIS + #include + double logb(double <[x]>); + float logbf(float <[x]>); + +DESCRIPTION +The <> 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]>, <>(<[x]>) = <>(<>(<>(<[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 <> functions examine the argument <[x]>, and return <[p]>. +The <> functions are similar to the <> 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 <> 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 +<>, <>, <>, <>--round to integer +INDEX + lrint +INDEX + lrintf +INDEX + llrint +INDEX + llrintf + +ANSI_SYNOPSIS + #include + long int lrint(double <[x]>); + long int lrintf(float <[x]>); + long long int llrint(double <[x]>); + long long int llrintf(float <[x]>); + +DESCRIPTION +The <> and <> 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 +<> + +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 +<>, <>, <>, <>--round to integer, to nearest +INDEX + lround +INDEX + lroundf +INDEX + llround +INDEX + llroundf + +ANSI_SYNOPSIS + #include + long int lround(double <[x]>); + long int lroundf(float <[x]>); + long long int llround(double <[x]>); + long long int llroundf(float <[x]>); + +DESCRIPTION + The <> and <> 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 <> functions for the return being the same floating-point type +as the argument. <>, <>. + +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 +<>, <>--round to integer +INDEX + nearbyint +INDEX + nearbyintf + +ANSI_SYNOPSIS + #include + double nearbyint(double <[x]>); + float nearbyintf(float <[x]>); + +DESCRIPTION +The <> 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 <> 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 <> functions are identical to the <> +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 +<>, <> +*/ #include #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 +<>, <>--remainder and part of quotient +INDEX + remquo +INDEX + remquof -#include "fdlibm.h" +ANSI_SYNOPSIS + #include + double remquo(double <[x]>, double <[y]>, int *<[quo]>); + float remquof(float <[x]>, float <[y]>, int *<[quo]>); + +DESCRIPTION +The <> functions compute the same remainder as the <> +functions; this value is in the range -<[y]>/2 ... +<[y]>/2. In the object +pointed to by <> they store a value whose sign is the sign of <>/<> +and whose magnitude is congruent modulo 2**n to the magnitude of the integral +quotient of <>/<>. (That is, <> 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 <> functions return <[x]> REM <[y]>. -#ifndef _DOUBLE_IS_32BITS +BUGS +IEEE754-2008 calls for <>(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 +#include +#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>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<>(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<>(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>31); lx = lx+lx;} + else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;} + q <<= 1; + } + hz=hx-hy;lz=lx-ly; if(lx=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 +<>, <>--round to integer +INDEX + rint +INDEX + rintf + +ANSI_SYNOPSIS + #include + double rint(double <[x]>); + float rintf(float <[x]>); + +DESCRIPTION + The <> 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 <> functions for the + same function with the "inexact" floating-point exception never being + raised. Newlib does not directly support floating-point exceptions. + The <> 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 +<>, <> + +*/ /* * 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 to integer, to nearest +INDEX + round +INDEX + roundf + +ANSI_SYNOPSIS + #include + double round(double <[x]>); + float roundf(float <[x]>); + +DESCRIPTION + The <> 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 + <> 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 +<>, <> + +*/ #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 -<>, <>---scale by power of two +<>, <>, <>, <>--scale by power of FLT_RADIX (=2) INDEX scalbn INDEX scalbnf +INDEX + scalbln +INDEX + scalblnf ANSI_SYNOPSIS #include - double scalbn(double <[x]>, int <[y]>); - float scalbnf(float <[x]>, int <[y]>); - -TRAD_SYNOPSIS - #include - 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 -<> and <> 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 <> and <> 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 <> +functions equivalent to the <> functions. RETURNS -<[x]> times 2 to the power <[n]>. +<[x]> times 2 to the power <[n]>. A range error may occur. PORTABILITY -Neither <> nor <> is required by ANSI C or by the System V -Interface Definition (Issue 2). +ANSI C, POSIX + +SEEALSO +<> */ 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 +<>--Does floating-point number have negative sign? + +INDEX + signbit + +ANSI_SYNOPSIS + #include + int signbit(real-floating <[x]>); + +DESCRIPTION +The <> 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 <> in all cases. + +RETURNS +The <> 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 +<>, <>--round to integer, towards zero +INDEX + trunc +INDEX + truncf + +ANSI_SYNOPSIS + #include + double trunc(double <[x]>); + float truncf(float <[x]>); + +DESCRIPTION + The <> 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 <> 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 #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>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; +} -- cgit v1.2.3