+2011-02-17 Jakub Jelinek <jakub@redhat.com>
+
+ * libquadmath.texi (FLT128_DIG, FLT128_MIN_10_EXP,
+ FLT128_MAX_10_EXP): Document.
+ (strtoflt128): Remove obsolete comment.
+ * configure.ac (HAVE_STRTOULL): New check.
+ * printf/gmp-impl.h (mpn_construct_float128): New prototype,
+ define.
+ * printf/mul_n.c: Include <config.h>.
+ * printf/add_n.c: Likewise.
+ * printf/cmp.c: Likewise.
+ * printf/fpioconst.c: Likewise.
+ * printf/mul_1.c: Likewise.
+ * printf/rshift.c: Likewise.
+ * printf/lshift.c: Likewise.
+ * printf/submul_1.c: Likewise.
+ * printf/sub_n.c: Likewise.
+ * printf/divrem.c: Likewise.
+ * printf/addmul_1.c: Likewise.
+ * printf/mul.c: Likewise.
+ * printf/quadmath-printf.h (isupper, isdigit, tolower): Change
+ to avoid evaluating argument multiple times.
+ (isxdigit): Redefine.
+ * strtod/strtoflt128.c: New file.
+ * strtod/strtod_l.c: New file.
+ * strtod/mpn2flt128.c: New file.
+ * strtod/grouping.h: New file.
+ * strtod/tens_in_limb.c: New file.
+ * gdtoa/arith.h: Removed.
+ * gdtoa/gd_qnan.h: Removed.
+ * gdtoa/gdtoa_fltrnds.h: Removed.
+ * gdtoa/gdtoa.h: Removed.
+ * gdtoa/gdtoaimp.h: Removed.
+ * gdtoa/gethex.c: Removed.
+ * gdtoa/gmisc.c: Removed.
+ * gdtoa/hd_init.c: Removed.
+ * gdtoa/hexnan.c: Removed.
+ * gdtoa/makefile: Removed.
+ * gdtoa/misc.c: Removed.
+ * gdtoa/README.gdtoa: Removed.
+ * gdtoa/smisc.c: Removed.
+ * gdtoa/strtodg.c: Removed.
+ * gdtoa/strtopQ.c: Removed.
+ * gdtoa/sum.c: Removed.
+ * quadmath.h (FLT128_DIG, FLT128_MIN_10_EXP, FLT128_MAX_10_EXP):
+ Define.
+ * Makefile.am (libquadmath_la_SOURCES): Remove gdtoa/*, add
+ strtod/strtoflt128.c, strtod/mpn2flt128.c and strtod/tens_in_limb.c.
+ * config.h.in: Regenerated.
+ * configure: Regenerated.
+ * Makefile.in: Regenerated.
+
2011-02-16 Jakub Jelinek <jakub@redhat.com>
* printf/quadmath-printf.c: Also check __GLIBC__ when checking
libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
libquadmath_la_SOURCES = \
- gdtoa/arith.h gdtoa/gdtoa_fltrnds.h gdtoa/gd_qnan.h gdtoa/gdtoaimp.h \
- gdtoa/gdtoa.h quadmath-imp.h \
- gdtoa/hd_init.c gdtoa/smisc.c gdtoa/sum.c \
- gdtoa/gethex.c gdtoa/hexnan.c gdtoa/strtodg.c \
- gdtoa/gmisc.c gdtoa/misc.c gdtoa/strtopQ.c \
math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \
math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
printf/addmul_1.c printf/add_n.c printf/cmp.c printf/divrem.c \
printf/flt1282mpn.c printf/fpioconst.c printf/lshift.c printf/mul_1.c \
printf/mul_n.c printf/mul.c printf/printf_fphex.c printf/printf_fp.c \
- printf/quadmath-printf.c printf/rshift.c printf/submul_1.c printf/sub_n.c
+ printf/quadmath-printf.c printf/rshift.c printf/submul_1.c printf/sub_n.c \
+ strtod/strtoflt128.c strtod/mpn2flt128.c strtod/tens_in_limb.c
# Work around what appears to be a GNU make bug handling MAKEFLAGS
"$(DESTDIR)$(libsubincludedir)"
LTLIBRARIES = $(toolexeclib_LTLIBRARIES)
am__dirstamp = $(am__leading_dot)dirstamp
-@BUILD_LIBQUADMATH_TRUE@am_libquadmath_la_OBJECTS = gdtoa/hd_init.lo \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/smisc.lo gdtoa/sum.lo \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/gethex.lo gdtoa/hexnan.lo \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/strtodg.lo gdtoa/gmisc.lo \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/misc.lo gdtoa/strtopQ.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/acoshq.lo math/fmodq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/acosq.lo math/frexpq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/rem_pio2q.lo math/asinhq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/hypotq.lo math/remainderq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/asinq.lo math/rintq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/atan2q.lo math/isinfq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/roundq.lo math/atanhq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/isnanq.lo math/scalblnq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/atanq.lo math/j0q.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/scalbnq.lo math/cbrtq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/j1q.lo math/signbitq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/ceilq.lo math/jnq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/sincos_table.lo math/complex.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/ldexpq.lo math/sincosq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/copysignq.lo math/lgammaq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/sincosq_kernel.lo math/coshq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/llroundq.lo math/sinhq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/cosq.lo math/log10q.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/sinq.lo math/cosq_kernel.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/log1pq.lo math/sinq_kernel.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/erfq.lo math/logq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/sqrtq.lo math/expm1q.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/lroundq.lo math/tanhq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/expq.lo math/modfq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/tanq.lo math/fabsq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/nanq.lo math/tgammaq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/finiteq.lo math/nextafterq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/truncq.lo math/floorq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/powq.lo math/fmaq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/cacoshq.lo math/cacosq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/casinhq.lo math/casinq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/catanhq.lo math/catanq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/cimagq.lo math/conjq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/cprojq.lo math/crealq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/fdimq.lo math/fmaxq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/fminq.lo math/ilogbq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/llrintq.lo math/log2q.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/lrintq.lo math/nearbyintq.lo \
-@BUILD_LIBQUADMATH_TRUE@ math/remquoq.lo printf/addmul_1.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/add_n.lo printf/cmp.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/divrem.lo printf/flt1282mpn.lo \
+@BUILD_LIBQUADMATH_TRUE@am_libquadmath_la_OBJECTS = math/acoshq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/fmodq.lo math/acosq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/frexpq.lo math/rem_pio2q.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/asinhq.lo math/hypotq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/remainderq.lo math/asinq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/rintq.lo math/atan2q.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/isinfq.lo math/roundq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/atanhq.lo math/isnanq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/scalblnq.lo math/atanq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/j0q.lo math/scalbnq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/cbrtq.lo math/j1q.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/signbitq.lo math/ceilq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/jnq.lo math/sincos_table.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/complex.lo math/ldexpq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/sincosq.lo math/copysignq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/lgammaq.lo math/sincosq_kernel.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/coshq.lo math/llroundq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/sinhq.lo math/cosq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/log10q.lo math/sinq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/cosq_kernel.lo math/log1pq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/sinq_kernel.lo math/erfq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/logq.lo math/sqrtq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/expm1q.lo math/lroundq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/tanhq.lo math/expq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/modfq.lo math/tanq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/fabsq.lo math/nanq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/tgammaq.lo math/finiteq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/nextafterq.lo math/truncq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/floorq.lo math/powq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/fmaq.lo math/cacoshq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/cacosq.lo math/casinhq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/casinq.lo math/catanhq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/catanq.lo math/cimagq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/conjq.lo math/cprojq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/crealq.lo math/fdimq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/fmaxq.lo math/fminq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/ilogbq.lo math/llrintq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/log2q.lo math/lrintq.lo \
+@BUILD_LIBQUADMATH_TRUE@ math/nearbyintq.lo math/remquoq.lo \
+@BUILD_LIBQUADMATH_TRUE@ printf/addmul_1.lo printf/add_n.lo \
+@BUILD_LIBQUADMATH_TRUE@ printf/cmp.lo printf/divrem.lo \
+@BUILD_LIBQUADMATH_TRUE@ printf/flt1282mpn.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/fpioconst.lo printf/lshift.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/mul_1.lo printf/mul_n.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/mul.lo printf/printf_fphex.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/printf_fp.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/quadmath-printf.lo \
@BUILD_LIBQUADMATH_TRUE@ printf/rshift.lo printf/submul_1.lo \
-@BUILD_LIBQUADMATH_TRUE@ printf/sub_n.lo
+@BUILD_LIBQUADMATH_TRUE@ printf/sub_n.lo strtod/strtoflt128.lo \
+@BUILD_LIBQUADMATH_TRUE@ strtod/mpn2flt128.lo \
+@BUILD_LIBQUADMATH_TRUE@ strtod/tens_in_limb.lo
libquadmath_la_OBJECTS = $(am_libquadmath_la_OBJECTS)
libquadmath_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
@BUILD_LIBQUADMATH_TRUE@nodist_libsubinclude_HEADERS = quadmath.h quadmath_weak.h
@BUILD_LIBQUADMATH_TRUE@libsubincludedir = $(libdir)/gcc/$(target_alias)/$(gcc_version)/include
@BUILD_LIBQUADMATH_TRUE@libquadmath_la_SOURCES = \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/arith.h gdtoa/gdtoa_fltrnds.h gdtoa/gd_qnan.h gdtoa/gdtoaimp.h \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/gdtoa.h quadmath-imp.h \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/hd_init.c gdtoa/smisc.c gdtoa/sum.c \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/gethex.c gdtoa/hexnan.c gdtoa/strtodg.c \
-@BUILD_LIBQUADMATH_TRUE@ gdtoa/gmisc.c gdtoa/misc.c gdtoa/strtopQ.c \
@BUILD_LIBQUADMATH_TRUE@ math/acoshq.c math/fmodq.c math/acosq.c math/frexpq.c \
@BUILD_LIBQUADMATH_TRUE@ math/rem_pio2q.c math/asinhq.c math/hypotq.c math/remainderq.c \
@BUILD_LIBQUADMATH_TRUE@ math/asinq.c math/rintq.c math/atan2q.c math/isinfq.c \
@BUILD_LIBQUADMATH_TRUE@ printf/addmul_1.c printf/add_n.c printf/cmp.c printf/divrem.c \
@BUILD_LIBQUADMATH_TRUE@ printf/flt1282mpn.c printf/fpioconst.c printf/lshift.c printf/mul_1.c \
@BUILD_LIBQUADMATH_TRUE@ printf/mul_n.c printf/mul.c printf/printf_fphex.c printf/printf_fp.c \
-@BUILD_LIBQUADMATH_TRUE@ printf/quadmath-printf.c printf/rshift.c printf/submul_1.c printf/sub_n.c
+@BUILD_LIBQUADMATH_TRUE@ printf/quadmath-printf.c printf/rshift.c printf/submul_1.c printf/sub_n.c \
+@BUILD_LIBQUADMATH_TRUE@ strtod/strtoflt128.c strtod/mpn2flt128.c strtod/tens_in_limb.c
# Work around what appears to be a GNU make bug handling MAKEFLAGS
echo "rm -f \"$${dir}/so_locations\""; \
rm -f "$${dir}/so_locations"; \
done
-gdtoa/$(am__dirstamp):
- @$(MKDIR_P) gdtoa
- @: > gdtoa/$(am__dirstamp)
-gdtoa/$(DEPDIR)/$(am__dirstamp):
- @$(MKDIR_P) gdtoa/$(DEPDIR)
- @: > gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/hd_init.lo: gdtoa/$(am__dirstamp) \
- gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/smisc.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/sum.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/gethex.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/hexnan.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/strtodg.lo: gdtoa/$(am__dirstamp) \
- gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/gmisc.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/misc.lo: gdtoa/$(am__dirstamp) gdtoa/$(DEPDIR)/$(am__dirstamp)
-gdtoa/strtopQ.lo: gdtoa/$(am__dirstamp) \
- gdtoa/$(DEPDIR)/$(am__dirstamp)
math/$(am__dirstamp):
@$(MKDIR_P) math
@: > math/$(am__dirstamp)
printf/$(DEPDIR)/$(am__dirstamp)
printf/sub_n.lo: printf/$(am__dirstamp) \
printf/$(DEPDIR)/$(am__dirstamp)
+strtod/$(am__dirstamp):
+ @$(MKDIR_P) strtod
+ @: > strtod/$(am__dirstamp)
+strtod/$(DEPDIR)/$(am__dirstamp):
+ @$(MKDIR_P) strtod/$(DEPDIR)
+ @: > strtod/$(DEPDIR)/$(am__dirstamp)
+strtod/strtoflt128.lo: strtod/$(am__dirstamp) \
+ strtod/$(DEPDIR)/$(am__dirstamp)
+strtod/mpn2flt128.lo: strtod/$(am__dirstamp) \
+ strtod/$(DEPDIR)/$(am__dirstamp)
+strtod/tens_in_limb.lo: strtod/$(am__dirstamp) \
+ strtod/$(DEPDIR)/$(am__dirstamp)
libquadmath.la: $(libquadmath_la_OBJECTS) $(libquadmath_la_DEPENDENCIES)
$(libquadmath_la_LINK) $(am_libquadmath_la_rpath) $(libquadmath_la_OBJECTS) $(libquadmath_la_LIBADD) $(LIBS)
mostlyclean-compile:
-rm -f *.$(OBJEXT)
- -rm -f gdtoa/gethex.$(OBJEXT)
- -rm -f gdtoa/gethex.lo
- -rm -f gdtoa/gmisc.$(OBJEXT)
- -rm -f gdtoa/gmisc.lo
- -rm -f gdtoa/hd_init.$(OBJEXT)
- -rm -f gdtoa/hd_init.lo
- -rm -f gdtoa/hexnan.$(OBJEXT)
- -rm -f gdtoa/hexnan.lo
- -rm -f gdtoa/misc.$(OBJEXT)
- -rm -f gdtoa/misc.lo
- -rm -f gdtoa/smisc.$(OBJEXT)
- -rm -f gdtoa/smisc.lo
- -rm -f gdtoa/strtodg.$(OBJEXT)
- -rm -f gdtoa/strtodg.lo
- -rm -f gdtoa/strtopQ.$(OBJEXT)
- -rm -f gdtoa/strtopQ.lo
- -rm -f gdtoa/sum.$(OBJEXT)
- -rm -f gdtoa/sum.lo
-rm -f math/acoshq.$(OBJEXT)
-rm -f math/acoshq.lo
-rm -f math/acosq.$(OBJEXT)
-rm -f printf/sub_n.lo
-rm -f printf/submul_1.$(OBJEXT)
-rm -f printf/submul_1.lo
+ -rm -f strtod/mpn2flt128.$(OBJEXT)
+ -rm -f strtod/mpn2flt128.lo
+ -rm -f strtod/strtoflt128.$(OBJEXT)
+ -rm -f strtod/strtoflt128.lo
+ -rm -f strtod/tens_in_limb.$(OBJEXT)
+ -rm -f strtod/tens_in_limb.lo
distclean-compile:
-rm -f *.tab.c
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/gethex.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/gmisc.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/hd_init.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/hexnan.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/misc.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/smisc.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/strtodg.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/strtopQ.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@gdtoa/$(DEPDIR)/sum.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/acoshq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/acosq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@math/$(DEPDIR)/asinhq.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/rshift.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/sub_n.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@printf/$(DEPDIR)/submul_1.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@strtod/$(DEPDIR)/mpn2flt128.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@strtod/$(DEPDIR)/strtoflt128.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@strtod/$(DEPDIR)/tens_in_limb.Plo@am__quote@
.c.o:
@am__fastdepCC_TRUE@ depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
clean-libtool:
-rm -rf .libs _libs
- -rm -rf gdtoa/.libs gdtoa/_libs
-rm -rf math/.libs math/_libs
-rm -rf printf/.libs printf/_libs
+ -rm -rf strtod/.libs strtod/_libs
distclean-libtool:
-rm -f libtool config.lt
distclean-generic:
-test -z "$(CONFIG_CLEAN_FILES)" || rm -f $(CONFIG_CLEAN_FILES)
-test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES)
- -rm -f gdtoa/$(DEPDIR)/$(am__dirstamp)
- -rm -f gdtoa/$(am__dirstamp)
-rm -f math/$(DEPDIR)/$(am__dirstamp)
-rm -f math/$(am__dirstamp)
-rm -f printf/$(DEPDIR)/$(am__dirstamp)
-rm -f printf/$(am__dirstamp)
+ -rm -f strtod/$(DEPDIR)/$(am__dirstamp)
+ -rm -f strtod/$(am__dirstamp)
maintainer-clean-generic:
@echo "This command is intended for maintainers to use"
distclean: distclean-am distclean-multi
-rm -f $(am__CONFIG_DISTCLEAN_FILES)
- -rm -rf gdtoa/$(DEPDIR) math/$(DEPDIR) printf/$(DEPDIR)
+ -rm -rf math/$(DEPDIR) printf/$(DEPDIR) strtod/$(DEPDIR)
-rm -f Makefile
distclean-am: clean-am distclean-compile distclean-generic \
distclean-hdr distclean-libtool distclean-tags
maintainer-clean: maintainer-clean-am maintainer-clean-multi
-rm -f $(am__CONFIG_DISTCLEAN_FILES)
-rm -rf $(top_srcdir)/autom4te.cache
- -rm -rf gdtoa/$(DEPDIR) math/$(DEPDIR) printf/$(DEPDIR)
+ -rm -rf math/$(DEPDIR) printf/$(DEPDIR) strtod/$(DEPDIR)
-rm -f Makefile
maintainer-clean-am: distclean-am maintainer-clean-aminfo \
maintainer-clean-generic
/* Define to 1 if you have the <string.h> header file. */
#undef HAVE_STRING_H
+/* Define to 1 if you have the `strtoull' function. */
+#undef HAVE_STRTOULL
+
/* Define to 1 if you have the <sys/stat.h> header file. */
#undef HAVE_SYS_STAT_H
fi
+ for ac_func in strtoull
+do :
+ ac_fn_c_check_func "$LINENO" "strtoull" "ac_cv_func_strtoull"
+if test "x$ac_cv_func_strtoull" = x""yes; then :
+ cat >>confdefs.h <<_ACEOF
+#define HAVE_STRTOULL 1
+_ACEOF
+
+fi
+done
+
else
if test "x$ac_cv_lib_m_sqrtl" = x""yes; then
AC_CHECK_LIB([m],[feupdateenv],[AC_DEFINE([HAVE_FEUPDATEENV],[1],[libm includes feupdateenv])])
AC_CHECK_LIB([m],[fesetenv],[AC_DEFINE([HAVE_FESETENV],[1],[libm includes fesetenv])])
AC_CHECK_LIB([m],[fetestexcept],[AC_DEFINE([HAVE_FETESTEXCEPT],[1],[libm includes fetestexcept])])
+ AC_CHECK_FUNCS(strtoull)
else
if test "x$ac_cv_lib_m_sqrtl" = x""yes; then
AC_DEFINE([HAVE_SQRTL],[1],[libm includes sqrtl])
+++ /dev/null
-The content below is the README file of the gdtoa distribution, available
-from http://www.netlib.org/fp/
-
-----------------------------------------------------
-
-This directory contains source for a library of binary -> decimal
-and decimal -> binary conversion routines, for single-, double-,
-and extended-precision IEEE binary floating-point arithmetic, and
-other IEEE-like binary floating-point, including "double double",
-as in
-
- T. J. Dekker, "A Floating-Point Technique for Extending the
- Available Precision", Numer. Math. 18 (1971), pp. 224-242
-
-and
-
- "Inside Macintosh: PowerPC Numerics", Addison-Wesley, 1994
-
-The conversion routines use double-precision floating-point arithmetic
-and, where necessary, high precision integer arithmetic. The routines
-are generalizations of the strtod and dtoa routines described in
-
- David M. Gay, "Correctly Rounded Binary-Decimal and
- Decimal-Binary Conversions", Numerical Analysis Manuscript
- No. 90-10, Bell Labs, Murray Hill, 1990;
- http://cm.bell-labs.com/cm/cs/what/ampl/REFS/rounding.ps.gz
-
-(based in part on papers by Clinger and Steele & White: see the
-references in the above paper).
-
-The present conversion routines should be able to use any of IEEE binary,
-VAX, or IBM-mainframe double-precision arithmetic internally, but I (dmg)
-have so far only had a chance to test them with IEEE double precision
-arithmetic.
-
-The core conversion routines are strtodg for decimal -> binary conversions
-and gdtoa for binary -> decimal conversions. These routines operate
-on arrays of unsigned 32-bit integers of type ULong, a signed 32-bit
-exponent of type Long, and arithmetic characteristics described in
-struct FPI; FPI, Long, and ULong are defined in gdtoa.h. File arith.h
-is supposed to provide #defines that cause gdtoa.h to define its
-types correctly. File arithchk.c is source for a program that
-generates a suitable arith.h on all systems where I've been able to
-test it.
-
-The core conversion routines are meant to be called by helper routines
-that know details of the particular binary arithmetic of interest and
-convert. The present directory provides helper routines for 5 variants
-of IEEE binary floating-point arithmetic, each indicated by one or
-two letters:
-
- f IEEE single precision
- d IEEE double precision
- x IEEE extended precision, as on Intel 80x87
- and software emulations of Motorola 68xxx chips
- that do not pad the way the 68xxx does, but
- only store 80 bits
- xL IEEE extended precision, as on Motorola 68xxx chips
- Q quad precision, as on Sun Sparc chips
- dd double double, pairs of IEEE double numbers
- whose sum is the desired value
-
-For decimal -> binary conversions, there are three families of
-helper routines: one for round-nearest (or the current rounding
-mode on IEEE-arithmetic systems that provide the C99 fegetround()
-function, if compiled with -DHonor_FLT_ROUNDS):
-
- strtof
- strtod
- strtodd
- strtopd
- strtopf
- strtopx
- strtopxL
- strtopQ
-
-one with rounding direction specified:
-
- strtorf
- strtord
- strtordd
- strtorx
- strtorxL
- strtorQ
-
-and one for computing an interval (at most one bit wide) that contains
-the decimal number:
-
- strtoIf
- strtoId
- strtoIdd
- strtoIx
- strtoIxL
- strtoIQ
-
-The latter call strtoIg, which makes one call on strtodg and adjusts
-the result to provide the desired interval. On systems where native
-arithmetic can easily make one-ulp adjustments on values in the
-desired floating-point format, it might be more efficient to use the
-native arithmetic. Routine strtodI is a variant of strtoId that
-illustrates one way to do this for IEEE binary double-precision
-arithmetic -- but whether this is more efficient remains to be seen.
-
-Functions strtod and strtof have "natural" return types, float and
-double -- strtod is specified by the C standard, and strtof appears
-in the stdlib.h of some systems, such as (at least some) Linux systems.
-The other functions write their results to their final argument(s):
-to the final two argument for the strtoI... (interval) functions,
-and to the final argument for the others (strtop... and strtor...).
-Where possible, these arguments have "natural" return types (double*
-or float*), to permit at least some type checking. In reality, they
-are viewed as arrays of ULong (or, for the "x" functions, UShort)
-values. On systems where long double is the appropriate type, one can
-pass long double* final argument(s) to these routines. The int value
-that these routines return is the return value from the call they make
-on strtodg; see the enum of possible return values in gdtoa.h.
-
-Source files g_ddfmt.c, misc.c, smisc.c, strtod.c, strtodg.c, and ulp.c
-should use true IEEE double arithmetic (not, e.g., double extended),
-at least for storing (and viewing the bits of) the variables declared
-"double" within them.
-
-One detail indicated in struct FPI is whether the target binary
-arithmetic departs from the IEEE standard by flushing denormalized
-numbers to 0. On systems that do this, the helper routines for
-conversion to double-double format (when compiled with
-Sudden_Underflow #defined) penalize the bottom of the exponent
-range so that they return a nonzero result only when the least
-significant bit of the less significant member of the pair of
-double values returned can be expressed as a normalized double
-value. An alternative would be to drop to 53-bit precision near
-the bottom of the exponent range. To get correct rounding, this
-would (in general) require two calls on strtodg (one specifying
-126-bit arithmetic, then, if necessary, one specifying 53-bit
-arithmetic).
-
-By default, the core routine strtodg and strtod set errno to ERANGE
-if the result overflows to +Infinity or underflows to 0. Compile
-these routines with NO_ERRNO #defined to inhibit errno assignments.
-
-Routine strtod is based on netlib's "dtoa.c from fp", and
-(f = strtod(s,se)) is more efficient for some conversions than, say,
-strtord(s,se,1,&f). Parts of strtod require true IEEE double
-arithmetic with the default rounding mode (round-to-nearest) and, on
-systems with IEEE extended-precision registers, double-precision
-(53-bit) rounding precision. If the machine uses (the equivalent of)
-Intel 80x87 arithmetic, the call
- _control87(PC_53, MCW_PC);
-does this with many compilers. Whether this or another call is
-appropriate depends on the compiler; for this to work, it may be
-necessary to #include "float.h" or another system-dependent header
-file.
-
-Source file strtodnrp.c gives a strtod that does not require 53-bit
-rounding precision on systems (such as Intel IA32 systems) that may
-suffer double rounding due to use of extended-precision registers.
-For some conversions this variant of strtod is less efficient than the
-one in strtod.c when the latter is run with 53-bit rounding precision.
-
-The values that the strto* routines return for NaNs are determined by
-gd_qnan.h, which the makefile generates by running the program whose
-source is qnan.c. Note that the rules for distinguishing signaling
-from quiet NaNs are system-dependent. For cross-compilation, you need
-to determine arith.h and gd_qnan.h suitably, e.g., using the
-arithmetic of the target machine.
-
-C99's hexadecimal floating-point constants are recognized by the
-strto* routines (but this feature has not yet been heavily tested).
-Compiling with NO_HEX_FP #defined disables this feature.
-
-When compiled with -DINFNAN_CHECK, the strto* routines recognize C99's
-NaN and Infinity syntax. Moreover, unless No_Hex_NaN is #defined, the
-strto* routines also recognize C99's NaN(...) syntax: they accept
-(case insensitively) strings of the form NaN(x), where x is a string
-of hexadecimal digits and spaces; if there is only one string of
-hexadecimal digits, it is taken for the fraction bits of the resulting
-NaN; if there are two or more strings of hexadecimal digits, each
-string is assigned to the next available sequence of 32-bit words of
-fractions bits (starting with the most significant), right-aligned in
-each sequence.
-
-For binary -> decimal conversions, I've provided just one family
-of helper routines:
-
- g_ffmt
- g_dfmt
- g_ddfmt
- g_xfmt
- g_xLfmt
- g_Qfmt
-
-which do a "%g" style conversion either to a specified number of decimal
-places (if their ndig argument is positive), or to the shortest
-decimal string that rounds to the given binary floating-point value
-(if ndig <= 0). They write into a buffer supplied as an argument
-and return either a pointer to the end of the string (a null character)
-in the buffer, if the buffer was long enough, or 0. Other forms of
-conversion are easily done with the help of gdtoa(), such as %e or %f
-style and conversions with direction of rounding specified (so that, if
-desired, the decimal value is either >= or <= the binary value).
-On IEEE-arithmetic systems that provide the C99 fegetround() function,
-if compiled with -DHonor_FLT_ROUNDS, these routines honor the current
-rounding mode.
-
-For an example of more general conversions based on dtoa(), see
-netlib's "printf.c from ampl/solvers".
-
-For double-double -> decimal, g_ddfmt() assumes IEEE-like arithmetic
-of precision max(126, #bits(input)) bits, where #bits(input) is the
-number of mantissa bits needed to represent the sum of the two double
-values in the input.
-
-The makefile creates a library, gdtoa.a. To use the helper
-routines, a program only needs to include gdtoa.h. All the
-source files for gdtoa.a include a more extensive gdtoaimp.h;
-among other things, gdtoaimp.h has #defines that make "internal"
-names end in _D2A. To make a "system" library, one could modify
-these #defines to make the names start with __.
-
-Various comments about possible #defines appear in gdtoaimp.h,
-but for most purposes, arith.h should set suitable #defines.
-
-Systems with preemptive scheduling of multiple threads require some
-manual intervention. On such systems, it's necessary to compile
-dmisc.c, dtoa.c gdota.c, and misc.c with MULTIPLE_THREADS #defined,
-and to provide (or suitably #define) two locks, acquired by
-ACQUIRE_DTOA_LOCK(n) and freed by FREE_DTOA_LOCK(n) for n = 0 or 1.
-(The second lock, accessed in pow5mult, ensures lazy evaluation of
-only one copy of high powers of 5; omitting this lock would introduce
-a small probability of wasting memory, but would otherwise be harmless.)
-Routines that call dtoa or gdtoa directly must also invoke freedtoa(s)
-to free the value s returned by dtoa or gdtoa. It's OK to do so whether
-or not MULTIPLE_THREADS is #defined, and the helper g_*fmt routines
-listed above all do this indirectly (in gfmt_D2A(), which they all call).
-
-By default, there is a private pool of memory of length 2000 bytes
-for intermediate quantities, and MALLOC (see gdtoaimp.h) is called only
-if the private pool does not suffice. 2000 is large enough that MALLOC
-is called only under very unusual circumstances (decimal -> binary
-conversion of very long strings) for conversions to and from double
-precision. For systems with preemptively scheduled multiple threads
-or for conversions to extended or quad, it may be appropriate to
-#define PRIVATE_MEM nnnn, where nnnn is a suitable value > 2000.
-For extended and quad precisions, -DPRIVATE_MEM=20000 is probably
-plenty even for many digits at the ends of the exponent range.
-Use of the private pool avoids some overhead.
-
-Directory test provides some test routines. See its README.
-I've also tested this stuff (except double double conversions)
-with Vern Paxson's testbase program: see
-
- V. Paxson and W. Kahan, "A Program for Testing IEEE Binary-Decimal
- Conversion", manuscript, May 1991,
- ftp://ftp.ee.lbl.gov/testbase-report.ps.Z .
-
-(The same ftp directory has source for testbase.)
-
-Some system-dependent additions to CFLAGS in the makefile:
-
- HU-UX: -Aa -Ae
- OSF (DEC Unix): -ieee_with_no_inexact
- SunOS 4.1x: -DKR_headers -DBad_float_h
-
-If you want to put this stuff into a shared library and your
-operating system requires export lists for shared libraries,
-the following would be an appropriate export list:
-
- dtoa
- freedtoa
- g_Qfmt
- g_ddfmt
- g_dfmt
- g_ffmt
- g_xLfmt
- g_xfmt
- gdtoa
- strtoIQ
- strtoId
- strtoIdd
- strtoIf
- strtoIx
- strtoIxL
- strtod
- strtodI
- strtodg
- strtof
- strtopQ
- strtopd
- strtopdd
- strtopf
- strtopx
- strtopxL
- strtorQ
- strtord
- strtordd
- strtorf
- strtorx
- strtorxL
-
-When time permits, I (dmg) hope to write in more detail about the
-present conversion routines; for now, this README file must suffice.
-Meanwhile, if you wish to write helper functions for other kinds of
-IEEE-like arithmetic, some explanation of struct FPI and the bits
-array may be helpful. Both gdtoa and strtodg operate on a bits array
-described by FPI *fpi. The bits array is of type ULong, a 32-bit
-unsigned integer type. Floating-point numbers have fpi->nbits bits,
-with the least significant 32 bits in bits[0], the next 32 bits in
-bits[1], etc. These numbers are regarded as integers multiplied by
-2^e (i.e., 2 to the power of the exponent e), where e is the second
-argument (be) to gdtoa and is stored in *exp by strtodg. The minimum
-and maximum exponent values fpi->emin and fpi->emax for normalized
-floating-point numbers reflect this arrangement. For example, the
-P754 standard for binary IEEE arithmetic specifies doubles as having
-53 bits, with normalized values of the form 1.xxxxx... times 2^(b-1023),
-with 52 bits (the x's) and the biased exponent b represented explicitly;
-b is an unsigned integer in the range 1 <= b <= 2046 for normalized
-finite doubles, b = 0 for denormals, and b = 2047 for Infinities and NaNs.
-To turn an IEEE double into the representation used by strtodg and gdtoa,
-we multiply 1.xxxx... by 2^52 (to make it an integer) and reduce the
-exponent e = (b-1023) by 52:
-
- fpi->emin = 1 - 1023 - 52
- fpi->emax = 1046 - 1023 - 52
-
-In various wrappers for IEEE double, we actually write -53 + 1 rather
-than -52, to emphasize that there are 53 bits including one implicit bit.
-Field fpi->rounding indicates the desired rounding direction, with
-possible values
- FPI_Round_zero = toward 0,
- FPI_Round_near = unbiased rounding -- the IEEE default,
- FPI_Round_up = toward +Infinity, and
- FPI_Round_down = toward -Infinity
-given in gdtoa.h.
-
-Field fpi->sudden_underflow indicates whether strtodg should return
-denormals or flush them to zero. Normal floating-point numbers have
-bit fpi->nbits in the bits array on. Denormals have it off, with
-exponent = fpi->emin. Strtodg provides distinct return values for normals
-and denormals; see gdtoa.h.
-
-Compiling g__fmt.c, strtod.c, and strtodg.c with -DUSE_LOCALE causes
-the decimal-point character to be taken from the current locale; otherwise
-it is '.'.
-
-Source files dtoa.c and strtod.c in this directory are derived from
-netlib's "dtoa.c from fp" and are meant to function equivalently.
-When compiled with Honor_FLT_ROUNDS #defined (on systems that provide
-FLT_ROUNDS and fegetround() as specified in the C99 standard), they
-honor the current rounding mode. Because FLT_ROUNDS is buggy on some
-(Linux) systems -- not reflecting calls on fesetround(), as the C99
-standard says it should -- when Honor_FLT_ROUNDS is #defined, the
-current rounding mode is obtained from fegetround() rather than from
-FLT_ROUNDS, unless Trust_FLT_ROUNDS is also #defined.
-
-Compile with -DUSE_LOCALE to use the current locale; otherwise
-decimal points are assumed to be '.'. With -DUSE_LOCALE, unless
-you also compile with -DNO_LOCALE_CACHE, the details about the
-current "decimal point" character string are cached and assumed not
-to change during the program's execution.
-
-On machines with a 64-bit long double and perhaps a 113-bit "quad"
-type, you can invoke "make Printf" to add Printf (and variants, such
-as Fprintf) to gdtoa.a. These are analogs, declared in stdio1.h, of
-printf and fprintf, etc. in which %La, %Le, %Lf, and %Lg are for long
-double and (if appropriate) %Lqa, %Lqe, %Lqf, and %Lqg are for quad
-precision printing.
-
-Please send comments to David M. Gay (dmg at acm dot org, with " at "
-changed at "@" and " dot " changed to ".").
+++ /dev/null
-#include <stdint.h>
-
-#if __LITTLE_ENDIAN__
-# define IEEE_8087
-#elif __BIG_ENDIAN__
-# define IEEE_MC68k
-#else
-// Because quad math is typically performed on little-endian hardware
-//# error "Unknown endianness"
-# define IEEE_8087
-#endif
-
-// This should be 32 bit integer type
-#define Long int
+++ /dev/null
-#define f_QNAN 0x7fc00000
-#define d_QNAN0 0x0
-#define d_QNAN1 0x7ff80000
-#define ld_QNAN0 0x0
-#define ld_QNAN1 0x0
-#define ld_QNAN2 0x0
-#define ld_QNAN3 0x7fff8000
-#define ldus_QNAN0 0x0
-#define ldus_QNAN1 0x0
-#define ldus_QNAN2 0x0
-#define ldus_QNAN3 0xc000
-#define ldus_QNAN4 0x7fff
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#ifndef GDTOA_H_INCLUDED
-#define GDTOA_H_INCLUDED
-
-#include "arith.h"
-#include <stddef.h> /* for size_t */
-#include <stdlib.h> /* for strtod */
-
-#ifndef Long
-#define Long long
-#endif
-#ifndef ULong
-typedef unsigned Long ULong;
-#endif
-#ifndef UShort
-typedef unsigned short UShort;
-#endif
-
-#ifndef ANSI
-#ifdef KR_headers
-#define ANSI(x) ()
-#define Void /*nothing*/
-#else
-#define ANSI(x) x
-#define Void void
-#endif
-#endif /* ANSI */
-
-#ifndef CONST
-#ifdef KR_headers
-#define CONST /* blank */
-#else
-#define CONST const
-#endif
-#endif /* CONST */
-
- enum { /* return values from strtodg */
- STRTOG_Zero = 0,
- STRTOG_Normal = 1,
- STRTOG_Denormal = 2,
- STRTOG_Infinite = 3,
- STRTOG_NaN = 4,
- STRTOG_NaNbits = 5,
- STRTOG_NoNumber = 6,
- STRTOG_Retmask = 7,
-
- /* The following may be or-ed into one of the above values. */
-
- STRTOG_Neg = 0x08, /* does not affect STRTOG_Inexlo or STRTOG_Inexhi */
- STRTOG_Inexlo = 0x10, /* returned result rounded toward zero */
- STRTOG_Inexhi = 0x20, /* returned result rounded away from zero */
- STRTOG_Inexact = 0x30,
- STRTOG_Underflow= 0x40,
- STRTOG_Overflow = 0x80
- };
-
- typedef struct
-FPI {
- int nbits;
- int emin;
- int emax;
- int rounding;
- int sudden_underflow;
- } FPI;
-
-enum { /* FPI.rounding values: same as FLT_ROUNDS */
- FPI_Round_zero = 0,
- FPI_Round_near = 1,
- FPI_Round_up = 2,
- FPI_Round_down = 3
- };
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-extern char* dtoa ANSI((double d, int mode, int ndigits, int *decpt,
- int *sign, char **rve));
-extern char* gdtoa ANSI((FPI *fpi, int be, ULong *bits, int *kindp,
- int mode, int ndigits, int *decpt, char **rve));
-extern void freedtoa ANSI((char*));
-extern int strtodg ANSI((CONST char*, char**, FPI*, Long*, ULong*));
-
-extern char* g_ddfmt ANSI((char*, double*, int, size_t));
-extern char* g_dfmt ANSI((char*, double*, int, size_t));
-extern char* g_ffmt ANSI((char*, float*, int, size_t));
-extern char* g_Qfmt ANSI((char*, void*, int, size_t));
-extern char* g_xfmt ANSI((char*, void*, int, size_t));
-extern char* g_xLfmt ANSI((char*, void*, int, size_t));
-
-extern int strtoId ANSI((CONST char*, char**, double*, double*));
-extern int strtoIdd ANSI((CONST char*, char**, double*, double*));
-extern int strtoIf ANSI((CONST char*, char**, float*, float*));
-extern int strtoIQ ANSI((CONST char*, char**, void*, void*));
-extern int strtoIx ANSI((CONST char*, char**, void*, void*));
-extern int strtoIxL ANSI((CONST char*, char**, void*, void*));
-extern int strtord ANSI((CONST char*, char**, int, double*));
-extern int strtordd ANSI((CONST char*, char**, int, double*));
-extern int strtorf ANSI((CONST char*, char**, int, float*));
-extern int strtorQ ANSI((CONST char*, char**, int, void*));
-extern int strtorx ANSI((CONST char*, char**, int, void*));
-extern int strtorxL ANSI((CONST char*, char**, int, void*));
-#if 1
-extern int strtodI ANSI((CONST char*, char**, double*));
-extern int strtopd ANSI((CONST char*, char**, double*));
-extern int strtopdd ANSI((CONST char*, char**, double*));
-extern int strtopf ANSI((CONST char*, char**, float*));
-extern int quadmath_strtopQ ANSI((CONST char*, char**, void*));
-extern int strtopx ANSI((CONST char*, char**, void*));
-extern int strtopxL ANSI((CONST char*, char**, void*));
-#else
-#define strtopd(s,se,x) strtord(s,se,1,x)
-#define strtopdd(s,se,x) strtordd(s,se,1,x)
-#define strtopf(s,se,x) strtorf(s,se,1,x)
-#define strtopQ(s,se,x) strtorQ(s,se,1,x)
-#define strtopx(s,se,x) strtorx(s,se,1,x)
-#define strtopxL(s,se,x) strtorxL(s,se,1,x)
-#endif
-
-#ifdef __cplusplus
-}
-#endif
-#endif /* GDTOA_H_INCLUDED */
+++ /dev/null
- FPI *fpi, fpi1;
- int Rounding;
-#ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
- Rounding = Flt_Rounds;
-#else /*}{*/
- Rounding = 1;
- switch(fegetround()) {
- case FE_TOWARDZERO: Rounding = 0; break;
- case FE_UPWARD: Rounding = 2; break;
- case FE_DOWNWARD: Rounding = 3;
- }
-#endif /*}}*/
- fpi = &fpi0;
- if (Rounding != 1) {
- fpi1 = fpi0;
- fpi = &fpi1;
- fpi1.rounding = Rounding;
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998-2000 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* This is a variation on dtoa.c that converts arbitary binary
- floating-point formats to and from decimal notation. It uses
- double-precision arithmetic internally, so there are still
- various #ifdefs that adapt the calculations to the native
- double-precision arithmetic (any of IEEE, VAX D_floating,
- or IBM mainframe arithmetic).
-
- Please send bug reports to David M. Gay (dmg at acm dot org,
- with " at " changed at "@" and " dot " changed to ".").
- */
-
-/* On a machine with IEEE extended-precision registers, it is
- * necessary to specify double-precision (53-bit) rounding precision
- * before invoking strtod or dtoa. If the machine uses (the equivalent
- * of) Intel 80x87 arithmetic, the call
- * _control87(PC_53, MCW_PC);
- * does this with many compilers. Whether this or another call is
- * appropriate depends on the compiler; for this to work, it may be
- * necessary to #include "float.h" or another system-dependent header
- * file.
- */
-
-/* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
- *
- * This strtod returns a nearest machine number to the input decimal
- * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
- * broken by the IEEE round-even rule. Otherwise ties are broken by
- * biased rounding (add half and chop).
- *
- * Inspired loosely by William D. Clinger's paper "How to Read Floating
- * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 112-126].
- *
- * Modifications:
- *
- * 1. We only require IEEE, IBM, or VAX double-precision
- * arithmetic (not IEEE double-extended).
- * 2. We get by with floating-point arithmetic in a case that
- * Clinger missed -- when we're computing d * 10^n
- * for a small integer d and the integer n is not too
- * much larger than 22 (the maximum integer k for which
- * we can represent 10^k exactly), we may be able to
- * compute (d*10^k) * 10^(e-k) with just one roundoff.
- * 3. Rather than a bit-at-a-time adjustment of the binary
- * result in the hard case, we use floating-point
- * arithmetic to determine the adjustment to within
- * one bit; only in really hard cases do we need to
- * compute a second residual.
- * 4. Because of 3., we don't need a large table of powers of 10
- * for ten-to-e (just some small tables, e.g. of 10^k
- * for 0 <= k <= 22).
- */
-
-/*
- * #define IEEE_8087 for IEEE-arithmetic machines where the least
- * significant byte has the lowest address.
- * #define IEEE_MC68k for IEEE-arithmetic machines where the most
- * significant byte has the lowest address.
- * #define Long int on machines with 32-bit ints and 64-bit longs.
- * #define Sudden_Underflow for IEEE-format machines without gradual
- * underflow (i.e., that flush to zero on underflow).
- * #define IBM for IBM mainframe-style floating-point arithmetic.
- * #define VAX for VAX-style floating-point arithmetic (D_floating).
- * #define No_leftright to omit left-right logic in fast floating-point
- * computation of dtoa.
- * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3.
- * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
- * that use extended-precision instructions to compute rounded
- * products and quotients) with IBM.
- * #define ROUND_BIASED for IEEE-format with biased rounding.
- * #define Inaccurate_Divide for IEEE-format with correctly rounded
- * products but inaccurate quotients, e.g., for Intel i860.
- * #define NO_LONG_LONG on machines that do not have a "long long"
- * integer type (of >= 64 bits). On such machines, you can
- * #define Just_16 to store 16 bits per 32-bit Long when doing
- * high-precision integer arithmetic. Whether this speeds things
- * up or slows things down depends on the machine and the number
- * being converted. If long long is available and the name is
- * something other than "long long", #define Llong to be the name,
- * and if "unsigned Llong" does not work as an unsigned version of
- * Llong, #define #ULLong to be the corresponding unsigned type.
- * #define KR_headers for old-style C function headers.
- * #define Bad_float_h if your system lacks a float.h or if it does not
- * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
- * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
- * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
- * if memory is available and otherwise does something you deem
- * appropriate. If MALLOC is undefined, malloc will be invoked
- * directly -- and assumed always to succeed. Similarly, if you
- * want something other than the system's free() to be called to
- * recycle memory acquired from MALLOC, #define FREE to be the
- * name of the alternate routine. (FREE or free is only called in
- * pathological cases, e.g., in a gdtoa call after a gdtoa return in
- * mode 3 with thousands of digits requested.)
- * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
- * memory allocations from a private pool of memory when possible.
- * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
- * unless #defined to be a different length. This default length
- * suffices to get rid of MALLOC calls except for unusual cases,
- * such as decimal-to-binary conversion of a very long string of
- * digits. When converting IEEE double precision values, the
- * longest string gdtoa can return is about 751 bytes long. For
- * conversions by strtod of strings of 800 digits and all gdtoa
- * conversions of IEEE doubles in single-threaded executions with
- * 8-byte pointers, PRIVATE_MEM >= 7400 appears to suffice; with
- * 4-byte pointers, PRIVATE_MEM >= 7112 appears adequate.
- * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
- * #defined automatically on IEEE systems. On such systems,
- * when INFNAN_CHECK is #defined, strtod checks
- * for Infinity and NaN (case insensitively).
- * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
- * strtodg also accepts (case insensitively) strings of the form
- * NaN(x), where x is a string of hexadecimal digits (optionally
- * preceded by 0x or 0X) and spaces; if there is only one string
- * of hexadecimal digits, it is taken for the fraction bits of the
- * resulting NaN; if there are two or more strings of hexadecimal
- * digits, each string is assigned to the next available sequence
- * of 32-bit words of fractions bits (starting with the most
- * significant), right-aligned in each sequence.
- * Unless GDTOA_NON_PEDANTIC_NANCHECK is #defined, input "NaN(...)"
- * is consumed even when ... has the wrong form (in which case the
- * "(...)" is consumed but ignored).
- * #define MULTIPLE_THREADS if the system offers preemptively scheduled
- * multiple threads. In this case, you must provide (or suitably
- * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
- * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
- * in pow5mult, ensures lazy evaluation of only one copy of high
- * powers of 5; omitting this lock would introduce a small
- * probability of wasting memory, but would otherwise be harmless.)
- * You must also invoke freedtoa(s) to free the value s returned by
- * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
- * #define IMPRECISE_INEXACT if you do not care about the setting of
- * the STRTOG_Inexact bits in the special case of doing IEEE double
- * precision conversions (which could also be done by the strtod in
- * dtoa.c).
- * #define NO_HEX_FP to disable recognition of C9x's hexadecimal
- * floating-point constants.
- * #define -DNO_ERRNO to suppress setting errno (in strtod.c and
- * strtodg.c).
- * #define NO_STRING_H to use private versions of memcpy.
- * On some K&R systems, it may also be necessary to
- * #define DECLARE_SIZE_T in this case.
- * #define USE_LOCALE to use the current locale's decimal_point value.
- */
-
-#ifndef GDTOAIMP_H_INCLUDED
-#define GDTOAIMP_H_INCLUDED
-#include "gdtoa.h"
-#include "gd_qnan.h"
-#ifdef Honor_FLT_ROUNDS
-#include <fenv.h>
-#endif
-
-#ifdef DEBUG
-#include "stdio.h"
-#define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
-#endif
-
-#include "stdlib.h"
-#include "string.h"
-
-#ifdef KR_headers
-#define Char char
-#else
-#define Char void
-#endif
-
-#ifdef MALLOC
-extern Char *MALLOC ANSI((size_t));
-#else
-#define MALLOC malloc
-#endif
-
-#undef IEEE_Arith
-#undef Avoid_Underflow
-#ifdef IEEE_MC68k
-#define IEEE_Arith
-#endif
-#ifdef IEEE_8087
-#define IEEE_Arith
-#endif
-
-#include "errno.h"
-#ifdef Bad_float_h
-
-#ifdef IEEE_Arith
-#define DBL_DIG 15
-#define DBL_MAX_10_EXP 308
-#define DBL_MAX_EXP 1024
-#define FLT_RADIX 2
-#define DBL_MAX 1.7976931348623157e+308
-#endif
-
-#ifdef IBM
-#define DBL_DIG 16
-#define DBL_MAX_10_EXP 75
-#define DBL_MAX_EXP 63
-#define FLT_RADIX 16
-#define DBL_MAX 7.2370055773322621e+75
-#endif
-
-#ifdef VAX
-#define DBL_DIG 16
-#define DBL_MAX_10_EXP 38
-#define DBL_MAX_EXP 127
-#define FLT_RADIX 2
-#define DBL_MAX 1.7014118346046923e+38
-#define n_bigtens 2
-#endif
-
-#ifndef LONG_MAX
-#define LONG_MAX 2147483647
-#endif
-
-#else /* ifndef Bad_float_h */
-#include "float.h"
-#endif /* Bad_float_h */
-
-#ifdef IEEE_Arith
-#define Scale_Bit 0x10
-#define n_bigtens 5
-#endif
-
-#ifdef IBM
-#define n_bigtens 3
-#endif
-
-#ifdef VAX
-#define n_bigtens 2
-#endif
-
-#ifndef __MATH_H__
-#include "math.h"
-#endif
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
-Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
-#endif
-
-typedef union { double d; ULong L[2]; } U;
-
-#ifdef IEEE_8087
-#define word0(x) (x)->L[1]
-#define word1(x) (x)->L[0]
-#else
-#define word0(x) (x)->L[0]
-#define word1(x) (x)->L[1]
-#endif
-#define dval(x) (x)->d
-
-/* The following definition of Storeinc is appropriate for MIPS processors.
- * An alternative that might be better on some machines is
- * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
- */
-#if defined(IEEE_8087) + defined(VAX)
-#define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
-((unsigned short *)a)[0] = (unsigned short)c, a++)
-#else
-#define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
-((unsigned short *)a)[1] = (unsigned short)c, a++)
-#endif
-
-/* #define P DBL_MANT_DIG */
-/* Ten_pmax = floor(P*log(2)/log(5)) */
-/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
-/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
-/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
-
-#ifdef IEEE_Arith
-#define Exp_shift 20
-#define Exp_shift1 20
-#define Exp_msk1 0x100000
-#define Exp_msk11 0x100000
-#define Exp_mask 0x7ff00000
-#define P 53
-#define Bias 1023
-#define Emin (-1022)
-#define Exp_1 0x3ff00000
-#define Exp_11 0x3ff00000
-#define Ebits 11
-#define Frac_mask 0xfffff
-#define Frac_mask1 0xfffff
-#define Ten_pmax 22
-#define Bletch 0x10
-#define Bndry_mask 0xfffff
-#define Bndry_mask1 0xfffff
-#define LSB 1
-#define Sign_bit 0x80000000
-#define Log2P 1
-#define Tiny0 0
-#define Tiny1 1
-#define Quick_max 14
-#define Int_max 14
-
-#ifndef Flt_Rounds
-#ifdef FLT_ROUNDS
-#define Flt_Rounds FLT_ROUNDS
-#else
-#define Flt_Rounds 1
-#endif
-#endif /*Flt_Rounds*/
-
-#else /* ifndef IEEE_Arith */
-#undef Sudden_Underflow
-#define Sudden_Underflow
-#ifdef IBM
-#undef Flt_Rounds
-#define Flt_Rounds 0
-#define Exp_shift 24
-#define Exp_shift1 24
-#define Exp_msk1 0x1000000
-#define Exp_msk11 0x1000000
-#define Exp_mask 0x7f000000
-#define P 14
-#define Bias 65
-#define Exp_1 0x41000000
-#define Exp_11 0x41000000
-#define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
-#define Frac_mask 0xffffff
-#define Frac_mask1 0xffffff
-#define Bletch 4
-#define Ten_pmax 22
-#define Bndry_mask 0xefffff
-#define Bndry_mask1 0xffffff
-#define LSB 1
-#define Sign_bit 0x80000000
-#define Log2P 4
-#define Tiny0 0x100000
-#define Tiny1 0
-#define Quick_max 14
-#define Int_max 15
-#else /* VAX */
-#undef Flt_Rounds
-#define Flt_Rounds 1
-#define Exp_shift 23
-#define Exp_shift1 7
-#define Exp_msk1 0x80
-#define Exp_msk11 0x800000
-#define Exp_mask 0x7f80
-#define P 56
-#define Bias 129
-#define Exp_1 0x40800000
-#define Exp_11 0x4080
-#define Ebits 8
-#define Frac_mask 0x7fffff
-#define Frac_mask1 0xffff007f
-#define Ten_pmax 24
-#define Bletch 2
-#define Bndry_mask 0xffff007f
-#define Bndry_mask1 0xffff007f
-#define LSB 0x10000
-#define Sign_bit 0x8000
-#define Log2P 1
-#define Tiny0 0x80
-#define Tiny1 0
-#define Quick_max 15
-#define Int_max 15
-#endif /* IBM, VAX */
-#endif /* IEEE_Arith */
-
-#ifndef IEEE_Arith
-#define ROUND_BIASED
-#endif
-
-#ifdef RND_PRODQUOT
-#define rounded_product(a,b) a = rnd_prod(a, b)
-#define rounded_quotient(a,b) a = rnd_quot(a, b)
-#ifdef KR_headers
-extern double rnd_prod(), rnd_quot();
-#else
-extern double rnd_prod(double, double), rnd_quot(double, double);
-#endif
-#else
-#define rounded_product(a,b) a *= b
-#define rounded_quotient(a,b) a /= b
-#endif
-
-#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
-#define Big1 0xffffffff
-
-#undef Pack_16
-#ifndef Pack_32
-#define Pack_32
-#endif
-
-#ifdef NO_LONG_LONG
-#undef ULLong
-#ifdef Just_16
-#undef Pack_32
-#define Pack_16
-/* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
- * This makes some inner loops simpler and sometimes saves work
- * during multiplications, but it often seems to make things slightly
- * slower. Hence the default is now to store 32 bits per Long.
- */
-#endif
-#else /* long long available */
-#ifndef Llong
-#define Llong long long
-#endif
-#ifndef ULLong
-#define ULLong unsigned Llong
-#endif
-#endif /* NO_LONG_LONG */
-
-#ifdef Pack_32
-#define ULbits 32
-#define kshift 5
-#define kmask 31
-#define ALL_ON 0xffffffff
-#else
-#define ULbits 16
-#define kshift 4
-#define kmask 15
-#define ALL_ON 0xffff
-#endif
-
-#ifndef MULTIPLE_THREADS
-#define ACQUIRE_DTOA_LOCK(n) /*nothing*/
-#define FREE_DTOA_LOCK(n) /*nothing*/
-#endif
-
-#define Kmax 9
-
- struct
-Bigint {
- struct Bigint *next;
- int k, maxwds, sign, wds;
- ULong x[1];
- };
-
- typedef struct Bigint Bigint;
-
-#ifdef NO_STRING_H
-#ifdef DECLARE_SIZE_T
-typedef unsigned int size_t;
-#endif
-extern void memcpy_D2A ANSI((void*, const void*, size_t));
-#define Bcopy(x,y) memcpy_D2A(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
-#else /* !NO_STRING_H */
-#define Bcopy(x,y) memcpy(&x->sign,&y->sign,y->wds*sizeof(ULong) + 2*sizeof(int))
-#endif /* NO_STRING_H */
-
-#define Balloc __quadmath_Balloc_D2A
-#define Bfree __quadmath_Bfree_D2A
-#define ULtoQ __quadmath_ULtoQ_D2A
-#define ULtof __quadmath_ULtof_D2A
-#define ULtod __quadmath_ULtod_D2A
-#define ULtodd __quadmath_ULtodd_D2A
-#define ULtox __quadmath_ULtox_D2A
-#define ULtoxL __quadmath_ULtoxL_D2A
-#define any_on __quadmath_any_on_D2A
-#define b2d __quadmath_b2d_D2A
-#define bigtens __quadmath_bigtens_D2A
-#define cmp __quadmath_cmp_D2A
-#define copybits __quadmath_copybits_D2A
-#define d2b __quadmath_d2b_D2A
-#define decrement __quadmath_decrement_D2A
-#define diff __quadmath_diff_D2A
-#define dtoa_result __quadmath_dtoa_result_D2A
-#define g__fmt __quadmath_g__fmt_D2A
-#define gethex __quadmath_gethex_D2A
-#define hexdig __quadmath_hexdig_D2A
-#define hexnan __quadmath_hexnan_D2A
-#define hi0bits_D2A __quadmath_hi0bits_D2A
-#define hi0bits(x) __quadmath_hi0bits_D2A((ULong)(x))
-#define i2b __quadmath_i2b_D2A
-#define increment __quadmath_increment_D2A
-#define lo0bits __quadmath_lo0bits_D2A
-#define lshift __quadmath_lshift_D2A
-#define match __quadmath_match_D2A
-#define mult __quadmath_mult_D2A
-#define multadd __quadmath_multadd_D2A
-#define nrv_alloc __quadmath_nrv_alloc_D2A
-#define pow5mult __quadmath_pow5mult_D2A
-#define quorem __quadmath_quorem_D2A
-#define ratio __quadmath_ratio_D2A
-#define rshift __quadmath_rshift_D2A
-#define rv_alloc __quadmath_rv_alloc_D2A
-#define s2b __quadmath_s2b_D2A
-#define set_ones __quadmath_set_ones_D2A
-#ifdef NO_STRING_H
-#define strcp __quadmath_strcp_D2A
-#endif
-#define strtoIg __quadmath_strtoIg_D2A
-#define sum __quadmath_sum_D2A
-#define tens __quadmath_tens_D2A
-#define tinytens __quadmath_tinytens_D2A
-#define tinytens __quadmath_tinytens_D2A
-#define trailz __quadmath_trailz_D2A
-#define ulp __quadmath_ulp_D2A
-
- extern char *dtoa_result;
- extern CONST double bigtens[], tens[], tinytens[];
- extern unsigned char hexdig[];
-
- extern Bigint *Balloc ANSI((int));
- extern void Bfree ANSI((Bigint*));
- extern void ULtof ANSI((ULong*, ULong*, Long, int));
- extern void ULtod ANSI((ULong*, ULong*, Long, int));
- extern void ULtodd ANSI((ULong*, ULong*, Long, int));
- extern void ULtoQ ANSI((ULong*, ULong*, Long, int));
- extern void ULtox ANSI((UShort*, ULong*, Long, int));
- extern void ULtoxL ANSI((ULong*, ULong*, Long, int));
- extern ULong any_on ANSI((Bigint*, int));
- extern double b2d ANSI((Bigint*, int*));
- extern int cmp ANSI((Bigint*, Bigint*));
- extern void copybits ANSI((ULong*, int, Bigint*));
- extern Bigint *d2b ANSI((double, int*, int*));
- extern void decrement ANSI((Bigint*));
- extern Bigint *diff ANSI((Bigint*, Bigint*));
- extern char *dtoa ANSI((double d, int mode, int ndigits,
- int *decpt, int *sign, char **rve));
- extern char *g__fmt ANSI((char*, char*, char*, int, ULong, size_t));
- extern int gethex ANSI((CONST char**, FPI*, Long*, Bigint**, int));
- extern void hexdig_init_D2A(Void);
- extern int hexnan ANSI((CONST char**, FPI*, ULong*));
- extern int hi0bits_D2A ANSI((ULong));
- extern Bigint *i2b ANSI((int));
- extern Bigint *increment ANSI((Bigint*));
- extern int lo0bits ANSI((ULong*));
- extern Bigint *lshift ANSI((Bigint*, int));
- extern int match ANSI((CONST char**, char*));
- extern Bigint *mult ANSI((Bigint*, Bigint*));
- extern Bigint *multadd ANSI((Bigint*, int, int));
- extern char *nrv_alloc ANSI((char*, char **, int));
- extern Bigint *pow5mult ANSI((Bigint*, int));
- extern int quorem ANSI((Bigint*, Bigint*));
- extern double ratio ANSI((Bigint*, Bigint*));
- extern void rshift ANSI((Bigint*, int));
- extern char *rv_alloc ANSI((int));
- extern Bigint *s2b ANSI((CONST char*, int, int, ULong, int));
- extern Bigint *set_ones ANSI((Bigint*, int));
-#ifdef NO_STRING_H
- extern char *strcpy ANSI((char*, const char*));
-#endif
- extern int strtoIg ANSI((CONST char*, char**, FPI*, Long*, Bigint**, int*));
- extern Bigint *sum ANSI((Bigint*, Bigint*));
- extern int trailz ANSI((Bigint*));
- extern double ulp ANSI((U*));
-
-#ifdef __cplusplus
-}
-#endif
-/*
- * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c. Prior to
- * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0,
- * respectively), but now are determined by compiling and running
- * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1.
- * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=...
- * and -DNAN_WORD1=... values if necessary. This should still work.
- * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
- */
-#ifdef IEEE_Arith
-#ifndef NO_INFNAN_CHECK
-#undef INFNAN_CHECK
-#define INFNAN_CHECK
-#endif
-#ifdef IEEE_MC68k
-#define _0 0
-#define _1 1
-#ifndef NAN_WORD0
-#define NAN_WORD0 d_QNAN0
-#endif
-#ifndef NAN_WORD1
-#define NAN_WORD1 d_QNAN1
-#endif
-#else
-#define _0 1
-#define _1 0
-#ifndef NAN_WORD0
-#define NAN_WORD0 d_QNAN1
-#endif
-#ifndef NAN_WORD1
-#define NAN_WORD1 d_QNAN0
-#endif
-#endif
-#else
-#undef INFNAN_CHECK
-#endif
-
-#undef SI
-#ifdef Sudden_Underflow
-#define SI 1
-#else
-#define SI 0
-#endif
-
-#endif /* GDTOAIMP_H_INCLUDED */
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
-#ifdef USE_LOCALE
-#include "locale.h"
-#endif
-
- int
-#ifdef KR_headers
-gethex(sp, fpi, exp, bp, sign)
- CONST char **sp; FPI *fpi; Long *exp; Bigint **bp; int sign;
-#else
-gethex( CONST char **sp, FPI *fpi, Long *exp, Bigint **bp, int sign)
-#endif
-{
- Bigint *b;
- CONST unsigned char *decpt, *s0, *s, *s1;
- int big, esign, havedig, irv, j, k, n, n0, nbits, up, zret;
- ULong L, lostbits, *x;
- Long e, e1;
-#ifdef USE_LOCALE
- int i;
-#ifdef NO_LOCALE_CACHE
- const unsigned char *decimalpoint = (unsigned char*)localeconv()->decimal_point;
-#else
- const unsigned char *decimalpoint;
- static unsigned char *decimalpoint_cache;
- if (!(s0 = decimalpoint_cache)) {
- s0 = (unsigned char*)localeconv()->decimal_point;
- if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
- strcpy(decimalpoint_cache, s0);
- s0 = decimalpoint_cache;
- }
- }
- decimalpoint = s0;
-#endif
-#endif
-
- if (!hexdig['0'])
- hexdig_init_D2A();
- *bp = 0;
- havedig = 0;
- s0 = *(CONST unsigned char **)sp + 2;
- while(s0[havedig] == '0')
- havedig++;
- s0 += havedig;
- s = s0;
- decpt = 0;
- zret = 0;
- e = 0;
- if (hexdig[*s])
- havedig++;
- else {
- zret = 1;
-#ifdef USE_LOCALE
- for(i = 0; decimalpoint[i]; ++i) {
- if (s[i] != decimalpoint[i])
- goto pcheck;
- }
- decpt = s += i;
-#else
- if (*s != '.')
- goto pcheck;
- decpt = ++s;
-#endif
- if (!hexdig[*s])
- goto pcheck;
- while(*s == '0')
- s++;
- if (hexdig[*s])
- zret = 0;
- havedig = 1;
- s0 = s;
- }
- while(hexdig[*s])
- s++;
-#ifdef USE_LOCALE
- if (*s == *decimalpoint && !decpt) {
- for(i = 1; decimalpoint[i]; ++i) {
- if (s[i] != decimalpoint[i])
- goto pcheck;
- }
- decpt = s += i;
-#else
- if (*s == '.' && !decpt) {
- decpt = ++s;
-#endif
- while(hexdig[*s])
- s++;
- }/*}*/
- if (decpt)
- e = -(((Long)(s-decpt)) << 2);
- pcheck:
- s1 = s;
- big = esign = 0;
- switch(*s) {
- case 'p':
- case 'P':
- switch(*++s) {
- case '-':
- esign = 1;
- /* no break */
- case '+':
- s++;
- }
- if ((n = hexdig[*s]) == 0 || n > 0x19) {
- s = s1;
- break;
- }
- e1 = n - 0x10;
- while((n = hexdig[*++s]) !=0 && n <= 0x19) {
- if (e1 & 0xf8000000)
- big = 1;
- e1 = 10*e1 + n - 0x10;
- }
- if (esign)
- e1 = -e1;
- e += e1;
- }
- *sp = (char*)s;
- if (!havedig)
- *sp = (char*)s0 - 1;
- if (zret)
- return STRTOG_Zero;
- if (big) {
- if (esign) {
- switch(fpi->rounding) {
- case FPI_Round_up:
- if (sign)
- break;
- goto ret_tiny;
- case FPI_Round_down:
- if (!sign)
- break;
- goto ret_tiny;
- }
- goto retz;
- ret_tiny:
- b = Balloc(0);
- b->wds = 1;
- b->x[0] = 1;
- goto dret;
- }
- switch(fpi->rounding) {
- case FPI_Round_near:
- goto ovfl1;
- case FPI_Round_up:
- if (!sign)
- goto ovfl1;
- goto ret_big;
- case FPI_Round_down:
- if (sign)
- goto ovfl1;
- goto ret_big;
- }
- ret_big:
- nbits = fpi->nbits;
- n0 = n = nbits >> kshift;
- if (nbits & kmask)
- ++n;
- for(j = n, k = 0; j >>= 1; ++k);
- *bp = b = Balloc(k);
- b->wds = n;
- for(j = 0; j < n0; ++j)
- b->x[j] = ALL_ON;
- if (n > n0)
- b->x[j] = ULbits >> (ULbits - (nbits & kmask));
- *exp = fpi->emin;
- return STRTOG_Normal | STRTOG_Inexlo;
- }
- n = s1 - s0 - 1;
- for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
- k++;
- b = Balloc(k);
- x = b->x;
- n = 0;
- L = 0;
-#ifdef USE_LOCALE
- for(i = 0; decimalpoint[i+1]; ++i);
-#endif
- while(s1 > s0) {
-#ifdef USE_LOCALE
- if (*--s1 == decimalpoint[i]) {
- s1 -= i;
- continue;
- }
-#else
- if (*--s1 == '.')
- continue;
-#endif
- if (n == ULbits) {
- *x++ = L;
- L = 0;
- n = 0;
- }
- L |= (hexdig[*s1] & 0x0f) << n;
- n += 4;
- }
- *x++ = L;
- b->wds = n = x - b->x;
- n = ULbits*n - hi0bits(L);
- nbits = fpi->nbits;
- lostbits = 0;
- x = b->x;
- if (n > nbits) {
- n -= nbits;
- if (any_on(b,n)) {
- lostbits = 1;
- k = n - 1;
- if (x[k>>kshift] & 1 << (k & kmask)) {
- lostbits = 2;
- if (k > 0 && any_on(b,k))
- lostbits = 3;
- }
- }
- rshift(b, n);
- e += n;
- }
- else if (n < nbits) {
- n = nbits - n;
- b = lshift(b, n);
- e -= n;
- x = b->x;
- }
- if (e > fpi->emax) {
- ovfl:
- Bfree(b);
- ovfl1:
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- return STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
- }
- irv = STRTOG_Normal;
- if (e < fpi->emin) {
- irv = STRTOG_Denormal;
- n = fpi->emin - e;
- if (n >= nbits) {
- switch (fpi->rounding) {
- case FPI_Round_near:
- if (n == nbits && (n < 2 || any_on(b,n-1)))
- goto one_bit;
- break;
- case FPI_Round_up:
- if (!sign)
- goto one_bit;
- break;
- case FPI_Round_down:
- if (sign) {
- one_bit:
- x[0] = b->wds = 1;
- dret:
- *bp = b;
- *exp = fpi->emin;
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- return STRTOG_Denormal | STRTOG_Inexhi
- | STRTOG_Underflow;
- }
- }
- Bfree(b);
- retz:
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- return STRTOG_Zero | STRTOG_Inexlo | STRTOG_Underflow;
- }
- k = n - 1;
- if (lostbits)
- lostbits = 1;
- else if (k > 0)
- lostbits = any_on(b,k);
- if (x[k>>kshift] & 1 << (k & kmask))
- lostbits |= 2;
- nbits -= n;
- rshift(b,n);
- e = fpi->emin;
- }
- if (lostbits) {
- up = 0;
- switch(fpi->rounding) {
- case FPI_Round_zero:
- break;
- case FPI_Round_near:
- if (lostbits & 2
- && (lostbits | x[0]) & 1)
- up = 1;
- break;
- case FPI_Round_up:
- up = 1 - sign;
- break;
- case FPI_Round_down:
- up = sign;
- }
- if (up) {
- k = b->wds;
- b = increment(b);
- x = b->x;
- if (irv == STRTOG_Denormal) {
- if (nbits == fpi->nbits - 1
- && x[nbits >> kshift] & 1 << (nbits & kmask))
- irv = STRTOG_Normal;
- }
- else if (b->wds > k
- || ((n = nbits & kmask) !=0
- && hi0bits(x[k-1]) < 32-n)) {
- rshift(b,1);
- if (++e > fpi->emax)
- goto ovfl;
- }
- irv |= STRTOG_Inexhi;
- }
- else
- irv |= STRTOG_Inexlo;
- }
- *bp = b;
- *exp = e;
- return irv;
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- void
-#ifdef KR_headers
-rshift(b, k) Bigint *b; int k;
-#else
-rshift(Bigint *b, int k)
-#endif
-{
- ULong *x, *x1, *xe, y;
- int n;
-
- x = x1 = b->x;
- n = k >> kshift;
- if (n < b->wds) {
- xe = x + b->wds;
- x += n;
- if (k &= kmask) {
- n = ULbits - k;
- y = *x++ >> k;
- while(x < xe) {
- *x1++ = (y | (*x << n)) & ALL_ON;
- y = *x++ >> k;
- }
- if ((*x1 = y) !=0)
- x1++;
- }
- else
- while(x < xe)
- *x1++ = *x++;
- }
- if ((b->wds = x1 - b->x) == 0)
- b->x[0] = 0;
- }
-
- int
-#ifdef KR_headers
-trailz(b) Bigint *b;
-#else
-trailz(Bigint *b)
-#endif
-{
- ULong L, *x, *xe;
- int n = 0;
-
- x = b->x;
- xe = x + b->wds;
- for(n = 0; x < xe && !*x; x++)
- n += ULbits;
- if (x < xe) {
- L = *x;
- n += lo0bits(&L);
- }
- return n;
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 2000 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- unsigned char hexdig[256];
-
- static void
-#ifdef KR_headers
-htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
-#else
-htinit(unsigned char *h, unsigned char *s, int inc)
-#endif
-{
- int i, j;
- for(i = 0; (j = s[i]) !=0; i++)
- h[j] = i + inc;
- }
-
- void
-hexdig_init_D2A(Void)
-{
-#define USC (unsigned char *)
- htinit(hexdig, USC "0123456789", 0x10);
- htinit(hexdig, USC "abcdef", 0x10 + 10);
- htinit(hexdig, USC "ABCDEF", 0x10 + 10);
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 2000 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- static void
-#ifdef KR_headers
-L_shift(x, x1, i) ULong *x; ULong *x1; int i;
-#else
-L_shift(ULong *x, ULong *x1, int i)
-#endif
-{
- int j;
-
- i = 8 - i;
- i <<= 2;
- j = ULbits - i;
- do {
- *x |= x[1] << j;
- x[1] >>= i;
- } while(++x < x1);
- }
-
- int
-#ifdef KR_headers
-hexnan(sp, fpi, x0)
- CONST char **sp; FPI *fpi; ULong *x0;
-#else
-hexnan( CONST char **sp, FPI *fpi, ULong *x0)
-#endif
-{
- ULong c, h, *x, *x1, *xe;
- CONST char *s;
- int havedig, hd0, i, nbits;
-
- if (!hexdig['0'])
- hexdig_init_D2A();
- nbits = fpi->nbits;
- x = x0 + (nbits >> kshift);
- if (nbits & kmask)
- x++;
- *--x = 0;
- x1 = xe = x;
- havedig = hd0 = i = 0;
- s = *sp;
- /* allow optional initial 0x or 0X */
- while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
- ++s;
- if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
- && *(CONST unsigned char*)(s+3) > ' ')
- s += 2;
- while((c = *(CONST unsigned char*)++s)) {
- if (!(h = hexdig[c])) {
- if (c <= ' ') {
- if (hd0 < havedig) {
- if (x < x1 && i < 8)
- L_shift(x, x1, i);
- if (x <= x0) {
- i = 8;
- continue;
- }
- hd0 = havedig;
- *--x = 0;
- x1 = x;
- i = 0;
- }
- while(*(CONST unsigned char*)(s+1) <= ' ')
- ++s;
- if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X')
- && *(CONST unsigned char*)(s+3) > ' ')
- s += 2;
- continue;
- }
- if (/*(*/ c == ')' && havedig) {
- *sp = s + 1;
- break;
- }
-#ifndef GDTOA_NON_PEDANTIC_NANCHECK
- do {
- if (/*(*/ c == ')') {
- *sp = s + 1;
- break;
- }
- } while((c = *++s));
-#endif
- return STRTOG_NaN;
- }
- havedig++;
- if (++i > 8) {
- if (x <= x0)
- continue;
- i = 1;
- *--x = 0;
- }
- *x = (*x << 4) | (h & 0xf);
- }
- if (!havedig)
- return STRTOG_NaN;
- if (x < x1 && i < 8)
- L_shift(x, x1, i);
- if (x > x0) {
- x1 = x0;
- do *x1++ = *x++;
- while(x <= xe);
- do *x1++ = 0;
- while(x1 <= xe);
- }
- else {
- /* truncate high-order word if necessary */
- if ( (i = nbits & (ULbits-1)) !=0)
- *xe &= ((ULong)0xffffffff) >> (ULbits - i);
- }
- for(x1 = xe;; --x1) {
- if (*x1 != 0)
- break;
- if (x1 == x0) {
- *x1 = 1;
- break;
- }
- }
- return STRTOG_NaNbits;
- }
+++ /dev/null
-# /****************************************************************
-# Copyright (C) 1998 by Lucent Technologies
-# All Rights Reserved
-#
-# Permission to use, copy, modify, and distribute this software and
-# its documentation for any purpose and without fee is hereby
-# granted, provided that the above copyright notice appear in all
-# copies and that both that the copyright notice and this
-# permission notice and warranty disclaimer appear in supporting
-# documentation, and that the name of Lucent or any of its entities
-# not be used in advertising or publicity pertaining to
-# distribution of the software without specific, written prior
-# permission.
-#
-# LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-# INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-# IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-# SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-# IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-# ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-# THIS SOFTWARE.
-#
-# ****************************************************************/
-
-.SUFFIXES: .c .o
-CC = cc
-CFLAGS = -g
-
-.c.o:
- $(CC) -c $(CFLAGS) $*.c
-
-all: gdtoa.a
-
-gdtoa.a: dmisc.c g_Qfmt.c g__fmt.c gdtoa.c gethex.c gmisc.c hd_init.c\
- hexnan.c misc.c smisc.c strtopQ.c strtodg.c sum.c ulp.c
- $(CC) -c $(CFLAGS) $?
- x=`echo $? | sed 's/\.c/.o/g'` && ar ruv gdtoa.a $$x && rm $$x
- ranlib gdtoa.a || true
-
-clean:
- rm -f *.[ao]
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998, 1999 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- static Bigint *freelist[Kmax+1];
-#ifndef Omit_Private_Memory
-#ifndef PRIVATE_MEM
-#define PRIVATE_MEM 2304
-#endif
-#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
-static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
-#endif
-
- Bigint *
-Balloc
-#ifdef KR_headers
- (k) int k;
-#else
- (int k)
-#endif
-{
- int x;
- Bigint *rv;
-#ifndef Omit_Private_Memory
- unsigned int len;
-#endif
-
- ACQUIRE_DTOA_LOCK(0);
- /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
- /* but this case seems very unlikely. */
- if (k <= Kmax && (rv = freelist[k]) !=0) {
- freelist[k] = rv->next;
- }
- else {
- x = 1 << k;
-#ifdef Omit_Private_Memory
- rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
-#else
- len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
- /sizeof(double);
- if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
- rv = (Bigint*)pmem_next;
- pmem_next += len;
- }
- else
- rv = (Bigint*)MALLOC(len*sizeof(double));
-#endif
- rv->k = k;
- rv->maxwds = x;
- }
- FREE_DTOA_LOCK(0);
- rv->sign = rv->wds = 0;
- return rv;
- }
-
- void
-Bfree
-#ifdef KR_headers
- (v) Bigint *v;
-#else
- (Bigint *v)
-#endif
-{
- if (v) {
- if (v->k > Kmax)
-#ifdef FREE
- FREE((void*)v);
-#else
- free((void*)v);
-#endif
- else {
- ACQUIRE_DTOA_LOCK(0);
- v->next = freelist[v->k];
- freelist[v->k] = v;
- FREE_DTOA_LOCK(0);
- }
- }
- }
-
- int
-lo0bits
-#ifdef KR_headers
- (y) ULong *y;
-#else
- (ULong *y)
-#endif
-{
- int k;
- ULong x = *y;
-
- if (x & 7) {
- if (x & 1)
- return 0;
- if (x & 2) {
- *y = x >> 1;
- return 1;
- }
- *y = x >> 2;
- return 2;
- }
- k = 0;
- if (!(x & 0xffff)) {
- k = 16;
- x >>= 16;
- }
- if (!(x & 0xff)) {
- k += 8;
- x >>= 8;
- }
- if (!(x & 0xf)) {
- k += 4;
- x >>= 4;
- }
- if (!(x & 0x3)) {
- k += 2;
- x >>= 2;
- }
- if (!(x & 1)) {
- k++;
- x >>= 1;
- if (!x)
- return 32;
- }
- *y = x;
- return k;
- }
-
- Bigint *
-multadd
-#ifdef KR_headers
- (b, m, a) Bigint *b; int m, a;
-#else
- (Bigint *b, int m, int a) /* multiply by m and add a */
-#endif
-{
- int i, wds;
-#ifdef ULLong
- ULong *x;
- ULLong carry, y;
-#else
- ULong carry, *x, y;
-#ifdef Pack_32
- ULong xi, z;
-#endif
-#endif
- Bigint *b1;
-
- wds = b->wds;
- x = b->x;
- i = 0;
- carry = a;
- do {
-#ifdef ULLong
- y = *x * (ULLong)m + carry;
- carry = y >> 32;
- *x++ = y & 0xffffffffUL;
-#else
-#ifdef Pack_32
- xi = *x;
- y = (xi & 0xffff) * m + carry;
- z = (xi >> 16) * m + (y >> 16);
- carry = z >> 16;
- *x++ = (z << 16) + (y & 0xffff);
-#else
- y = *x * m + carry;
- carry = y >> 16;
- *x++ = y & 0xffff;
-#endif
-#endif
- }
- while(++i < wds);
- if (carry) {
- if (wds >= b->maxwds) {
- b1 = Balloc(b->k+1);
- Bcopy(b1, b);
- Bfree(b);
- b = b1;
- }
- b->x[wds++] = carry;
- b->wds = wds;
- }
- return b;
- }
-
- int
-hi0bits_D2A
-#ifdef KR_headers
- (x) ULong x;
-#else
- (ULong x)
-#endif
-{
- int k = 0;
-
- if (!(x & 0xffff0000)) {
- k = 16;
- x <<= 16;
- }
- if (!(x & 0xff000000)) {
- k += 8;
- x <<= 8;
- }
- if (!(x & 0xf0000000)) {
- k += 4;
- x <<= 4;
- }
- if (!(x & 0xc0000000)) {
- k += 2;
- x <<= 2;
- }
- if (!(x & 0x80000000)) {
- k++;
- if (!(x & 0x40000000))
- return 32;
- }
- return k;
- }
-
- Bigint *
-i2b
-#ifdef KR_headers
- (i) int i;
-#else
- (int i)
-#endif
-{
- Bigint *b;
-
- b = Balloc(1);
- b->x[0] = i;
- b->wds = 1;
- return b;
- }
-
- Bigint *
-mult
-#ifdef KR_headers
- (a, b) Bigint *a, *b;
-#else
- (Bigint *a, Bigint *b)
-#endif
-{
- Bigint *c;
- int k, wa, wb, wc;
- ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
- ULong y;
-#ifdef ULLong
- ULLong carry, z;
-#else
- ULong carry, z;
-#ifdef Pack_32
- ULong z2;
-#endif
-#endif
-
- if (a->wds < b->wds) {
- c = a;
- a = b;
- b = c;
- }
- k = a->k;
- wa = a->wds;
- wb = b->wds;
- wc = wa + wb;
- if (wc > a->maxwds)
- k++;
- c = Balloc(k);
- for(x = c->x, xa = x + wc; x < xa; x++)
- *x = 0;
- xa = a->x;
- xae = xa + wa;
- xb = b->x;
- xbe = xb + wb;
- xc0 = c->x;
-#ifdef ULLong
- for(; xb < xbe; xc0++) {
- if ( (y = *xb++) !=0) {
- x = xa;
- xc = xc0;
- carry = 0;
- do {
- z = *x++ * (ULLong)y + *xc + carry;
- carry = z >> 32;
- *xc++ = z & 0xffffffffUL;
- }
- while(x < xae);
- *xc = carry;
- }
- }
-#else
-#ifdef Pack_32
- for(; xb < xbe; xb++, xc0++) {
- if ( (y = *xb & 0xffff) !=0) {
- x = xa;
- xc = xc0;
- carry = 0;
- do {
- z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
- carry = z >> 16;
- z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
- carry = z2 >> 16;
- Storeinc(xc, z2, z);
- }
- while(x < xae);
- *xc = carry;
- }
- if ( (y = *xb >> 16) !=0) {
- x = xa;
- xc = xc0;
- carry = 0;
- z2 = *xc;
- do {
- z = (*x & 0xffff) * y + (*xc >> 16) + carry;
- carry = z >> 16;
- Storeinc(xc, z, z2);
- z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
- carry = z2 >> 16;
- }
- while(x < xae);
- *xc = z2;
- }
- }
-#else
- for(; xb < xbe; xc0++) {
- if ( (y = *xb++) !=0) {
- x = xa;
- xc = xc0;
- carry = 0;
- do {
- z = *x++ * y + *xc + carry;
- carry = z >> 16;
- *xc++ = z & 0xffff;
- }
- while(x < xae);
- *xc = carry;
- }
- }
-#endif
-#endif
- for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
- c->wds = wc;
- return c;
- }
-
- static Bigint *p5s;
-
- Bigint *
-pow5mult
-#ifdef KR_headers
- (b, k) Bigint *b; int k;
-#else
- (Bigint *b, int k)
-#endif
-{
- Bigint *b1, *p5, *p51;
- int i;
- static int p05[3] = { 5, 25, 125 };
-
- if ( (i = k & 3) !=0)
- b = multadd(b, p05[i-1], 0);
-
- if (!(k >>= 2))
- return b;
- if ((p5 = p5s) == 0) {
- /* first time */
-#ifdef MULTIPLE_THREADS
- ACQUIRE_DTOA_LOCK(1);
- if (!(p5 = p5s)) {
- p5 = p5s = i2b(625);
- p5->next = 0;
- }
- FREE_DTOA_LOCK(1);
-#else
- p5 = p5s = i2b(625);
- p5->next = 0;
-#endif
- }
- for(;;) {
- if (k & 1) {
- b1 = mult(b, p5);
- Bfree(b);
- b = b1;
- }
- if (!(k >>= 1))
- break;
- if ((p51 = p5->next) == 0) {
-#ifdef MULTIPLE_THREADS
- ACQUIRE_DTOA_LOCK(1);
- if (!(p51 = p5->next)) {
- p51 = p5->next = mult(p5,p5);
- p51->next = 0;
- }
- FREE_DTOA_LOCK(1);
-#else
- p51 = p5->next = mult(p5,p5);
- p51->next = 0;
-#endif
- }
- p5 = p51;
- }
- return b;
- }
-
- Bigint *
-lshift
-#ifdef KR_headers
- (b, k) Bigint *b; int k;
-#else
- (Bigint *b, int k)
-#endif
-{
- int i, k1, n, n1;
- Bigint *b1;
- ULong *x, *x1, *xe, z;
-
- n = k >> kshift;
- k1 = b->k;
- n1 = n + b->wds + 1;
- for(i = b->maxwds; n1 > i; i <<= 1)
- k1++;
- b1 = Balloc(k1);
- x1 = b1->x;
- for(i = 0; i < n; i++)
- *x1++ = 0;
- x = b->x;
- xe = x + b->wds;
- if (k &= kmask) {
-#ifdef Pack_32
- k1 = 32 - k;
- z = 0;
- do {
- *x1++ = *x << k | z;
- z = *x++ >> k1;
- }
- while(x < xe);
- if ((*x1 = z) !=0)
- ++n1;
-#else
- k1 = 16 - k;
- z = 0;
- do {
- *x1++ = *x << k & 0xffff | z;
- z = *x++ >> k1;
- }
- while(x < xe);
- if (*x1 = z)
- ++n1;
-#endif
- }
- else do
- *x1++ = *x++;
- while(x < xe);
- b1->wds = n1 - 1;
- Bfree(b);
- return b1;
- }
-
- int
-cmp
-#ifdef KR_headers
- (a, b) Bigint *a, *b;
-#else
- (Bigint *a, Bigint *b)
-#endif
-{
- ULong *xa, *xa0, *xb, *xb0;
- int i, j;
-
- i = a->wds;
- j = b->wds;
-#ifdef DEBUG
- if (i > 1 && !a->x[i-1])
- Bug("cmp called with a->x[a->wds-1] == 0");
- if (j > 1 && !b->x[j-1])
- Bug("cmp called with b->x[b->wds-1] == 0");
-#endif
- if (i -= j)
- return i;
- xa0 = a->x;
- xa = xa0 + j;
- xb0 = b->x;
- xb = xb0 + j;
- for(;;) {
- if (*--xa != *--xb)
- return *xa < *xb ? -1 : 1;
- if (xa <= xa0)
- break;
- }
- return 0;
- }
-
- Bigint *
-diff
-#ifdef KR_headers
- (a, b) Bigint *a, *b;
-#else
- (Bigint *a, Bigint *b)
-#endif
-{
- Bigint *c;
- int i, wa, wb;
- ULong *xa, *xae, *xb, *xbe, *xc;
-#ifdef ULLong
- ULLong borrow, y;
-#else
- ULong borrow, y;
-#ifdef Pack_32
- ULong z;
-#endif
-#endif
-
- i = cmp(a,b);
- if (!i) {
- c = Balloc(0);
- c->wds = 1;
- c->x[0] = 0;
- return c;
- }
- if (i < 0) {
- c = a;
- a = b;
- b = c;
- i = 1;
- }
- else
- i = 0;
- c = Balloc(a->k);
- c->sign = i;
- wa = a->wds;
- xa = a->x;
- xae = xa + wa;
- wb = b->wds;
- xb = b->x;
- xbe = xb + wb;
- xc = c->x;
- borrow = 0;
-#ifdef ULLong
- do {
- y = (ULLong)*xa++ - *xb++ - borrow;
- borrow = y >> 32 & 1UL;
- *xc++ = y & 0xffffffffUL;
- }
- while(xb < xbe);
- while(xa < xae) {
- y = *xa++ - borrow;
- borrow = y >> 32 & 1UL;
- *xc++ = y & 0xffffffffUL;
- }
-#else
-#ifdef Pack_32
- do {
- y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
- while(xb < xbe);
- while(xa < xae) {
- y = (*xa & 0xffff) - borrow;
- borrow = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) - borrow;
- borrow = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
-#else
- do {
- y = *xa++ - *xb++ - borrow;
- borrow = (y & 0x10000) >> 16;
- *xc++ = y & 0xffff;
- }
- while(xb < xbe);
- while(xa < xae) {
- y = *xa++ - borrow;
- borrow = (y & 0x10000) >> 16;
- *xc++ = y & 0xffff;
- }
-#endif
-#endif
- while(!*--xc)
- wa--;
- c->wds = wa;
- return c;
- }
-
- double
-b2d
-#ifdef KR_headers
- (a, e) Bigint *a; int *e;
-#else
- (Bigint *a, int *e)
-#endif
-{
- ULong *xa, *xa0, w, y, z;
- int k;
- U d;
-#ifdef VAX
- ULong d0, d1;
-#else
-#define d0 word0(&d)
-#define d1 word1(&d)
-#endif
-
- xa0 = a->x;
- xa = xa0 + a->wds;
- y = *--xa;
-#ifdef DEBUG
- if (!y) Bug("zero y in b2d");
-#endif
- k = hi0bits(y);
- *e = 32 - k;
-#ifdef Pack_32
- if (k < Ebits) {
- d0 = Exp_1 | y >> (Ebits - k);
- w = xa > xa0 ? *--xa : 0;
- d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
- goto ret_d;
- }
- z = xa > xa0 ? *--xa : 0;
- if (k -= Ebits) {
- d0 = Exp_1 | y << k | z >> (32 - k);
- y = xa > xa0 ? *--xa : 0;
- d1 = z << k | y >> (32 - k);
- }
- else {
- d0 = Exp_1 | y;
- d1 = z;
- }
-#else
- if (k < Ebits + 16) {
- z = xa > xa0 ? *--xa : 0;
- d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
- w = xa > xa0 ? *--xa : 0;
- y = xa > xa0 ? *--xa : 0;
- d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
- goto ret_d;
- }
- z = xa > xa0 ? *--xa : 0;
- w = xa > xa0 ? *--xa : 0;
- k -= Ebits + 16;
- d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
- y = xa > xa0 ? *--xa : 0;
- d1 = w << k + 16 | y << k;
-#endif
- ret_d:
-#ifdef VAX
- word0(&d) = d0 >> 16 | d0 << 16;
- word1(&d) = d1 >> 16 | d1 << 16;
-#endif
- return dval(&d);
- }
-#undef d0
-#undef d1
-
- Bigint *
-d2b
-#ifdef KR_headers
- (dd, e, bits) double dd; int *e, *bits;
-#else
- (double dd, int *e, int *bits)
-#endif
-{
- Bigint *b;
- U d;
-#ifndef Sudden_Underflow
- int i;
-#endif
- int de, k;
- ULong *x, y, z;
-#ifdef VAX
- ULong d0, d1;
-#else
-#define d0 word0(&d)
-#define d1 word1(&d)
-#endif
- d.d = dd;
-#ifdef VAX
- d0 = word0(&d) >> 16 | word0(&d) << 16;
- d1 = word1(&d) >> 16 | word1(&d) << 16;
-#endif
-
-#ifdef Pack_32
- b = Balloc(1);
-#else
- b = Balloc(2);
-#endif
- x = b->x;
-
- z = d0 & Frac_mask;
- d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
-#ifdef Sudden_Underflow
- de = (int)(d0 >> Exp_shift);
-#ifndef IBM
- z |= Exp_msk11;
-#endif
-#else
- if ( (de = (int)(d0 >> Exp_shift)) !=0)
- z |= Exp_msk1;
-#endif
-#ifdef Pack_32
- if ( (y = d1) !=0) {
- if ( (k = lo0bits(&y)) !=0) {
- x[0] = y | z << (32 - k);
- z >>= k;
- }
- else
- x[0] = y;
-#ifndef Sudden_Underflow
- i =
-#endif
- b->wds = (x[1] = z) !=0 ? 2 : 1;
- }
- else {
- k = lo0bits(&z);
- x[0] = z;
-#ifndef Sudden_Underflow
- i =
-#endif
- b->wds = 1;
- k += 32;
- }
-#else
- if ( (y = d1) !=0) {
- if ( (k = lo0bits(&y)) !=0)
- if (k >= 16) {
- x[0] = y | z << 32 - k & 0xffff;
- x[1] = z >> k - 16 & 0xffff;
- x[2] = z >> k;
- i = 2;
- }
- else {
- x[0] = y & 0xffff;
- x[1] = y >> 16 | z << 16 - k & 0xffff;
- x[2] = z >> k & 0xffff;
- x[3] = z >> k+16;
- i = 3;
- }
- else {
- x[0] = y & 0xffff;
- x[1] = y >> 16;
- x[2] = z & 0xffff;
- x[3] = z >> 16;
- i = 3;
- }
- }
- else {
-#ifdef DEBUG
- if (!z)
- Bug("Zero passed to d2b");
-#endif
- k = lo0bits(&z);
- if (k >= 16) {
- x[0] = z;
- i = 0;
- }
- else {
- x[0] = z & 0xffff;
- x[1] = z >> 16;
- i = 1;
- }
- k += 32;
- }
- while(!x[i])
- --i;
- b->wds = i + 1;
-#endif
-#ifndef Sudden_Underflow
- if (de) {
-#endif
-#ifdef IBM
- *e = (de - Bias - (P-1) << 2) + k;
- *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
-#else
- *e = de - Bias - (P-1) + k;
- *bits = P - k;
-#endif
-#ifndef Sudden_Underflow
- }
- else {
- *e = de - Bias - (P-1) + 1 + k;
-#ifdef Pack_32
- *bits = 32*i - hi0bits(x[i-1]);
-#else
- *bits = (i+2)*16 - hi0bits(x[i]);
-#endif
- }
-#endif
- return b;
- }
-#undef d0
-#undef d1
-
- CONST double
-#ifdef IEEE_Arith
-bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
-CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
- };
-#else
-#ifdef IBM
-bigtens[] = { 1e16, 1e32, 1e64 };
-CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
-#else
-bigtens[] = { 1e16, 1e32 };
-CONST double tinytens[] = { 1e-16, 1e-32 };
-#endif
-#endif
-
- CONST double
-tens[] = {
- 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
- 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
- 1e20, 1e21, 1e22
-#ifdef VAX
- , 1e23, 1e24
-#endif
- };
-
-#ifdef NO_STRING_H
- char *
-#ifdef KR_headers
-strcpy_D2A(a, b) char *a; char *b;
-#else
-strcpy_D2A(char *a, CONST char *b)
-#endif
-{
- while((*a = *b++))
- a++;
- return a;
- }
-
-
- Char *
-#ifdef KR_headers
-memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
-#else
-memcpy_D2A(void *a1, void *b1, size_t len)
-#endif
-{
- char *a = (char*)a1, *ae = a + len;
- char *b = (char*)b1, *a0 = a;
- while(a < ae)
- *a++ = *b++;
- return a0;
- }
-
-#endif /* NO_STRING_H */
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998, 1999 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- Bigint *
-s2b
-#ifdef KR_headers
- (s, nd0, nd, y9, dplen) CONST char *s; int dplen, nd0, nd; ULong y9;
-#else
- (CONST char *s, int nd0, int nd, ULong y9, int dplen)
-#endif
-{
- Bigint *b;
- int i, k;
- Long x, y;
-
- x = (nd + 8) / 9;
- for(k = 0, y = 1; x > y; y <<= 1, k++) ;
-#ifdef Pack_32
- b = Balloc(k);
- b->x[0] = y9;
- b->wds = 1;
-#else
- b = Balloc(k+1);
- b->x[0] = y9 & 0xffff;
- b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
-#endif
-
- i = 9;
- if (9 < nd0) {
- s += 9;
- do b = multadd(b, 10, *s++ - '0');
- while(++i < nd0);
- s += dplen;
- }
- else
- s += dplen + 9;
- for(; i < nd; i++)
- b = multadd(b, 10, *s++ - '0');
- return b;
- }
-
- double
-ratio
-#ifdef KR_headers
- (a, b) Bigint *a, *b;
-#else
- (Bigint *a, Bigint *b)
-#endif
-{
- U da, db;
- int k, ka, kb;
-
- dval(&da) = b2d(a, &ka);
- dval(&db) = b2d(b, &kb);
- k = ka - kb + ULbits*(a->wds - b->wds);
-#ifdef IBM
- if (k > 0) {
- word0(&da) += (k >> 2)*Exp_msk1;
- if (k &= 3)
- dval(&da) *= 1 << k;
- }
- else {
- k = -k;
- word0(&db) += (k >> 2)*Exp_msk1;
- if (k &= 3)
- dval(&db) *= 1 << k;
- }
-#else
- if (k > 0)
- word0(&da) += k*Exp_msk1;
- else {
- k = -k;
- word0(&db) += k*Exp_msk1;
- }
-#endif
- return dval(&da) / dval(&db);
- }
-
-#ifdef INFNAN_CHECK
-
- int
-match
-#ifdef KR_headers
- (sp, t) char **sp, *t;
-#else
- (CONST char **sp, char *t)
-#endif
-{
- int c, d;
- CONST char *s = *sp;
-
- while( (d = *t++) !=0) {
- if ((c = *++s) >= 'A' && c <= 'Z')
- c += 'a' - 'A';
- if (c != d)
- return 0;
- }
- *sp = s + 1;
- return 1;
- }
-#endif /* INFNAN_CHECK */
-
- void
-#ifdef KR_headers
-copybits(c, n, b) ULong *c; int n; Bigint *b;
-#else
-copybits(ULong *c, int n, Bigint *b)
-#endif
-{
- ULong *ce, *x, *xe;
-#ifdef Pack_16
- int nw, nw1;
-#endif
-
- ce = c + ((n-1) >> kshift) + 1;
- x = b->x;
-#ifdef Pack_32
- xe = x + b->wds;
- while(x < xe)
- *c++ = *x++;
-#else
- nw = b->wds;
- nw1 = nw & 1;
- for(xe = x + (nw - nw1); x < xe; x += 2)
- Storeinc(c, x[1], x[0]);
- if (nw1)
- *c++ = *x;
-#endif
- while(c < ce)
- *c++ = 0;
- }
-
- ULong
-#ifdef KR_headers
-any_on(b, k) Bigint *b; int k;
-#else
-any_on(Bigint *b, int k)
-#endif
-{
- int n, nwds;
- ULong *x, *x0, x1, x2;
-
- x = b->x;
- nwds = b->wds;
- n = k >> kshift;
- if (n > nwds)
- n = nwds;
- else if (n < nwds && (k &= kmask)) {
- x1 = x2 = x[n];
- x1 >>= k;
- x1 <<= k;
- if (x1 != x2)
- return 1;
- }
- x0 = x;
- x += n;
- while(x > x0)
- if (*--x)
- return 1;
- return 0;
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998-2001 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
-#ifdef USE_LOCALE
-#include "locale.h"
-#endif
-
- static CONST int
-fivesbits[] = { 0, 3, 5, 7, 10, 12, 14, 17, 19, 21,
- 24, 26, 28, 31, 33, 35, 38, 40, 42, 45,
- 47, 49, 52
-#ifdef VAX
- , 54, 56
-#endif
- };
-
- Bigint *
-#ifdef KR_headers
-increment(b) Bigint *b;
-#else
-increment(Bigint *b)
-#endif
-{
- ULong *x, *xe;
- Bigint *b1;
-#ifdef Pack_16
- ULong carry = 1, y;
-#endif
-
- x = b->x;
- xe = x + b->wds;
-#ifdef Pack_32
- do {
- if (*x < (ULong)0xffffffffL) {
- ++*x;
- return b;
- }
- *x++ = 0;
- } while(x < xe);
-#else
- do {
- y = *x + carry;
- carry = y >> 16;
- *x++ = y & 0xffff;
- if (!carry)
- return b;
- } while(x < xe);
- if (carry)
-#endif
- {
- if (b->wds >= b->maxwds) {
- b1 = Balloc(b->k+1);
- Bcopy(b1,b);
- Bfree(b);
- b = b1;
- }
- b->x[b->wds++] = 1;
- }
- return b;
- }
-
- void
-#ifdef KR_headers
-decrement(b) Bigint *b;
-#else
-decrement(Bigint *b)
-#endif
-{
- ULong *x, *xe;
-#ifdef Pack_16
- ULong borrow = 1, y;
-#endif
-
- x = b->x;
- xe = x + b->wds;
-#ifdef Pack_32
- do {
- if (*x) {
- --*x;
- break;
- }
- *x++ = 0xffffffffL;
- }
- while(x < xe);
-#else
- do {
- y = *x - borrow;
- borrow = (y & 0x10000) >> 16;
- *x++ = y & 0xffff;
- } while(borrow && x < xe);
-#endif
- }
-
- static int
-#ifdef KR_headers
-all_on(b, n) Bigint *b; int n;
-#else
-all_on(Bigint *b, int n)
-#endif
-{
- ULong *x, *xe;
-
- x = b->x;
- xe = x + (n >> kshift);
- while(x < xe)
- if ((*x++ & ALL_ON) != ALL_ON)
- return 0;
- if (n &= kmask)
- return ((*x | (ALL_ON << n)) & ALL_ON) == ALL_ON;
- return 1;
- }
-
- Bigint *
-#ifdef KR_headers
-set_ones(b, n) Bigint *b; int n;
-#else
-set_ones(Bigint *b, int n)
-#endif
-{
- int k;
- ULong *x, *xe;
-
- k = (n + ((1 << kshift) - 1)) >> kshift;
- if (b->k < k) {
- Bfree(b);
- b = Balloc(k);
- }
- k = n >> kshift;
- if (n &= kmask)
- k++;
- b->wds = k;
- x = b->x;
- xe = x + k;
- while(x < xe)
- *x++ = ALL_ON;
- if (n)
- x[-1] >>= ULbits - n;
- return b;
- }
-
- static int
-rvOK
-#ifdef KR_headers
- (d, fpi, exp, bits, exact, rd, irv)
- U *d; FPI *fpi; Long *exp; ULong *bits; int exact, rd, *irv;
-#else
- (U *d, FPI *fpi, Long *exp, ULong *bits, int exact, int rd, int *irv)
-#endif
-{
- Bigint *b;
- ULong carry, inex, lostbits;
- int bdif, e, j, k, k1, nb, rv;
-
- carry = rv = 0;
- b = d2b(dval(d), &e, &bdif);
- bdif -= nb = fpi->nbits;
- e += bdif;
- if (bdif <= 0) {
- if (exact)
- goto trunc;
- goto ret;
- }
- if (P == nb) {
- if (
-#ifndef IMPRECISE_INEXACT
- exact &&
-#endif
- fpi->rounding ==
-#ifdef RND_PRODQUOT
- FPI_Round_near
-#else
- Flt_Rounds
-#endif
- ) goto trunc;
- goto ret;
- }
- switch(rd) {
- case 1: /* round down (toward -Infinity) */
- goto trunc;
- case 2: /* round up (toward +Infinity) */
- break;
- default: /* round near */
- k = bdif - 1;
- if (k < 0)
- goto trunc;
- if (!k) {
- if (!exact)
- goto ret;
- if (b->x[0] & 2)
- break;
- goto trunc;
- }
- if (b->x[k>>kshift] & ((ULong)1 << (k & kmask)))
- break;
- goto trunc;
- }
- /* "break" cases: round up 1 bit, then truncate; bdif > 0 */
- carry = 1;
- trunc:
- inex = lostbits = 0;
- if (bdif > 0) {
- if ( (lostbits = any_on(b, bdif)) !=0)
- inex = STRTOG_Inexlo;
- rshift(b, bdif);
- if (carry) {
- inex = STRTOG_Inexhi;
- b = increment(b);
- if ( (j = nb & kmask) !=0)
- j = ULbits - j;
- if (hi0bits(b->x[b->wds - 1]) != j) {
- if (!lostbits)
- lostbits = b->x[0] & 1;
- rshift(b, 1);
- e++;
- }
- }
- }
- else if (bdif < 0)
- b = lshift(b, -bdif);
- if (e < fpi->emin) {
- k = fpi->emin - e;
- e = fpi->emin;
- if (k > nb || fpi->sudden_underflow) {
- b->wds = inex = 0;
- *irv = STRTOG_Underflow | STRTOG_Inexlo;
- }
- else {
- k1 = k - 1;
- if (k1 > 0 && !lostbits)
- lostbits = any_on(b, k1);
- if (!lostbits && !exact)
- goto ret;
- lostbits |=
- carry = b->x[k1>>kshift] & (1 << (k1 & kmask));
- rshift(b, k);
- *irv = STRTOG_Denormal;
- if (carry) {
- b = increment(b);
- inex = STRTOG_Inexhi | STRTOG_Underflow;
- }
- else if (lostbits)
- inex = STRTOG_Inexlo | STRTOG_Underflow;
- }
- }
- else if (e > fpi->emax) {
- e = fpi->emax + 1;
- *irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- b->wds = inex = 0;
- }
- *exp = e;
- copybits(bits, nb, b);
- *irv |= inex;
- rv = 1;
- ret:
- Bfree(b);
- return rv;
- }
-
- static int
-#ifdef KR_headers
-mantbits(d) U *d;
-#else
-mantbits(U *d)
-#endif
-{
- ULong L;
-#ifdef VAX
- L = word1(d) << 16 | word1(d) >> 16;
- if (L)
-#else
- if ( (L = word1(d)) !=0)
-#endif
- return P - lo0bits(&L);
-#ifdef VAX
- L = word0(d) << 16 | word0(d) >> 16 | Exp_msk11;
-#else
- L = word0(d) | Exp_msk1;
-#endif
- return P - 32 - lo0bits(&L);
- }
-
- int
-strtodg
-#ifdef KR_headers
- (s00, se, fpi, exp, bits)
- CONST char *s00; char **se; FPI *fpi; Long *exp; ULong *bits;
-#else
- (CONST char *s00, char **se, FPI *fpi, Long *exp, ULong *bits)
-#endif
-{
- int abe, abits, asub;
- int bb0, bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, denorm;
- int dsign, e, e1, e2, emin, esign, finished, i, inex, irv;
- int j, k, nbits, nd, nd0, nf, nz, nz0, rd, rvbits, rve, rve1, sign;
- int sudden_underflow;
- CONST char *s, *s0, *s1;
- double adj0, tol;
- Long L;
- U adj, rv;
- ULong *b, *be, y, z;
- Bigint *ab, *bb, *bb1, *bd, *bd0, *bs, *delta, *rvb, *rvb0;
-#ifdef USE_LOCALE /*{{*/
-#ifdef NO_LOCALE_CACHE
- char *decimalpoint = localeconv()->decimal_point;
- int dplen = strlen(decimalpoint);
-#else
- char *decimalpoint;
- static char *decimalpoint_cache;
- static int dplen;
- if (!(s0 = decimalpoint_cache)) {
- s0 = localeconv()->decimal_point;
- if ((decimalpoint_cache = (char*)MALLOC(strlen(s0) + 1))) {
- strcpy(decimalpoint_cache, s0);
- s0 = decimalpoint_cache;
- }
- dplen = strlen(s0);
- }
- decimalpoint = (char*)s0;
-#endif /*NO_LOCALE_CACHE*/
-#else /*USE_LOCALE}{*/
-#define dplen 1
-#endif /*USE_LOCALE}}*/
-
- irv = STRTOG_Zero;
- denorm = sign = nz0 = nz = 0;
- dval(&rv) = 0.;
- rvb = 0;
- nbits = fpi->nbits;
- for(s = s00;;s++) switch(*s) {
- case '-':
- sign = 1;
- /* no break */
- case '+':
- if (*++s)
- goto break2;
- /* no break */
- case 0:
- sign = 0;
- irv = STRTOG_NoNumber;
- s = s00;
- goto ret;
- case '\t':
- case '\n':
- case '\v':
- case '\f':
- case '\r':
- case ' ':
- continue;
- default:
- goto break2;
- }
- break2:
- if (*s == '0') {
-#ifndef NO_HEX_FP
- switch(s[1]) {
- case 'x':
- case 'X':
- irv = gethex(&s, fpi, exp, &rvb, sign);
- if (irv == STRTOG_NoNumber) {
- s = s00;
- sign = 0;
- }
- goto ret;
- }
-#endif
- nz0 = 1;
- while(*++s == '0') ;
- if (!*s)
- goto ret;
- }
- sudden_underflow = fpi->sudden_underflow;
- s0 = s;
- y = z = 0;
- for(decpt = nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
- if (nd < 9)
- y = 10*y + c - '0';
- else if (nd < 16)
- z = 10*z + c - '0';
- nd0 = nd;
-#ifdef USE_LOCALE
- if (c == *decimalpoint) {
- for(i = 1; decimalpoint[i]; ++i)
- if (s[i] != decimalpoint[i])
- goto dig_done;
- s += i;
- c = *s;
-#else
- if (c == '.') {
- c = *++s;
-#endif
- decpt = 1;
- if (!nd) {
- for(; c == '0'; c = *++s)
- nz++;
- if (c > '0' && c <= '9') {
- s0 = s;
- nf += nz;
- nz = 0;
- goto have_dig;
- }
- goto dig_done;
- }
- for(; c >= '0' && c <= '9'; c = *++s) {
- have_dig:
- nz++;
- if (c -= '0') {
- nf += nz;
- for(i = 1; i < nz; i++)
- if (nd++ < 9)
- y *= 10;
- else if (nd <= DBL_DIG + 1)
- z *= 10;
- if (nd++ < 9)
- y = 10*y + c;
- else if (nd <= DBL_DIG + 1)
- z = 10*z + c;
- nz = 0;
- }
- }
- }/*}*/
- dig_done:
- e = 0;
- if (c == 'e' || c == 'E') {
- if (!nd && !nz && !nz0) {
- irv = STRTOG_NoNumber;
- s = s00;
- goto ret;
- }
- s00 = s;
- esign = 0;
- switch(c = *++s) {
- case '-':
- esign = 1;
- case '+':
- c = *++s;
- }
- if (c >= '0' && c <= '9') {
- while(c == '0')
- c = *++s;
- if (c > '0' && c <= '9') {
- L = c - '0';
- s1 = s;
- while((c = *++s) >= '0' && c <= '9')
- L = 10*L + c - '0';
- if (s - s1 > 8 || L > 19999)
- /* Avoid confusion from exponents
- * so large that e might overflow.
- */
- e = 19999; /* safe for 16 bit ints */
- else
- e = (int)L;
- if (esign)
- e = -e;
- }
- else
- e = 0;
- }
- else
- s = s00;
- }
- if (!nd) {
- if (!nz && !nz0) {
-#ifdef INFNAN_CHECK
- /* Check for Nan and Infinity */
- if (!decpt)
- switch(c) {
- case 'i':
- case 'I':
- if (match(&s,"nf")) {
- --s;
- if (!match(&s,"inity"))
- ++s;
- irv = STRTOG_Infinite;
- goto infnanexp;
- }
- break;
- case 'n':
- case 'N':
- if (match(&s, "an")) {
- irv = STRTOG_NaN;
- *exp = fpi->emax + 1;
-#ifndef No_Hex_NaN
- if (*s == '(') /*)*/
- irv = hexnan(&s, fpi, bits);
-#endif
- goto infnanexp;
- }
- }
-#endif /* INFNAN_CHECK */
- irv = STRTOG_NoNumber;
- s = s00;
- }
- goto ret;
- }
-
- irv = STRTOG_Normal;
- e1 = e -= nf;
- rd = 0;
- switch(fpi->rounding & 3) {
- case FPI_Round_up:
- rd = 2 - sign;
- break;
- case FPI_Round_zero:
- rd = 1;
- break;
- case FPI_Round_down:
- rd = 1 + sign;
- }
-
- /* Now we have nd0 digits, starting at s0, followed by a
- * decimal point, followed by nd-nd0 digits. The number we're
- * after is the integer represented by those digits times
- * 10**e */
-
- if (!nd0)
- nd0 = nd;
- k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
- dval(&rv) = y;
- if (k > 9)
- dval(&rv) = tens[k - 9] * dval(&rv) + z;
- bd0 = 0;
- if (nbits <= P && nd <= DBL_DIG) {
- if (!e) {
- if (rvOK(&rv, fpi, exp, bits, 1, rd, &irv))
- goto ret;
- }
- else if (e > 0) {
- if (e <= Ten_pmax) {
-#ifdef VAX
- goto vax_ovfl_check;
-#else
- i = fivesbits[e] + mantbits(&rv) <= P;
- /* rv = */ rounded_product(dval(&rv), tens[e]);
- if (rvOK(&rv, fpi, exp, bits, i, rd, &irv))
- goto ret;
- e1 -= e;
- goto rv_notOK;
-#endif
- }
- i = DBL_DIG - nd;
- if (e <= Ten_pmax + i) {
- /* A fancier test would sometimes let us do
- * this for larger i values.
- */
- e2 = e - i;
- e1 -= i;
- dval(&rv) *= tens[i];
-#ifdef VAX
- /* VAX exponent range is so narrow we must
- * worry about overflow here...
- */
- vax_ovfl_check:
- dval(&adj) = dval(&rv);
- word0(&adj) -= P*Exp_msk1;
- /* adj = */ rounded_product(dval(&adj), tens[e2]);
- if ((word0(&adj) & Exp_mask)
- > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
- goto rv_notOK;
- word0(&adj) += P*Exp_msk1;
- dval(&rv) = dval(&adj);
-#else
- /* rv = */ rounded_product(dval(&rv), tens[e2]);
-#endif
- if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
- goto ret;
- e1 -= e2;
- }
- }
-#ifndef Inaccurate_Divide
- else if (e >= -Ten_pmax) {
- /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
- if (rvOK(&rv, fpi, exp, bits, 0, rd, &irv))
- goto ret;
- e1 -= e;
- }
-#endif
- }
- rv_notOK:
- e1 += nd - k;
-
- /* Get starting approximation = rv * 10**e1 */
-
- e2 = 0;
- if (e1 > 0) {
- if ( (i = e1 & 15) !=0)
- dval(&rv) *= tens[i];
- if (e1 &= ~15) {
- e1 >>= 4;
- while(e1 >= (1 << (n_bigtens-1))) {
- e2 += ((word0(&rv) & Exp_mask)
- >> Exp_shift1) - Bias;
- word0(&rv) &= ~Exp_mask;
- word0(&rv) |= Bias << Exp_shift1;
- dval(&rv) *= bigtens[n_bigtens-1];
- e1 -= 1 << (n_bigtens-1);
- }
- e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
- word0(&rv) &= ~Exp_mask;
- word0(&rv) |= Bias << Exp_shift1;
- for(j = 0; e1 > 0; j++, e1 >>= 1)
- if (e1 & 1)
- dval(&rv) *= bigtens[j];
- }
- }
- else if (e1 < 0) {
- e1 = -e1;
- if ( (i = e1 & 15) !=0)
- dval(&rv) /= tens[i];
- if (e1 &= ~15) {
- e1 >>= 4;
- while(e1 >= (1 << (n_bigtens-1))) {
- e2 += ((word0(&rv) & Exp_mask)
- >> Exp_shift1) - Bias;
- word0(&rv) &= ~Exp_mask;
- word0(&rv) |= Bias << Exp_shift1;
- dval(&rv) *= tinytens[n_bigtens-1];
- e1 -= 1 << (n_bigtens-1);
- }
- e2 += ((word0(&rv) & Exp_mask) >> Exp_shift1) - Bias;
- word0(&rv) &= ~Exp_mask;
- word0(&rv) |= Bias << Exp_shift1;
- for(j = 0; e1 > 0; j++, e1 >>= 1)
- if (e1 & 1)
- dval(&rv) *= tinytens[j];
- }
- }
-#ifdef IBM
- /* e2 is a correction to the (base 2) exponent of the return
- * value, reflecting adjustments above to avoid overflow in the
- * native arithmetic. For native IBM (base 16) arithmetic, we
- * must multiply e2 by 4 to change from base 16 to 2.
- */
- e2 <<= 2;
-#endif
- rvb = d2b(dval(&rv), &rve, &rvbits); /* rv = rvb * 2^rve */
- rve += e2;
- if ((j = rvbits - nbits) > 0) {
- rshift(rvb, j);
- rvbits = nbits;
- rve += j;
- }
- bb0 = 0; /* trailing zero bits in rvb */
- e2 = rve + rvbits - nbits;
- if (e2 > fpi->emax + 1)
- goto huge;
- rve1 = rve + rvbits - nbits;
- if (e2 < (emin = fpi->emin)) {
- denorm = 1;
- j = rve - emin;
- if (j > 0) {
- rvb = lshift(rvb, j);
- rvbits += j;
- }
- else if (j < 0) {
- rvbits += j;
- if (rvbits <= 0) {
- if (rvbits < -1) {
- ufl:
- rvb->wds = 0;
- rvb->x[0] = 0;
- *exp = emin;
- irv = STRTOG_Underflow | STRTOG_Inexlo;
- goto ret;
- }
- rvb->x[0] = rvb->wds = rvbits = 1;
- }
- else
- rshift(rvb, -j);
- }
- rve = rve1 = emin;
- if (sudden_underflow && e2 + 1 < emin)
- goto ufl;
- }
-
- /* Now the hard part -- adjusting rv to the correct value.*/
-
- /* Put digits into bd: true value = bd * 10^e */
-
- bd0 = s2b(s0, nd0, nd, y, dplen);
-
- for(;;) {
- bd = Balloc(bd0->k);
- Bcopy(bd, bd0);
- bb = Balloc(rvb->k);
- Bcopy(bb, rvb);
- bbbits = rvbits - bb0;
- bbe = rve + bb0;
- bs = i2b(1);
-
- if (e >= 0) {
- bb2 = bb5 = 0;
- bd2 = bd5 = e;
- }
- else {
- bb2 = bb5 = -e;
- bd2 = bd5 = 0;
- }
- if (bbe >= 0)
- bb2 += bbe;
- else
- bd2 -= bbe;
- bs2 = bb2;
- j = nbits + 1 - bbbits;
- i = bbe + bbbits - nbits;
- if (i < emin) /* denormal */
- j += i - emin;
- bb2 += j;
- bd2 += j;
- i = bb2 < bd2 ? bb2 : bd2;
- if (i > bs2)
- i = bs2;
- if (i > 0) {
- bb2 -= i;
- bd2 -= i;
- bs2 -= i;
- }
- if (bb5 > 0) {
- bs = pow5mult(bs, bb5);
- bb1 = mult(bs, bb);
- Bfree(bb);
- bb = bb1;
- }
- bb2 -= bb0;
- if (bb2 > 0)
- bb = lshift(bb, bb2);
- else if (bb2 < 0)
- rshift(bb, -bb2);
- if (bd5 > 0)
- bd = pow5mult(bd, bd5);
- if (bd2 > 0)
- bd = lshift(bd, bd2);
- if (bs2 > 0)
- bs = lshift(bs, bs2);
- asub = 1;
- inex = STRTOG_Inexhi;
- delta = diff(bb, bd);
- if (delta->wds <= 1 && !delta->x[0])
- break;
- dsign = delta->sign;
- delta->sign = finished = 0;
- L = 0;
- i = cmp(delta, bs);
- if (rd && i <= 0) {
- irv = STRTOG_Normal;
- if ( (finished = dsign ^ (rd&1)) !=0) {
- if (dsign != 0) {
- irv |= STRTOG_Inexhi;
- goto adj1;
- }
- irv |= STRTOG_Inexlo;
- if (rve1 == emin)
- goto adj1;
- for(i = 0, j = nbits; j >= ULbits;
- i++, j -= ULbits) {
- if (rvb->x[i] & ALL_ON)
- goto adj1;
- }
- if (j > 1 && lo0bits(rvb->x + i) < j - 1)
- goto adj1;
- rve = rve1 - 1;
- rvb = set_ones(rvb, rvbits = nbits);
- break;
- }
- irv |= dsign ? STRTOG_Inexlo : STRTOG_Inexhi;
- break;
- }
- if (i < 0) {
- /* Error is less than half an ulp -- check for
- * special case of mantissa a power of two.
- */
- irv = dsign
- ? STRTOG_Normal | STRTOG_Inexlo
- : STRTOG_Normal | STRTOG_Inexhi;
- if (dsign || bbbits > 1 || denorm || rve1 == emin)
- break;
- delta = lshift(delta,1);
- if (cmp(delta, bs) > 0) {
- irv = STRTOG_Normal | STRTOG_Inexlo;
- goto drop_down;
- }
- break;
- }
- if (i == 0) {
- /* exactly half-way between */
- if (dsign) {
- if (denorm && all_on(rvb, rvbits)) {
- /*boundary case -- increment exponent*/
- rvb->wds = 1;
- rvb->x[0] = 1;
- rve = emin + nbits - (rvbits = 1);
- irv = STRTOG_Normal | STRTOG_Inexhi;
- denorm = 0;
- break;
- }
- irv = STRTOG_Normal | STRTOG_Inexlo;
- }
- else if (bbbits == 1) {
- irv = STRTOG_Normal;
- drop_down:
- /* boundary case -- decrement exponent */
- if (rve1 == emin) {
- irv = STRTOG_Normal | STRTOG_Inexhi;
- if (rvb->wds == 1 && rvb->x[0] == 1)
- sudden_underflow = 1;
- break;
- }
- rve -= nbits;
- rvb = set_ones(rvb, rvbits = nbits);
- break;
- }
- else
- irv = STRTOG_Normal | STRTOG_Inexhi;
- if ((bbbits < nbits && !denorm) || !(rvb->x[0] & 1))
- break;
- if (dsign) {
- rvb = increment(rvb);
- j = kmask & (ULbits - (rvbits & kmask));
- if (hi0bits(rvb->x[rvb->wds - 1]) != j)
- rvbits++;
- irv = STRTOG_Normal | STRTOG_Inexhi;
- }
- else {
- if (bbbits == 1)
- goto undfl;
- decrement(rvb);
- irv = STRTOG_Normal | STRTOG_Inexlo;
- }
- break;
- }
- if ((dval(&adj) = ratio(delta, bs)) <= 2.) {
- adj1:
- inex = STRTOG_Inexlo;
- if (dsign) {
- asub = 0;
- inex = STRTOG_Inexhi;
- }
- else if (denorm && bbbits <= 1) {
- undfl:
- rvb->wds = 0;
- rve = emin;
- irv = STRTOG_Underflow | STRTOG_Inexlo;
- break;
- }
- adj0 = dval(&adj) = 1.;
- }
- else {
- adj0 = dval(&adj) *= 0.5;
- if (dsign) {
- asub = 0;
- inex = STRTOG_Inexlo;
- }
- if (dval(&adj) < 2147483647.) {
- L = adj0;
- adj0 -= L;
- switch(rd) {
- case 0:
- if (adj0 >= .5)
- goto inc_L;
- break;
- case 1:
- if (asub && adj0 > 0.)
- goto inc_L;
- break;
- case 2:
- if (!asub && adj0 > 0.) {
- inc_L:
- L++;
- inex = STRTOG_Inexact - inex;
- }
- }
- dval(&adj) = L;
- }
- }
- y = rve + rvbits;
-
- /* adj *= ulp(dval(&rv)); */
- /* if (asub) rv -= adj; else rv += adj; */
-
- if (!denorm && rvbits < nbits) {
- rvb = lshift(rvb, j = nbits - rvbits);
- rve -= j;
- rvbits = nbits;
- }
- ab = d2b(dval(&adj), &abe, &abits);
- if (abe < 0)
- rshift(ab, -abe);
- else if (abe > 0)
- ab = lshift(ab, abe);
- rvb0 = rvb;
- if (asub) {
- /* rv -= adj; */
- j = hi0bits(rvb->x[rvb->wds-1]);
- rvb = diff(rvb, ab);
- k = rvb0->wds - 1;
- if (denorm)
- /* do nothing */;
- else if (rvb->wds <= k
- || hi0bits( rvb->x[k]) >
- hi0bits(rvb0->x[k])) {
- /* unlikely; can only have lost 1 high bit */
- if (rve1 == emin) {
- --rvbits;
- denorm = 1;
- }
- else {
- rvb = lshift(rvb, 1);
- --rve;
- --rve1;
- L = finished = 0;
- }
- }
- }
- else {
- rvb = sum(rvb, ab);
- k = rvb->wds - 1;
- if (k >= rvb0->wds
- || hi0bits(rvb->x[k]) < hi0bits(rvb0->x[k])) {
- if (denorm) {
- if (++rvbits == nbits)
- denorm = 0;
- }
- else {
- rshift(rvb, 1);
- rve++;
- rve1++;
- L = 0;
- }
- }
- }
- Bfree(ab);
- Bfree(rvb0);
- if (finished)
- break;
-
- z = rve + rvbits;
- if (y == z && L) {
- /* Can we stop now? */
- tol = dval(&adj) * 5e-16; /* > max rel error */
- dval(&adj) = adj0 - .5;
- if (dval(&adj) < -tol) {
- if (adj0 > tol) {
- irv |= inex;
- break;
- }
- }
- else if (dval(&adj) > tol && adj0 < 1. - tol) {
- irv |= inex;
- break;
- }
- }
- bb0 = denorm ? 0 : trailz(rvb);
- Bfree(bb);
- Bfree(bd);
- Bfree(bs);
- Bfree(delta);
- }
- if (!denorm && (j = nbits - rvbits)) {
- if (j > 0)
- rvb = lshift(rvb, j);
- else
- rshift(rvb, -j);
- rve -= j;
- }
- *exp = rve;
- Bfree(bb);
- Bfree(bd);
- Bfree(bs);
- Bfree(bd0);
- Bfree(delta);
- if (rve > fpi->emax) {
- switch(fpi->rounding & 3) {
- case FPI_Round_near:
- goto huge;
- case FPI_Round_up:
- if (!sign)
- goto huge;
- break;
- case FPI_Round_down:
- if (sign)
- goto huge;
- }
- /* Round to largest representable magnitude */
- Bfree(rvb);
- rvb = 0;
- irv = STRTOG_Normal | STRTOG_Inexlo;
- *exp = fpi->emax;
- b = bits;
- be = b + ((fpi->nbits + 31) >> 5);
- while(b < be)
- *b++ = -1;
- if ((j = fpi->nbits & 0x1f))
- *--be >>= (32 - j);
- goto ret;
- huge:
- rvb->wds = 0;
- irv = STRTOG_Infinite | STRTOG_Overflow | STRTOG_Inexhi;
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- infnanexp:
- *exp = fpi->emax + 1;
- }
- ret:
- if (denorm) {
- if (sudden_underflow) {
- rvb->wds = 0;
- irv = STRTOG_Underflow | STRTOG_Inexlo;
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- }
- else {
- irv = (irv & ~STRTOG_Retmask) |
- (rvb->wds > 0 ? STRTOG_Denormal : STRTOG_Zero);
- if (irv & STRTOG_Inexact) {
- irv |= STRTOG_Underflow;
-#ifndef NO_ERRNO
- errno = ERANGE;
-#endif
- }
- }
- }
- if (se)
- *se = (char *)s;
- if (sign)
- irv |= STRTOG_Neg;
- if (rvb) {
- copybits(bits, nbits, rvb);
- Bfree(rvb);
- }
- return irv;
- }
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998, 2000 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
-#undef _0
-#undef _1
-
-/* one or the other of IEEE_MC68k or IEEE_8087 should be #defined */
-
-#ifdef IEEE_MC68k
-#define _0 0
-#define _1 1
-#define _2 2
-#define _3 3
-#endif
-#ifdef IEEE_8087
-#define _0 3
-#define _1 2
-#define _2 1
-#define _3 0
-#endif
-
-__float128
-strtoflt128(CONST char *s, char **sp)
-{
- static FPI fpi0 = { 113, 1-16383-113+1, 32766 - 16383 - 113 + 1, 1, SI };
- ULong bits[4];
- Long exp;
- int k;
- union { __float128 f; ULong L[4]; } u;
- ULong *L = &u.L[0];
-#ifdef Honor_FLT_ROUNDS
-#include "gdtoa_fltrnds.h"
-#else
-#define fpi &fpi0
-#endif
-
- k = strtodg(s, sp, fpi, &exp, bits);
- switch(k & STRTOG_Retmask) {
- case STRTOG_NoNumber:
- case STRTOG_Zero:
- L[0] = L[1] = L[2] = L[3] = 0;
- break;
-
- case STRTOG_Normal:
- case STRTOG_NaNbits:
- L[_3] = bits[0];
- L[_2] = bits[1];
- L[_1] = bits[2];
- L[_0] = (bits[3] & ~0x10000) | ((exp + 0x3fff + 112) << 16);
- break;
-
- case STRTOG_Denormal:
- L[_3] = bits[0];
- L[_2] = bits[1];
- L[_1] = bits[2];
- L[_0] = bits[3];
- break;
-
- case STRTOG_Infinite:
- L[_0] = 0x7fff0000;
- L[_1] = L[_2] = L[_3] = 0;
- break;
-
- case STRTOG_NaN:
- L[_0] = ld_QNAN3;
- L[_1] = ld_QNAN2;
- L[_2] = ld_QNAN1;
- L[_3] = ld_QNAN0;
- }
- if (k & STRTOG_Neg)
- L[_0] |= 0x80000000L;
- return u.f;
-}
+++ /dev/null
-/****************************************************************
-
-The author of this software is David M. Gay.
-
-Copyright (C) 1998 by Lucent Technologies
-All Rights Reserved
-
-Permission to use, copy, modify, and distribute this software and
-its documentation for any purpose and without fee is hereby
-granted, provided that the above copyright notice appear in all
-copies and that both that the copyright notice and this
-permission notice and warranty disclaimer appear in supporting
-documentation, and that the name of Lucent or any of its entities
-not be used in advertising or publicity pertaining to
-distribution of the software without specific, written prior
-permission.
-
-LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
-INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
-IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
-SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
-WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
-IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
-THIS SOFTWARE.
-
-****************************************************************/
-
-/* Please send bug reports to David M. Gay (dmg at acm dot org,
- * with " at " changed at "@" and " dot " changed to "."). */
-
-#include "gdtoaimp.h"
-
- Bigint *
-#ifdef KR_headers
-sum(a, b) Bigint *a; Bigint *b;
-#else
-sum(Bigint *a, Bigint *b)
-#endif
-{
- Bigint *c;
- ULong carry, *xc, *xa, *xb, *xe, y;
-#ifdef Pack_32
- ULong z;
-#endif
-
- if (a->wds < b->wds) {
- c = b; b = a; a = c;
- }
- c = Balloc(a->k);
- c->wds = a->wds;
- carry = 0;
- xa = a->x;
- xb = b->x;
- xc = c->x;
- xe = xc + b->wds;
-#ifdef Pack_32
- do {
- y = (*xa & 0xffff) + (*xb & 0xffff) + carry;
- carry = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) + (*xb++ >> 16) + carry;
- carry = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
- while(xc < xe);
- xe += a->wds - b->wds;
- while(xc < xe) {
- y = (*xa & 0xffff) + carry;
- carry = (y & 0x10000) >> 16;
- z = (*xa++ >> 16) + carry;
- carry = (z & 0x10000) >> 16;
- Storeinc(xc, z, y);
- }
-#else
- do {
- y = *xa++ + *xb++ + carry;
- carry = (y & 0x10000) >> 16;
- *xc++ = y & 0xffff;
- }
- while(xc < xe);
- xe += a->wds - b->wds;
- while(xc < xe) {
- y = *xa++ + carry;
- carry = (y & 0x10000) >> 16;
- *xc++ = y & 0xffff;
- }
-#endif
- if (carry) {
- if (c->wds == c->maxwds) {
- b = Balloc(c->k + 1);
- Bcopy(b, c);
- Bfree(c);
- c = b;
- }
- c->x[c->wds++] = 1;
- }
- return c;
- }
@item @code{FLT128_MANT_DIG}: number of digits in the mantissa (bit precision)
@item @code{FLT128_MIN_EXP}: maximal negative exponent
@item @code{FLT128_MAX_EXP}: maximal positive exponent
+@item @code{FLT128_DIG}: number of decimal digits in the mantissa
+@item @code{FLT128_MIN_10_EXP}: maximal negative decimal exponent
+@item @code{FLT128_MAX_10_EXP}: maximal positive decimal exponent
@end table
The following mathematical constants of type @code{__float128} are defined.
@item Syntax
@code{__float128 strtoflt128 (const char *s, char **sp)}
-@c The return values are defined in gdtoa/gdtoa.h STRTOG_*
-@c However, the values are currently not exported - thus we
-@c do not define them here, either.
-
@item @emph{Arguments}:
@multitable @columnfractions .15 .70
@item @var{s} @tab input string
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
mp_limb_t
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
mp_limb_t
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Compare OP1_PTR/OP1_SIZE with OP2_PTR/OP2_SIZE.
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA. */
+#include <config.h>
#include "gmp-impl.h" /* This defines BITS_PER_MP_LIMB. */
#include "fpioconst.h"
mp_size_t mpn_extract_flt128 (mp_ptr res_ptr, mp_size_t size, int *expt,
int *is_neg, __float128 value) attribute_hidden;
+#define mpn_construct_float128 __MPN(construct_float128)
+__float128 mpn_construct_float128 (mp_srcptr frac_ptr, int expt, int sign)
+ attribute_hidden;
+
#define mpn_divmod(qp,np,nsize,dp,dsize) mpn_divrem (qp,0,np,nsize,dp,dsize)
static inline mp_limb_t
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE digits long) CNT bits to the left
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
mp_limb_t
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
/* Won't work for EBCDIC. */
#undef isupper
#undef isdigit
+#undef isxdigit
#undef tolower
-#define isupper(x) ((x) >= 'A' && (x) <= 'Z')
-#define isdigit(x) ((x) >= '0' && (x) <= '9')
-#define tolower(x) (isupper (x) ? (x) - 'A' + 'a' : (x))
+#define isupper(x) \
+ ({__typeof(x) __is_x = (x); __is_x >= 'A' && __is_x <= 'Z'; })
+#define isdigit(x) \
+ ({__typeof(x) __is_x = (x); __is_x >= '0' && __is_x <= '9'; })
+#define isxdigit(x) \
+ ({__typeof(x) __is_x = (x); \
+ (__is_x >= '0' && __is_x <= '9') \
+ || ((x) >= 'A' && (x) <= 'F') \
+ || ((x) >= 'a' && (x) <= 'f'); })
+#define tolower(x) \
+ ({__typeof(x) __is_x = (x); \
+ (__is_x >= 'A' && __is_x <= 'Z') ? __is_x - 'A' + 'a' : __is_x; })
#endif
#ifndef CHAR_MAX
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
/* Shift U (pointed to by UP and USIZE limbs long) CNT bits to the right
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
mp_limb_t
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */
+#include <config.h>
#include "gmp-impl.h"
mp_limb_t
#define FLT128_MANT_DIG 113
#define FLT128_MIN_EXP (-16381)
#define FLT128_MAX_EXP 16384
-/* TODO -- One day, we need to add the following macros:
- FLT128_DIG, FLT128_MIN_10_EXP, FLT128_MAX_10_EXP */
+#define FLT128_DIG 33
+#define FLT128_MIN_10_EXP (-4931)
+#define FLT128_MAX_10_EXP 4932
#define HUGE_VALQ __builtin_huge_valq()
--- /dev/null
+/* Internal header for proving correct grouping in strings of numbers.
+ Copyright (C) 1995,1996,1997,1998,2000,2003 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@gnu.ai.mit.edu>, 1995.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+/* Find the maximum prefix of the string between BEGIN and END which
+ satisfies the grouping rules. It is assumed that at least one digit
+ follows BEGIN directly. */
+extern const wchar_t *__correctly_grouped_prefixwc (const wchar_t *begin,
+ const wchar_t *end,
+ wchar_t thousands,
+ const char *grouping)
+ attribute_hidden;
+
+extern const char *__correctly_grouped_prefixmb (const char *begin,
+ const char *end,
+ const char *thousands,
+ const char *grouping)
+ attribute_hidden;
+
+/* Disable grouping support for now. */
+#define __correctly_grouped_prefixmb(b,e,t,g) e
--- /dev/null
+/* Copyright (C) 1995,1996,1997,1998,1999,2002,2003
+ Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+#include <config.h>
+#include <float.h>
+#include <math.h>
+#include "../printf/gmp-impl.h"
+
+/* Convert a multi-precision integer of the needed number of bits (113 for
+ long double) and an integral power of two to a `long double' in IEEE854
+ quad-precision format. */
+
+__float128
+mpn_construct_float128 (mp_srcptr frac_ptr, int expt, int sign)
+{
+ ieee854_float128 u;
+
+ u.ieee.negative = sign;
+ u.ieee.exponent = expt + IEEE854_FLOAT128_BIAS;
+#if BITS_PER_MP_LIMB == 32
+ u.ieee.mant_low = (((uint64_t) frac_ptr[1]) << 32)
+ | (frac_ptr[0] & 0xffffffff);
+ u.ieee.mant_high = (((uint64_t) frac_ptr[3]
+ & (((mp_limb_t) 1 << (FLT128_MANT_DIG - 96)) - 1))
+ << 32) | (frac_ptr[2] & 0xffffffff);
+#elif BITS_PER_MP_LIMB == 64
+ u.ieee.mant_low = frac_ptr[0];
+ u.ieee.mant_high = frac_ptr[1]
+ & (((mp_limb_t) 1 << (FLT128_MANT_DIG - 64)) - 1);
+#else
+ #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+ return u.value;
+}
--- /dev/null
+/* Convert string representing a number to float value, using given locale.
+ Copyright (C) 1997,1998,2002,2004,2005,2006,2007,2008,2009,2010
+ Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+ Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+#include <config.h>
+#include <stdarg.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#define NDEBUG 1
+#include <assert.h>
+#ifdef HAVE_ERRNO_H
+#include <errno.h>
+#endif
+#include "../printf/quadmath-printf.h"
+#include "../printf/fpioconst.h"
+
+
+#undef L_
+#ifdef USE_WIDE_CHAR
+# define STRING_TYPE wchar_t
+# define CHAR_TYPE wint_t
+# define L_(Ch) L##Ch
+# define ISSPACE(Ch) __iswspace_l ((Ch), loc)
+# define ISDIGIT(Ch) __iswdigit_l ((Ch), loc)
+# define ISXDIGIT(Ch) __iswxdigit_l ((Ch), loc)
+# define TOLOWER(Ch) __towlower_l ((Ch), loc)
+# define TOLOWER_C(Ch) __towlower_l ((Ch), _nl_C_locobj_ptr)
+# define STRNCASECMP(S1, S2, N) \
+ __wcsncasecmp_l ((S1), (S2), (N), _nl_C_locobj_ptr)
+# define STRTOULL(S, E, B) ____wcstoull_l_internal ((S), (E), (B), 0, loc)
+#else
+# define STRING_TYPE char
+# define CHAR_TYPE char
+# define L_(Ch) Ch
+# define ISSPACE(Ch) isspace (Ch)
+# define ISDIGIT(Ch) isdigit (Ch)
+# define ISXDIGIT(Ch) isxdigit (Ch)
+# define TOLOWER(Ch) tolower (Ch)
+# define TOLOWER_C(Ch) \
+ ({__typeof(Ch) __tlc = (Ch); \
+ (__tlc >= 'A' && __tlc <= 'Z') ? __tlc - 'A' + 'a' : __tlc; })
+# define STRNCASECMP(S1, S2, N) \
+ __quadmath_strncasecmp_c (S1, S2, N)
+# ifdef HAVE_STRTOULL
+# define STRTOULL(S, E, B) strtoull (S, E, B)
+# else
+# define STRTOULL(S, E, B) strtoul (S, E, B)
+# endif
+
+static inline int
+__quadmath_strncasecmp_c (const char *s1, const char *s2, size_t n)
+{
+ const unsigned char *p1 = (const unsigned char *) s1;
+ const unsigned char *p2 = (const unsigned char *) s2;
+ int result;
+ if (p1 == p2 || n == 0)
+ return 0;
+ while ((result = TOLOWER_C (*p1) - TOLOWER_C (*p2++)) == 0)
+ if (*p1++ == '\0' || --n == 0)
+ break;
+
+ return result;
+}
+#endif
+
+
+/* Constants we need from float.h; select the set for the FLOAT precision. */
+#define MANT_DIG PASTE(FLT,_MANT_DIG)
+#define DIG PASTE(FLT,_DIG)
+#define MAX_EXP PASTE(FLT,_MAX_EXP)
+#define MIN_EXP PASTE(FLT,_MIN_EXP)
+#define MAX_10_EXP PASTE(FLT,_MAX_10_EXP)
+#define MIN_10_EXP PASTE(FLT,_MIN_10_EXP)
+
+/* Extra macros required to get FLT expanded before the pasting. */
+#define PASTE(a,b) PASTE1(a,b)
+#define PASTE1(a,b) a##b
+
+/* Function to construct a floating point number from an MP integer
+ containing the fraction bits, a base 2 exponent, and a sign flag. */
+extern FLOAT MPN2FLOAT (mp_srcptr mpn, int exponent, int negative);
+\f
+/* Definitions according to limb size used. */
+#if BITS_PER_MP_LIMB == 32
+# define MAX_DIG_PER_LIMB 9
+# define MAX_FAC_PER_LIMB 1000000000UL
+#elif BITS_PER_MP_LIMB == 64
+# define MAX_DIG_PER_LIMB 19
+# define MAX_FAC_PER_LIMB 10000000000000000000ULL
+#else
+# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+#define _tens_in_limb __quadmath_tens_in_limb
+extern const mp_limb_t _tens_in_limb[MAX_DIG_PER_LIMB + 1] attribute_hidden;
+\f
+#ifndef howmany
+#define howmany(x,y) (((x)+((y)-1))/(y))
+#endif
+#define SWAP(x, y) ({ typeof(x) _tmp = x; x = y; y = _tmp; })
+
+#define NDIG (MAX_10_EXP - MIN_10_EXP + 2 * MANT_DIG)
+#define HEXNDIG ((MAX_EXP - MIN_EXP + 7) / 8 + 2 * MANT_DIG)
+#define RETURN_LIMB_SIZE howmany (MANT_DIG, BITS_PER_MP_LIMB)
+
+#define RETURN(val,end) \
+ do { if (endptr != NULL) *endptr = (STRING_TYPE *) (end); \
+ return val; } while (0)
+
+/* Maximum size necessary for mpn integers to hold floating point numbers. */
+#define MPNSIZE (howmany (MAX_EXP + 2 * MANT_DIG, BITS_PER_MP_LIMB) \
+ + 2)
+/* Declare an mpn integer variable that big. */
+#define MPN_VAR(name) mp_limb_t name[MPNSIZE]; mp_size_t name##size
+/* Copy an mpn integer value. */
+#define MPN_ASSIGN(dst, src) \
+ memcpy (dst, src, (dst##size = src##size) * sizeof (mp_limb_t))
+
+
+/* Return a floating point number of the needed type according to the given
+ multi-precision number after possible rounding. */
+static FLOAT
+round_and_return (mp_limb_t *retval, int exponent, int negative,
+ mp_limb_t round_limb, mp_size_t round_bit, int more_bits)
+{
+ if (exponent < MIN_EXP - 1)
+ {
+ mp_size_t shift = MIN_EXP - 1 - exponent;
+
+ if (shift > MANT_DIG)
+ {
+#if defined HAVE_ERRNO_H && defined EDOM
+ errno = EDOM;
+#endif
+ return 0.0;
+ }
+
+ more_bits |= (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0;
+ if (shift == MANT_DIG)
+ /* This is a special case to handle the very seldom case where
+ the mantissa will be empty after the shift. */
+ {
+ int i;
+
+ round_limb = retval[RETURN_LIMB_SIZE - 1];
+ round_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
+ for (i = 0; i < RETURN_LIMB_SIZE; ++i)
+ more_bits |= retval[i] != 0;
+ MPN_ZERO (retval, RETURN_LIMB_SIZE);
+ }
+ else if (shift >= BITS_PER_MP_LIMB)
+ {
+ int i;
+
+ round_limb = retval[(shift - 1) / BITS_PER_MP_LIMB];
+ round_bit = (shift - 1) % BITS_PER_MP_LIMB;
+ for (i = 0; i < (shift - 1) / BITS_PER_MP_LIMB; ++i)
+ more_bits |= retval[i] != 0;
+ more_bits |= ((round_limb & ((((mp_limb_t) 1) << round_bit) - 1))
+ != 0);
+
+ (void) mpn_rshift (retval, &retval[shift / BITS_PER_MP_LIMB],
+ RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB),
+ shift % BITS_PER_MP_LIMB);
+ MPN_ZERO (&retval[RETURN_LIMB_SIZE - (shift / BITS_PER_MP_LIMB)],
+ shift / BITS_PER_MP_LIMB);
+ }
+ else if (shift > 0)
+ {
+ round_limb = retval[0];
+ round_bit = shift - 1;
+ (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, shift);
+ }
+ /* This is a hook for the m68k long double format, where the
+ exponent bias is the same for normalized and denormalized
+ numbers. */
+#ifndef DENORM_EXP
+# define DENORM_EXP (MIN_EXP - 2)
+#endif
+ exponent = DENORM_EXP;
+#if defined HAVE_ERRNO_H && defined ERANGE
+ errno = ERANGE;
+#endif
+ }
+
+ if ((round_limb & (((mp_limb_t) 1) << round_bit)) != 0
+ && (more_bits || (retval[0] & 1) != 0
+ || (round_limb & ((((mp_limb_t) 1) << round_bit) - 1)) != 0))
+ {
+ mp_limb_t cy = mpn_add_1 (retval, retval, RETURN_LIMB_SIZE, 1);
+
+ if (((MANT_DIG % BITS_PER_MP_LIMB) == 0 && cy) ||
+ ((MANT_DIG % BITS_PER_MP_LIMB) != 0 &&
+ (retval[RETURN_LIMB_SIZE - 1]
+ & (((mp_limb_t) 1) << (MANT_DIG % BITS_PER_MP_LIMB))) != 0))
+ {
+ ++exponent;
+ (void) mpn_rshift (retval, retval, RETURN_LIMB_SIZE, 1);
+ retval[RETURN_LIMB_SIZE - 1]
+ |= ((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB);
+ }
+ else if (exponent == DENORM_EXP
+ && (retval[RETURN_LIMB_SIZE - 1]
+ & (((mp_limb_t) 1) << ((MANT_DIG - 1) % BITS_PER_MP_LIMB)))
+ != 0)
+ /* The number was denormalized but now normalized. */
+ exponent = MIN_EXP - 1;
+ }
+
+ if (exponent > MAX_EXP)
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+
+ return MPN2FLOAT (retval, exponent, negative);
+}
+
+
+/* Read a multi-precision integer starting at STR with exactly DIGCNT digits
+ into N. Return the size of the number limbs in NSIZE at the first
+ character od the string that is not part of the integer as the function
+ value. If the EXPONENT is small enough to be taken as an additional
+ factor for the resulting number (see code) multiply by it. */
+static const STRING_TYPE *
+str_to_mpn (const STRING_TYPE *str, int digcnt, mp_limb_t *n, mp_size_t *nsize,
+ int *exponent
+#ifndef USE_WIDE_CHAR
+ , const char *decimal, size_t decimal_len, const char *thousands
+#endif
+
+ )
+{
+ /* Number of digits for actual limb. */
+ int cnt = 0;
+ mp_limb_t low = 0;
+ mp_limb_t start;
+
+ *nsize = 0;
+ assert (digcnt > 0);
+ do
+ {
+ if (cnt == MAX_DIG_PER_LIMB)
+ {
+ if (*nsize == 0)
+ {
+ n[0] = low;
+ *nsize = 1;
+ }
+ else
+ {
+ mp_limb_t cy;
+ cy = mpn_mul_1 (n, n, *nsize, MAX_FAC_PER_LIMB);
+ cy += mpn_add_1 (n, n, *nsize, low);
+ if (cy != 0)
+ {
+ n[*nsize] = cy;
+ ++(*nsize);
+ }
+ }
+ cnt = 0;
+ low = 0;
+ }
+
+ /* There might be thousands separators or radix characters in
+ the string. But these all can be ignored because we know the
+ format of the number is correct and we have an exact number
+ of characters to read. */
+#ifdef USE_WIDE_CHAR
+ if (*str < L'0' || *str > L'9')
+ ++str;
+#else
+ if (*str < '0' || *str > '9')
+ {
+ int inner = 0;
+ if (thousands != NULL && *str == *thousands
+ && ({ for (inner = 1; thousands[inner] != '\0'; ++inner)
+ if (thousands[inner] != str[inner])
+ break;
+ thousands[inner] == '\0'; }))
+ str += inner;
+ else
+ str += decimal_len;
+ }
+#endif
+ low = low * 10 + *str++ - L_('0');
+ ++cnt;
+ }
+ while (--digcnt > 0);
+
+ if (*exponent > 0 && cnt + *exponent <= MAX_DIG_PER_LIMB)
+ {
+ low *= _tens_in_limb[*exponent];
+ start = _tens_in_limb[cnt + *exponent];
+ *exponent = 0;
+ }
+ else
+ start = _tens_in_limb[cnt];
+
+ if (*nsize == 0)
+ {
+ n[0] = low;
+ *nsize = 1;
+ }
+ else
+ {
+ mp_limb_t cy;
+ cy = mpn_mul_1 (n, n, *nsize, start);
+ cy += mpn_add_1 (n, n, *nsize, low);
+ if (cy != 0)
+ n[(*nsize)++] = cy;
+ }
+
+ return str;
+}
+
+
+/* Shift {PTR, SIZE} COUNT bits to the left, and fill the vacated bits
+ with the COUNT most significant bits of LIMB.
+
+ Tege doesn't like this function so I have to write it here myself. :)
+ --drepper */
+static inline void
+__attribute ((always_inline))
+mpn_lshift_1 (mp_limb_t *ptr, mp_size_t size, unsigned int count,
+ mp_limb_t limb)
+{
+ if (__builtin_constant_p (count) && count == BITS_PER_MP_LIMB)
+ {
+ /* Optimize the case of shifting by exactly a word:
+ just copy words, with no actual bit-shifting. */
+ mp_size_t i;
+ for (i = size - 1; i > 0; --i)
+ ptr[i] = ptr[i - 1];
+ ptr[0] = limb;
+ }
+ else
+ {
+ (void) mpn_lshift (ptr, ptr, size, count);
+ ptr[0] |= limb >> (BITS_PER_MP_LIMB - count);
+ }
+}
+
+
+#define INTERNAL(x) INTERNAL1(x)
+#define INTERNAL1(x) __##x##_internal
+#ifndef ____STRTOF_INTERNAL
+# define ____STRTOF_INTERNAL INTERNAL (__STRTOF)
+#endif
+
+/* This file defines a function to check for correct grouping. */
+#include "grouping.h"
+
+
+/* Return a floating point number with the value of the given string NPTR.
+ Set *ENDPTR to the character after the last used one. If the number is
+ smaller than the smallest representable number, set `errno' to ERANGE and
+ return 0.0. If the number is too big to be represented, set `errno' to
+ ERANGE and return HUGE_VAL with the appropriate sign. */
+FLOAT
+____STRTOF_INTERNAL (nptr, endptr, group)
+ const STRING_TYPE *nptr;
+ STRING_TYPE **endptr;
+ int group;
+{
+ int negative; /* The sign of the number. */
+ MPN_VAR (num); /* MP representation of the number. */
+ int exponent; /* Exponent of the number. */
+
+ /* Numbers starting `0X' or `0x' have to be processed with base 16. */
+ int base = 10;
+
+ /* When we have to compute fractional digits we form a fraction with a
+ second multi-precision number (and we sometimes need a second for
+ temporary results). */
+ MPN_VAR (den);
+
+ /* Representation for the return value. */
+ mp_limb_t retval[RETURN_LIMB_SIZE];
+ /* Number of bits currently in result value. */
+ int bits;
+
+ /* Running pointer after the last character processed in the string. */
+ const STRING_TYPE *cp, *tp;
+ /* Start of significant part of the number. */
+ const STRING_TYPE *startp, *start_of_digits;
+ /* Points at the character following the integer and fractional digits. */
+ const STRING_TYPE *expp;
+ /* Total number of digit and number of digits in integer part. */
+ int dig_no, int_no, lead_zero;
+ /* Contains the last character read. */
+ CHAR_TYPE c;
+
+ /* The radix character of the current locale. */
+#ifdef USE_WIDE_CHAR
+ wchar_t decimal;
+#else
+ const char *decimal;
+ size_t decimal_len;
+#endif
+ /* The thousands character of the current locale. */
+#ifdef USE_WIDE_CHAR
+ wchar_t thousands = L'\0';
+#else
+ const char *thousands = NULL;
+#endif
+ /* The numeric grouping specification of the current locale,
+ in the format described in <locale.h>. */
+ const char *grouping;
+ /* Used in several places. */
+ int cnt;
+
+#if defined USE_LOCALECONV && !defined USE_NL_LANGINFO
+ const struct lconv *lc = localeconv ();
+#endif
+
+ if (__builtin_expect (group, 0))
+ {
+#ifdef USE_NL_LANGINFO
+ grouping = nl_langinfo (GROUPING);
+ if (*grouping <= 0 || *grouping == CHAR_MAX)
+ grouping = NULL;
+ else
+ {
+ /* Figure out the thousands separator character. */
+#ifdef USE_WIDE_CHAR
+ thousands = nl_langinfo_wc (_NL_NUMERIC_THOUSANDS_SEP_WC);
+ if (thousands == L'\0')
+ grouping = NULL;
+#else
+ thousands = nl_langinfo (THOUSANDS_SEP);
+ if (*thousands == '\0')
+ {
+ thousands = NULL;
+ grouping = NULL;
+ }
+#endif
+ }
+#elif defined USE_LOCALECONV
+ grouping = lc->grouping;
+ if (grouping == NULL || *grouping <= 0 || *grouping == CHAR_MAX)
+ grouping = NULL;
+ else
+ {
+ /* Figure out the thousands separator character. */
+ thousands = lc->thousands_sep;
+ if (thousands == NULL || *thousands == '\0')
+ {
+ thousands = NULL;
+ grouping = NULL;
+ }
+ }
+#else
+ grouping = NULL;
+#endif
+ }
+ else
+ grouping = NULL;
+
+ /* Find the locale's decimal point character. */
+#ifdef USE_WIDE_CHAR
+ decimal = nl_langinfo_wc (_NL_NUMERIC_DECIMAL_POINT_WC);
+ assert (decimal != L'\0');
+# define decimal_len 1
+#else
+#ifdef USE_NL_LANGINFO
+ decimal = nl_langinfo (DECIMAL_POINT);
+ decimal_len = strlen (decimal);
+ assert (decimal_len > 0);
+#elif defined USE_LOCALECONV
+ decimal = lc->decimal_point;
+ if (decimal == NULL || *decimal == '\0')
+ decimal = ".";
+ decimal_len = strlen (decimal);
+#else
+ decimal = ".";
+ decimal_len = 1;
+#endif
+#endif
+
+ /* Prepare number representation. */
+ exponent = 0;
+ negative = 0;
+ bits = 0;
+
+ /* Parse string to get maximal legal prefix. We need the number of
+ characters of the integer part, the fractional part and the exponent. */
+ cp = nptr - 1;
+ /* Ignore leading white space. */
+ do
+ c = *++cp;
+ while (ISSPACE (c));
+
+ /* Get sign of the result. */
+ if (c == L_('-'))
+ {
+ negative = 1;
+ c = *++cp;
+ }
+ else if (c == L_('+'))
+ c = *++cp;
+
+ /* Return 0.0 if no legal string is found.
+ No character is used even if a sign was found. */
+#ifdef USE_WIDE_CHAR
+ if (c == (wint_t) decimal
+ && (wint_t) cp[1] >= L'0' && (wint_t) cp[1] <= L'9')
+ {
+ /* We accept it. This funny construct is here only to indent
+ the code correctly. */
+ }
+#else
+ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
+ if (cp[cnt] != decimal[cnt])
+ break;
+ if (decimal[cnt] == '\0' && cp[cnt] >= '0' && cp[cnt] <= '9')
+ {
+ /* We accept it. This funny construct is here only to indent
+ the code correctly. */
+ }
+#endif
+ else if (c < L_('0') || c > L_('9'))
+ {
+ /* Check for `INF' or `INFINITY'. */
+ CHAR_TYPE lowc = TOLOWER_C (c);
+
+ if (lowc == L_('i') && STRNCASECMP (cp, L_("inf"), 3) == 0)
+ {
+ /* Return +/- infinity. */
+ if (endptr != NULL)
+ *endptr = (STRING_TYPE *)
+ (cp + (STRNCASECMP (cp + 3, L_("inity"), 5) == 0
+ ? 8 : 3));
+
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+ }
+
+ if (lowc == L_('n') && STRNCASECMP (cp, L_("nan"), 3) == 0)
+ {
+ /* Return NaN. */
+ FLOAT retval = NAN;
+
+ cp += 3;
+
+ /* Match `(n-char-sequence-digit)'. */
+ if (*cp == L_('('))
+ {
+ const STRING_TYPE *startp = cp;
+ do
+ ++cp;
+ while ((*cp >= L_('0') && *cp <= L_('9'))
+ || ({ CHAR_TYPE lo = TOLOWER (*cp);
+ lo >= L_('a') && lo <= L_('z'); })
+ || *cp == L_('_'));
+
+ if (*cp != L_(')'))
+ /* The closing brace is missing. Only match the NAN
+ part. */
+ cp = startp;
+ else
+ {
+ /* This is a system-dependent way to specify the
+ bitmask used for the NaN. We expect it to be
+ a number which is put in the mantissa of the
+ number. */
+ STRING_TYPE *endp;
+ unsigned long long int mant;
+
+ mant = STRTOULL (startp + 1, &endp, 0);
+ if (endp == cp)
+ SET_MANTISSA (retval, mant);
+
+ /* Consume the closing brace. */
+ ++cp;
+ }
+ }
+
+ if (endptr != NULL)
+ *endptr = (STRING_TYPE *) cp;
+
+ return retval;
+ }
+
+ /* It is really a text we do not recognize. */
+ RETURN (0.0, nptr);
+ }
+
+ /* First look whether we are faced with a hexadecimal number. */
+ if (c == L_('0') && TOLOWER (cp[1]) == L_('x'))
+ {
+ /* Okay, it is a hexa-decimal number. Remember this and skip
+ the characters. BTW: hexadecimal numbers must not be
+ grouped. */
+ base = 16;
+ cp += 2;
+ c = *cp;
+ grouping = NULL;
+ }
+
+ /* Record the start of the digits, in case we will check their grouping. */
+ start_of_digits = startp = cp;
+
+ /* Ignore leading zeroes. This helps us to avoid useless computations. */
+#ifdef USE_WIDE_CHAR
+ while (c == L'0' || ((wint_t) thousands != L'\0' && c == (wint_t) thousands))
+ c = *++cp;
+#else
+ if (__builtin_expect (thousands == NULL, 1))
+ while (c == '0')
+ c = *++cp;
+ else
+ {
+ /* We also have the multibyte thousands string. */
+ while (1)
+ {
+ if (c != '0')
+ {
+ for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
+ if (thousands[cnt] != cp[cnt])
+ break;
+ if (thousands[cnt] != '\0')
+ break;
+ cp += cnt - 1;
+ }
+ c = *++cp;
+ }
+ }
+#endif
+
+ /* If no other digit but a '0' is found the result is 0.0.
+ Return current read pointer. */
+ CHAR_TYPE lowc = TOLOWER (c);
+ if (!((c >= L_('0') && c <= L_('9'))
+ || (base == 16 && lowc >= L_('a') && lowc <= L_('f'))
+ || (
+#ifdef USE_WIDE_CHAR
+ c == (wint_t) decimal
+#else
+ ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
+ if (decimal[cnt] != cp[cnt])
+ break;
+ decimal[cnt] == '\0'; })
+#endif
+ /* '0x.' alone is not a valid hexadecimal number.
+ '.' alone is not valid either, but that has been checked
+ already earlier. */
+ && (base != 16
+ || cp != start_of_digits
+ || (cp[decimal_len] >= L_('0') && cp[decimal_len] <= L_('9'))
+ || ({ CHAR_TYPE lo = TOLOWER (cp[decimal_len]);
+ lo >= L_('a') && lo <= L_('f'); })))
+ || (base == 16 && (cp != start_of_digits
+ && lowc == L_('p')))
+ || (base != 16 && lowc == L_('e'))))
+ {
+#ifdef USE_WIDE_CHAR
+ tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
+ grouping);
+#else
+ tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
+ grouping);
+#endif
+ /* If TP is at the start of the digits, there was no correctly
+ grouped prefix of the string; so no number found. */
+ RETURN (negative ? -0.0 : 0.0,
+ tp == start_of_digits ? (base == 16 ? cp - 1 : nptr) : tp);
+ }
+
+ /* Remember first significant digit and read following characters until the
+ decimal point, exponent character or any non-FP number character. */
+ startp = cp;
+ dig_no = 0;
+ while (1)
+ {
+ if ((c >= L_('0') && c <= L_('9'))
+ || (base == 16
+ && ({ CHAR_TYPE lo = TOLOWER (c);
+ lo >= L_('a') && lo <= L_('f'); })))
+ ++dig_no;
+ else
+ {
+#ifdef USE_WIDE_CHAR
+ if (__builtin_expect ((wint_t) thousands == L'\0', 1)
+ || c != (wint_t) thousands)
+ /* Not a digit or separator: end of the integer part. */
+ break;
+#else
+ if (__builtin_expect (thousands == NULL, 1))
+ break;
+ else
+ {
+ for (cnt = 0; thousands[cnt] != '\0'; ++cnt)
+ if (thousands[cnt] != cp[cnt])
+ break;
+ if (thousands[cnt] != '\0')
+ break;
+ cp += cnt - 1;
+ }
+#endif
+ }
+ c = *++cp;
+ }
+
+ if (__builtin_expect (grouping != NULL, 0) && cp > start_of_digits)
+ {
+ /* Check the grouping of the digits. */
+#ifdef USE_WIDE_CHAR
+ tp = __correctly_grouped_prefixwc (start_of_digits, cp, thousands,
+ grouping);
+#else
+ tp = __correctly_grouped_prefixmb (start_of_digits, cp, thousands,
+ grouping);
+#endif
+ if (cp != tp)
+ {
+ /* Less than the entire string was correctly grouped. */
+
+ if (tp == start_of_digits)
+ /* No valid group of numbers at all: no valid number. */
+ RETURN (0.0, nptr);
+
+ if (tp < startp)
+ /* The number is validly grouped, but consists
+ only of zeroes. The whole value is zero. */
+ RETURN (negative ? -0.0 : 0.0, tp);
+
+ /* Recompute DIG_NO so we won't read more digits than
+ are properly grouped. */
+ cp = tp;
+ dig_no = 0;
+ for (tp = startp; tp < cp; ++tp)
+ if (*tp >= L_('0') && *tp <= L_('9'))
+ ++dig_no;
+
+ int_no = dig_no;
+ lead_zero = 0;
+
+ goto number_parsed;
+ }
+ }
+
+ /* We have the number of digits in the integer part. Whether these
+ are all or any is really a fractional digit will be decided
+ later. */
+ int_no = dig_no;
+ lead_zero = int_no == 0 ? -1 : 0;
+
+ /* Read the fractional digits. A special case are the 'american
+ style' numbers like `16.' i.e. with decimal point but without
+ trailing digits. */
+ if (
+#ifdef USE_WIDE_CHAR
+ c == (wint_t) decimal
+#else
+ ({ for (cnt = 0; decimal[cnt] != '\0'; ++cnt)
+ if (decimal[cnt] != cp[cnt])
+ break;
+ decimal[cnt] == '\0'; })
+#endif
+ )
+ {
+ cp += decimal_len;
+ c = *cp;
+ while ((c >= L_('0') && c <= L_('9')) ||
+ (base == 16 && ({ CHAR_TYPE lo = TOLOWER (c);
+ lo >= L_('a') && lo <= L_('f'); })))
+ {
+ if (c != L_('0') && lead_zero == -1)
+ lead_zero = dig_no - int_no;
+ ++dig_no;
+ c = *++cp;
+ }
+ }
+
+ /* Remember start of exponent (if any). */
+ expp = cp;
+
+ /* Read exponent. */
+ lowc = TOLOWER (c);
+ if ((base == 16 && lowc == L_('p'))
+ || (base != 16 && lowc == L_('e')))
+ {
+ int exp_negative = 0;
+
+ c = *++cp;
+ if (c == L_('-'))
+ {
+ exp_negative = 1;
+ c = *++cp;
+ }
+ else if (c == L_('+'))
+ c = *++cp;
+
+ if (c >= L_('0') && c <= L_('9'))
+ {
+ int exp_limit;
+
+ /* Get the exponent limit. */
+ if (base == 16)
+ exp_limit = (exp_negative ?
+ -MIN_EXP + MANT_DIG + 4 * int_no :
+ MAX_EXP - 4 * int_no + 4 * lead_zero + 3);
+ else
+ exp_limit = (exp_negative ?
+ -MIN_10_EXP + MANT_DIG + int_no :
+ MAX_10_EXP - int_no + lead_zero + 1);
+
+ do
+ {
+ exponent *= 10;
+ exponent += c - L_('0');
+
+ if (__builtin_expect (exponent > exp_limit, 0))
+ /* The exponent is too large/small to represent a valid
+ number. */
+ {
+ FLOAT result;
+
+ /* We have to take care for special situation: a joker
+ might have written "0.0e100000" which is in fact
+ zero. */
+ if (lead_zero == -1)
+ result = negative ? -0.0 : 0.0;
+ else
+ {
+ /* Overflow or underflow. */
+#if defined HAVE_ERRNO_H && defined ERANGE
+ errno = ERANGE;
+#endif
+ result = (exp_negative ? (negative ? -0.0 : 0.0) :
+ negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL);
+ }
+
+ /* Accept all following digits as part of the exponent. */
+ do
+ ++cp;
+ while (*cp >= L_('0') && *cp <= L_('9'));
+
+ RETURN (result, cp);
+ /* NOTREACHED */
+ }
+
+ c = *++cp;
+ }
+ while (c >= L_('0') && c <= L_('9'));
+
+ if (exp_negative)
+ exponent = -exponent;
+ }
+ else
+ cp = expp;
+ }
+
+ /* We don't want to have to work with trailing zeroes after the radix. */
+ if (dig_no > int_no)
+ {
+ while (expp[-1] == L_('0'))
+ {
+ --expp;
+ --dig_no;
+ }
+ assert (dig_no >= int_no);
+ }
+
+ if (dig_no == int_no && dig_no > 0 && exponent < 0)
+ do
+ {
+ while (! (base == 16 ? ISXDIGIT (expp[-1]) : ISDIGIT (expp[-1])))
+ --expp;
+
+ if (expp[-1] != L_('0'))
+ break;
+
+ --expp;
+ --dig_no;
+ --int_no;
+ exponent += base == 16 ? 4 : 1;
+ }
+ while (dig_no > 0 && exponent < 0);
+
+ number_parsed:
+
+ /* The whole string is parsed. Store the address of the next character. */
+ if (endptr)
+ *endptr = (STRING_TYPE *) cp;
+
+ if (dig_no == 0)
+ return negative ? -0.0 : 0.0;
+
+ if (lead_zero)
+ {
+ /* Find the decimal point */
+#ifdef USE_WIDE_CHAR
+ while (*startp != decimal)
+ ++startp;
+#else
+ while (1)
+ {
+ if (*startp == decimal[0])
+ {
+ for (cnt = 1; decimal[cnt] != '\0'; ++cnt)
+ if (decimal[cnt] != startp[cnt])
+ break;
+ if (decimal[cnt] == '\0')
+ break;
+ }
+ ++startp;
+ }
+#endif
+ startp += lead_zero + decimal_len;
+ exponent -= base == 16 ? 4 * lead_zero : lead_zero;
+ dig_no -= lead_zero;
+ }
+
+ /* If the BASE is 16 we can use a simpler algorithm. */
+ if (base == 16)
+ {
+ static const int nbits[16] = { 0, 1, 2, 2, 3, 3, 3, 3,
+ 4, 4, 4, 4, 4, 4, 4, 4 };
+ int idx = (MANT_DIG - 1) / BITS_PER_MP_LIMB;
+ int pos = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
+ mp_limb_t val;
+
+ while (!ISXDIGIT (*startp))
+ ++startp;
+ while (*startp == L_('0'))
+ ++startp;
+ if (ISDIGIT (*startp))
+ val = *startp++ - L_('0');
+ else
+ val = 10 + TOLOWER (*startp++) - L_('a');
+ bits = nbits[val];
+ /* We cannot have a leading zero. */
+ assert (bits != 0);
+
+ if (pos + 1 >= 4 || pos + 1 >= bits)
+ {
+ /* We don't have to care for wrapping. This is the normal
+ case so we add the first clause in the `if' expression as
+ an optimization. It is a compile-time constant and so does
+ not cost anything. */
+ retval[idx] = val << (pos - bits + 1);
+ pos -= bits;
+ }
+ else
+ {
+ retval[idx--] = val >> (bits - pos - 1);
+ retval[idx] = val << (BITS_PER_MP_LIMB - (bits - pos - 1));
+ pos = BITS_PER_MP_LIMB - 1 - (bits - pos - 1);
+ }
+
+ /* Adjust the exponent for the bits we are shifting in. */
+ exponent += bits - 1 + (int_no - 1) * 4;
+
+ while (--dig_no > 0 && idx >= 0)
+ {
+ if (!ISXDIGIT (*startp))
+ startp += decimal_len;
+ if (ISDIGIT (*startp))
+ val = *startp++ - L_('0');
+ else
+ val = 10 + TOLOWER (*startp++) - L_('a');
+
+ if (pos + 1 >= 4)
+ {
+ retval[idx] |= val << (pos - 4 + 1);
+ pos -= 4;
+ }
+ else
+ {
+ retval[idx--] |= val >> (4 - pos - 1);
+ val <<= BITS_PER_MP_LIMB - (4 - pos - 1);
+ if (idx < 0)
+ return round_and_return (retval, exponent, negative, val,
+ BITS_PER_MP_LIMB - 1, dig_no > 0);
+
+ retval[idx] = val;
+ pos = BITS_PER_MP_LIMB - 1 - (4 - pos - 1);
+ }
+ }
+
+ /* We ran out of digits. */
+ MPN_ZERO (retval, idx);
+
+ return round_and_return (retval, exponent, negative, 0, 0, 0);
+ }
+
+ /* Now we have the number of digits in total and the integer digits as well
+ as the exponent and its sign. We can decide whether the read digits are
+ really integer digits or belong to the fractional part; i.e. we normalize
+ 123e-2 to 1.23. */
+ {
+ register int incr = (exponent < 0 ? MAX (-int_no, exponent)
+ : MIN (dig_no - int_no, exponent));
+ int_no += incr;
+ exponent -= incr;
+ }
+
+ if (__builtin_expect (int_no + exponent > MAX_10_EXP + 1, 0))
+ {
+#if defined HAVE_ERRNO_H && defined ERANGE
+ errno = ERANGE;
+#endif
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+ }
+
+ if (__builtin_expect (exponent < MIN_10_EXP - (DIG + 1), 0))
+ {
+#if defined HAVE_ERRNO_H && defined ERANGE
+ errno = ERANGE;
+#endif
+ return negative ? -0.0 : 0.0;
+ }
+
+ if (int_no > 0)
+ {
+ /* Read the integer part as a multi-precision number to NUM. */
+ startp = str_to_mpn (startp, int_no, num, &numsize, &exponent
+#ifndef USE_WIDE_CHAR
+ , decimal, decimal_len, thousands
+#endif
+ );
+
+ if (exponent > 0)
+ {
+ /* We now multiply the gained number by the given power of ten. */
+ mp_limb_t *psrc = num;
+ mp_limb_t *pdest = den;
+ int expbit = 1;
+ const struct mp_power *ttab = &_fpioconst_pow10[0];
+
+ do
+ {
+ if ((exponent & expbit) != 0)
+ {
+ size_t size = ttab->arraysize - _FPIO_CONST_OFFSET;
+ mp_limb_t cy;
+ exponent ^= expbit;
+
+ /* FIXME: not the whole multiplication has to be
+ done. If we have the needed number of bits we
+ only need the information whether more non-zero
+ bits follow. */
+ if (numsize >= ttab->arraysize - _FPIO_CONST_OFFSET)
+ cy = mpn_mul (pdest, psrc, numsize,
+ &__tens[ttab->arrayoff
+ + _FPIO_CONST_OFFSET],
+ size);
+ else
+ cy = mpn_mul (pdest, &__tens[ttab->arrayoff
+ + _FPIO_CONST_OFFSET],
+ size, psrc, numsize);
+ numsize += size;
+ if (cy == 0)
+ --numsize;
+ (void) SWAP (psrc, pdest);
+ }
+ expbit <<= 1;
+ ++ttab;
+ }
+ while (exponent != 0);
+
+ if (psrc == den)
+ memcpy (num, den, numsize * sizeof (mp_limb_t));
+ }
+
+ /* Determine how many bits of the result we already have. */
+ count_leading_zeros (bits, num[numsize - 1]);
+ bits = numsize * BITS_PER_MP_LIMB - bits;
+
+ /* Now we know the exponent of the number in base two.
+ Check it against the maximum possible exponent. */
+ if (__builtin_expect (bits > MAX_EXP, 0))
+ {
+#if defined HAVE_ERRNO_H && defined ERANGE
+ errno = ERANGE;
+#endif
+ return negative ? -FLOAT_HUGE_VAL : FLOAT_HUGE_VAL;
+ }
+
+ /* We have already the first BITS bits of the result. Together with
+ the information whether more non-zero bits follow this is enough
+ to determine the result. */
+ if (bits > MANT_DIG)
+ {
+ int i;
+ const mp_size_t least_idx = (bits - MANT_DIG) / BITS_PER_MP_LIMB;
+ const mp_size_t least_bit = (bits - MANT_DIG) % BITS_PER_MP_LIMB;
+ const mp_size_t round_idx = least_bit == 0 ? least_idx - 1
+ : least_idx;
+ const mp_size_t round_bit = least_bit == 0 ? BITS_PER_MP_LIMB - 1
+ : least_bit - 1;
+
+ if (least_bit == 0)
+ memcpy (retval, &num[least_idx],
+ RETURN_LIMB_SIZE * sizeof (mp_limb_t));
+ else
+ {
+ for (i = least_idx; i < numsize - 1; ++i)
+ retval[i - least_idx] = (num[i] >> least_bit)
+ | (num[i + 1]
+ << (BITS_PER_MP_LIMB - least_bit));
+ if (i - least_idx < RETURN_LIMB_SIZE)
+ retval[RETURN_LIMB_SIZE - 1] = num[i] >> least_bit;
+ }
+
+ /* Check whether any limb beside the ones in RETVAL are non-zero. */
+ for (i = 0; num[i] == 0; ++i)
+ ;
+
+ return round_and_return (retval, bits - 1, negative,
+ num[round_idx], round_bit,
+ int_no < dig_no || i < round_idx);
+ /* NOTREACHED */
+ }
+ else if (dig_no == int_no)
+ {
+ const mp_size_t target_bit = (MANT_DIG - 1) % BITS_PER_MP_LIMB;
+ const mp_size_t is_bit = (bits - 1) % BITS_PER_MP_LIMB;
+
+ if (target_bit == is_bit)
+ {
+ memcpy (&retval[RETURN_LIMB_SIZE - numsize], num,
+ numsize * sizeof (mp_limb_t));
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
+ }
+ else if (target_bit > is_bit)
+ {
+ (void) mpn_lshift (&retval[RETURN_LIMB_SIZE - numsize],
+ num, numsize, target_bit - is_bit);
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize);
+ }
+ else
+ {
+ mp_limb_t cy;
+ assert (numsize < RETURN_LIMB_SIZE);
+
+ cy = mpn_rshift (&retval[RETURN_LIMB_SIZE - numsize],
+ num, numsize, is_bit - target_bit);
+ retval[RETURN_LIMB_SIZE - numsize - 1] = cy;
+ /* FIXME: the following loop can be avoided if we assume a
+ maximal MANT_DIG value. */
+ MPN_ZERO (retval, RETURN_LIMB_SIZE - numsize - 1);
+ }
+
+ return round_and_return (retval, bits - 1, negative, 0, 0, 0);
+ /* NOTREACHED */
+ }
+
+ /* Store the bits we already have. */
+ memcpy (retval, num, numsize * sizeof (mp_limb_t));
+#if RETURN_LIMB_SIZE > 1
+ if (numsize < RETURN_LIMB_SIZE)
+# if RETURN_LIMB_SIZE == 2
+ retval[numsize] = 0;
+# else
+ MPN_ZERO (retval + numsize, RETURN_LIMB_SIZE - numsize);
+# endif
+#endif
+ }
+
+ /* We have to compute at least some of the fractional digits. */
+ {
+ /* We construct a fraction and the result of the division gives us
+ the needed digits. The denominator is 1.0 multiplied by the
+ exponent of the lowest digit; i.e. 0.123 gives 123 / 1000 and
+ 123e-6 gives 123 / 1000000. */
+
+ int expbit;
+ int neg_exp;
+ int more_bits;
+ mp_limb_t cy;
+ mp_limb_t *psrc = den;
+ mp_limb_t *pdest = num;
+ const struct mp_power *ttab = &_fpioconst_pow10[0];
+
+ assert (dig_no > int_no && exponent <= 0);
+
+
+ /* For the fractional part we need not process too many digits. One
+ decimal digits gives us log_2(10) ~ 3.32 bits. If we now compute
+ ceil(BITS / 3) =: N
+ digits we should have enough bits for the result. The remaining
+ decimal digits give us the information that more bits are following.
+ This can be used while rounding. (Two added as a safety margin.) */
+ if (dig_no - int_no > (MANT_DIG - bits + 2) / 3 + 2)
+ {
+ dig_no = int_no + (MANT_DIG - bits + 2) / 3 + 2;
+ more_bits = 1;
+ }
+ else
+ more_bits = 0;
+
+ neg_exp = dig_no - int_no - exponent;
+
+ /* Construct the denominator. */
+ densize = 0;
+ expbit = 1;
+ do
+ {
+ if ((neg_exp & expbit) != 0)
+ {
+ mp_limb_t cy;
+ neg_exp ^= expbit;
+
+ if (densize == 0)
+ {
+ densize = ttab->arraysize - _FPIO_CONST_OFFSET;
+ memcpy (psrc, &__tens[ttab->arrayoff + _FPIO_CONST_OFFSET],
+ densize * sizeof (mp_limb_t));
+ }
+ else
+ {
+ cy = mpn_mul (pdest, &__tens[ttab->arrayoff
+ + _FPIO_CONST_OFFSET],
+ ttab->arraysize - _FPIO_CONST_OFFSET,
+ psrc, densize);
+ densize += ttab->arraysize - _FPIO_CONST_OFFSET;
+ if (cy == 0)
+ --densize;
+ (void) SWAP (psrc, pdest);
+ }
+ }
+ expbit <<= 1;
+ ++ttab;
+ }
+ while (neg_exp != 0);
+
+ if (psrc == num)
+ memcpy (den, num, densize * sizeof (mp_limb_t));
+
+ /* Read the fractional digits from the string. */
+ (void) str_to_mpn (startp, dig_no - int_no, num, &numsize, &exponent
+#ifndef USE_WIDE_CHAR
+ , decimal, decimal_len, thousands
+#endif
+ );
+
+ /* We now have to shift both numbers so that the highest bit in the
+ denominator is set. In the same process we copy the numerator to
+ a high place in the array so that the division constructs the wanted
+ digits. This is done by a "quasi fix point" number representation.
+
+ num: ddddddddddd . 0000000000000000000000
+ |--- m ---|
+ den: ddddddddddd n >= m
+ |--- n ---|
+ */
+
+ count_leading_zeros (cnt, den[densize - 1]);
+
+ if (cnt > 0)
+ {
+ /* Don't call `mpn_shift' with a count of zero since the specification
+ does not allow this. */
+ (void) mpn_lshift (den, den, densize, cnt);
+ cy = mpn_lshift (num, num, numsize, cnt);
+ if (cy != 0)
+ num[numsize++] = cy;
+ }
+
+ /* Now we are ready for the division. But it is not necessary to
+ do a full multi-precision division because we only need a small
+ number of bits for the result. So we do not use mpn_divmod
+ here but instead do the division here by hand and stop whenever
+ the needed number of bits is reached. The code itself comes
+ from the GNU MP Library by Torbj\"orn Granlund. */
+
+ exponent = bits;
+
+ switch (densize)
+ {
+ case 1:
+ {
+ mp_limb_t d, n, quot;
+ int used = 0;
+
+ n = num[0];
+ d = den[0];
+ assert (numsize == 1 && n < d);
+
+ do
+ {
+ udiv_qrnnd (quot, n, n, 0, d);
+
+#define got_limb \
+ if (bits == 0) \
+ { \
+ register int cnt; \
+ if (quot == 0) \
+ cnt = BITS_PER_MP_LIMB; \
+ else \
+ count_leading_zeros (cnt, quot); \
+ exponent -= cnt; \
+ if (BITS_PER_MP_LIMB - cnt > MANT_DIG) \
+ { \
+ used = MANT_DIG + cnt; \
+ retval[0] = quot >> (BITS_PER_MP_LIMB - used); \
+ bits = MANT_DIG + 1; \
+ } \
+ else \
+ { \
+ /* Note that we only clear the second element. */ \
+ /* The conditional is determined at compile time. */ \
+ if (RETURN_LIMB_SIZE > 1) \
+ retval[1] = 0; \
+ retval[0] = quot; \
+ bits = -cnt; \
+ } \
+ } \
+ else if (bits + BITS_PER_MP_LIMB <= MANT_DIG) \
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE, BITS_PER_MP_LIMB, \
+ quot); \
+ else \
+ { \
+ used = MANT_DIG - bits; \
+ if (used > 0) \
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, quot); \
+ } \
+ bits += BITS_PER_MP_LIMB
+
+ got_limb;
+ }
+ while (bits <= MANT_DIG);
+
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ more_bits || n != 0);
+ }
+ case 2:
+ {
+ mp_limb_t d0, d1, n0, n1;
+ mp_limb_t quot = 0;
+ int used = 0;
+
+ d0 = den[0];
+ d1 = den[1];
+
+ if (numsize < densize)
+ {
+ if (num[0] >= d1)
+ {
+ /* The numerator of the number occupies fewer bits than
+ the denominator but the one limb is bigger than the
+ high limb of the numerator. */
+ n1 = 0;
+ n0 = num[0];
+ }
+ else
+ {
+ if (bits <= 0)
+ exponent -= BITS_PER_MP_LIMB;
+ else
+ {
+ if (bits + BITS_PER_MP_LIMB <= MANT_DIG)
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
+ BITS_PER_MP_LIMB, 0);
+ else
+ {
+ used = MANT_DIG - bits;
+ if (used > 0)
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
+ }
+ bits += BITS_PER_MP_LIMB;
+ }
+ n1 = num[0];
+ n0 = 0;
+ }
+ }
+ else
+ {
+ n1 = num[1];
+ n0 = num[0];
+ }
+
+ while (bits <= MANT_DIG)
+ {
+ mp_limb_t r;
+
+ if (n1 == d1)
+ {
+ /* QUOT should be either 111..111 or 111..110. We need
+ special treatment of this rare case as normal division
+ would give overflow. */
+ quot = ~(mp_limb_t) 0;
+
+ r = n0 + d1;
+ if (r < d1) /* Carry in the addition? */
+ {
+ add_ssaaaa (n1, n0, r - d0, 0, 0, d0);
+ goto have_quot;
+ }
+ n1 = d0 - (d0 != 0);
+ n0 = -d0;
+ }
+ else
+ {
+ udiv_qrnnd (quot, r, n1, n0, d1);
+ umul_ppmm (n1, n0, d0, quot);
+ }
+
+ q_test:
+ if (n1 > r || (n1 == r && n0 > 0))
+ {
+ /* The estimated QUOT was too large. */
+ --quot;
+
+ sub_ddmmss (n1, n0, n1, n0, 0, d0);
+ r += d1;
+ if (r >= d1) /* If not carry, test QUOT again. */
+ goto q_test;
+ }
+ sub_ddmmss (n1, n0, r, 0, n1, n0);
+
+ have_quot:
+ got_limb;
+ }
+
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ more_bits || n1 != 0 || n0 != 0);
+ }
+ default:
+ {
+ int i;
+ mp_limb_t cy, dX, d1, n0, n1;
+ mp_limb_t quot = 0;
+ int used = 0;
+
+ dX = den[densize - 1];
+ d1 = den[densize - 2];
+
+ /* The division does not work if the upper limb of the two-limb
+ numerator is greater than the denominator. */
+ if (mpn_cmp (num, &den[densize - numsize], numsize) > 0)
+ num[numsize++] = 0;
+
+ if (numsize < densize)
+ {
+ mp_size_t empty = densize - numsize;
+ register int i;
+
+ if (bits <= 0)
+ exponent -= empty * BITS_PER_MP_LIMB;
+ else
+ {
+ if (bits + empty * BITS_PER_MP_LIMB <= MANT_DIG)
+ {
+ /* We make a difference here because the compiler
+ cannot optimize the `else' case that good and
+ this reflects all currently used FLOAT types
+ and GMP implementations. */
+#if RETURN_LIMB_SIZE <= 2
+ assert (empty == 1);
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE,
+ BITS_PER_MP_LIMB, 0);
+#else
+ for (i = RETURN_LIMB_SIZE - 1; i >= empty; --i)
+ retval[i] = retval[i - empty];
+ while (i >= 0)
+ retval[i--] = 0;
+#endif
+ }
+ else
+ {
+ used = MANT_DIG - bits;
+ if (used >= BITS_PER_MP_LIMB)
+ {
+ register int i;
+ (void) mpn_lshift (&retval[used
+ / BITS_PER_MP_LIMB],
+ retval,
+ (RETURN_LIMB_SIZE
+ - used / BITS_PER_MP_LIMB),
+ used % BITS_PER_MP_LIMB);
+ for (i = used / BITS_PER_MP_LIMB - 1; i >= 0; --i)
+ retval[i] = 0;
+ }
+ else if (used > 0)
+ mpn_lshift_1 (retval, RETURN_LIMB_SIZE, used, 0);
+ }
+ bits += empty * BITS_PER_MP_LIMB;
+ }
+ for (i = numsize; i > 0; --i)
+ num[i + empty] = num[i - 1];
+ MPN_ZERO (num, empty + 1);
+ }
+ else
+ {
+ int i;
+ assert (numsize == densize);
+ for (i = numsize; i > 0; --i)
+ num[i] = num[i - 1];
+ }
+
+ den[densize] = 0;
+ n0 = num[densize];
+
+ while (bits <= MANT_DIG)
+ {
+ if (n0 == dX)
+ /* This might over-estimate QUOT, but it's probably not
+ worth the extra code here to find out. */
+ quot = ~(mp_limb_t) 0;
+ else
+ {
+ mp_limb_t r;
+
+ udiv_qrnnd (quot, r, n0, num[densize - 1], dX);
+ umul_ppmm (n1, n0, d1, quot);
+
+ while (n1 > r || (n1 == r && n0 > num[densize - 2]))
+ {
+ --quot;
+ r += dX;
+ if (r < dX) /* I.e. "carry in previous addition?" */
+ break;
+ n1 -= n0 < d1;
+ n0 -= d1;
+ }
+ }
+
+ /* Possible optimization: We already have (q * n0) and (1 * n1)
+ after the calculation of QUOT. Taking advantage of this, we
+ could make this loop make two iterations less. */
+
+ cy = mpn_submul_1 (num, den, densize + 1, quot);
+
+ if (num[densize] != cy)
+ {
+ cy = mpn_add_n (num, num, den, densize);
+ assert (cy != 0);
+ --quot;
+ }
+ n0 = num[densize] = num[densize - 1];
+ for (i = densize - 1; i > 0; --i)
+ num[i] = num[i - 1];
+
+ got_limb;
+ }
+
+ for (i = densize; num[i] == 0 && i >= 0; --i)
+ ;
+ return round_and_return (retval, exponent - 1, negative,
+ quot, BITS_PER_MP_LIMB - 1 - used,
+ more_bits || i >= 0);
+ }
+ }
+ }
+
+ /* NOTREACHED */
+}
--- /dev/null
+/* Copyright (C) 1999, 2004 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, write to the Free
+ Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
+ 02111-1307 USA. */
+
+/* The actual implementation for all floating point sizes is in strtod.c.
+ These macros tell it to produce the `__float128' version, `strtold'. */
+
+#define FLOAT __float128
+#define FLT FLT128
+#ifdef USE_WIDE_CHAR
+# define STRTOF wcstoflt128
+# define __STRTOF __wcstoflt128
+#else
+# define STRTOF strtoflt128
+# define __STRTOF __strtoflt128
+#endif
+#define MPN2FLOAT mpn_construct_float128
+#define FLOAT_HUGE_VAL HUGE_VALQ
+#define SET_MANTISSA(flt, mant) \
+ do { ieee854_float128 u; \
+ u.value = (flt); \
+ u.ieee.mant_high = 0x800000000000ULL; \
+ u.ieee.mant_low = mant; \
+ (flt) = u.value; \
+ } while (0)
+
+static inline __attribute__((__always_inline__))
+__float128 ____strtoflt128_internal (const char *, char **, int);
+
+#include "strtod_l.c"
+
+__float128
+strtoflt128 (const char *nptr, char **endptr)
+{
+ return ____STRTOF_INTERNAL (nptr, endptr, 0);
+}
--- /dev/null
+#include <config.h>
+#include "../printf/gmp-impl.h"
+
+
+/* Definitions according to limb size used. */
+#if BITS_PER_MP_LIMB == 32
+# define MAX_DIG_PER_LIMB 9
+# define MAX_FAC_PER_LIMB 1000000000UL
+#elif BITS_PER_MP_LIMB == 64
+# define MAX_DIG_PER_LIMB 19
+# define MAX_FAC_PER_LIMB 10000000000000000000ULL
+#else
+# error "mp_limb_t size " BITS_PER_MP_LIMB "not accounted for"
+#endif
+
+
+/* Local data structure. */
+const mp_limb_t __quadmath_tens_in_limb[MAX_DIG_PER_LIMB + 1] attribute_hidden
+=
+{ 0, 10, 100,
+ 1000, 10000, 100000L,
+ 1000000L, 10000000L, 100000000L,
+ 1000000000L
+#if BITS_PER_MP_LIMB > 32
+ , 10000000000ULL, 100000000000ULL,
+ 1000000000000ULL, 10000000000000ULL, 100000000000000ULL,
+ 1000000000000000ULL, 10000000000000000ULL, 100000000000000000ULL,
+ 1000000000000000000ULL, 10000000000000000000ULL
+#endif
+#if BITS_PER_MP_LIMB > 64
+ #error "Need to expand tens_in_limb table to" MAX_DIG_PER_LIMB
+#endif
+};
projects / gcc.git / commitdiff