diff options
author | Javier Jardón <jjardon@gnome.org> | 2015-04-27 11:43:02 +0100 |
---|---|---|
committer | Pedro Alvarez <pedro.alvarez@codethink.co.uk> | 2016-05-27 16:24:35 +0100 |
commit | 3c6e2acb7a85497467c87410220ec34e275b44ed (patch) | |
tree | 9071b814720e96d21a3fed8221b2284d8b3debf3 /gmp/mpf | |
parent | 7b48bf2011b4020c4a5a2d5d4149b03983f72cc2 (diff) | |
download | gcc-tarball-3c6e2acb7a85497467c87410220ec34e275b44ed.tar.gz |
Add GMP 6.0.0a tarball
Diffstat (limited to 'gmp/mpf')
70 files changed, 6983 insertions, 0 deletions
diff --git a/gmp/mpf/Makefile.am b/gmp/mpf/Makefile.am new file mode 100644 index 0000000000..30b0ea8f5f --- /dev/null +++ b/gmp/mpf/Makefile.am @@ -0,0 +1,47 @@ +## Process this file with automake to generate Makefile.in + +# Copyright 1996, 1998-2002 Free Software Foundation, Inc. +# +# This file is part of the GNU MP Library. +# +# The GNU MP Library is free software; you can redistribute it and/or modify +# it under the terms of either: +# +# * the GNU Lesser General Public License as published by the Free +# Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# or +# +# * the GNU General Public License as published by the Free Software +# Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# or both in parallel, as here. +# +# The GNU MP 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 General Public License +# for more details. +# +# You should have received copies of the GNU General Public License and the +# GNU Lesser General Public License along with the GNU MP Library. If not, +# see https://www.gnu.org/licenses/. + + +INCLUDES = -D__GMP_WITHIN_GMP -I$(top_srcdir) + +noinst_LTLIBRARIES = libmpf.la +libmpf_la_SOURCES = \ + init.c init2.c inits.c set.c set_ui.c set_si.c set_str.c set_d.c set_z.c \ + set_q.c iset.c iset_ui.c iset_si.c iset_str.c iset_d.c clear.c clears.c \ + get_str.c dump.c size.c eq.c reldiff.c sqrt.c random2.c inp_str.c out_str.c \ + add.c add_ui.c sub.c sub_ui.c ui_sub.c mul.c mul_ui.c div.c div_ui.c \ + cmp.c cmp_d.c cmp_si.c cmp_ui.c mul_2exp.c div_2exp.c abs.c neg.c get_d.c \ + get_d_2exp.c set_dfl_prec.c set_prc.c set_prc_raw.c get_dfl_prec.c get_prc.c \ + ui_div.c sqrt_ui.c \ + pow_ui.c urandomb.c swap.c get_si.c get_ui.c int_p.c \ + ceilfloor.c trunc.c \ + fits_sint.c fits_slong.c fits_sshort.c \ + fits_uint.c fits_ulong.c fits_ushort.c \ + fits_s.h fits_u.h diff --git a/gmp/mpf/Makefile.in b/gmp/mpf/Makefile.in new file mode 100644 index 0000000000..2f5238be58 --- /dev/null +++ b/gmp/mpf/Makefile.in @@ -0,0 +1,575 @@ +# Makefile.in generated by automake 1.11.6 from Makefile.am. +# @configure_input@ + +# Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, +# 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software +# Foundation, Inc. +# This Makefile.in is free software; the Free Software Foundation +# gives unlimited permission to copy and/or distribute it, +# with or without modifications, as long as this notice is preserved. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY, to the extent permitted by law; without +# even the implied warranty of MERCHANTABILITY or FITNESS FOR A +# PARTICULAR PURPOSE. + +@SET_MAKE@ + +# Copyright 1996, 1998-2002 Free Software Foundation, Inc. +# +# This file is part of the GNU MP Library. +# +# The GNU MP Library is free software; you can redistribute it and/or modify +# it under the terms of either: +# +# * the GNU Lesser General Public License as published by the Free +# Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# or +# +# * the GNU General Public License as published by the Free Software +# Foundation; either version 2 of the License, or (at your option) any +# later version. +# +# or both in parallel, as here. +# +# The GNU MP 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 General Public License +# for more details. +# +# You should have received copies of the GNU General Public License and the +# GNU Lesser General Public License along with the GNU MP Library. If not, +# see https://www.gnu.org/licenses/. + +VPATH = @srcdir@ +am__make_dryrun = \ + { \ + am__dry=no; \ + case $$MAKEFLAGS in \ + *\\[\ \ ]*) \ + echo 'am--echo: ; @echo "AM" OK' | $(MAKE) -f - 2>/dev/null \ + | grep '^AM OK$$' >/dev/null || am__dry=yes;; \ + *) \ + for am__flg in $$MAKEFLAGS; do \ + case $$am__flg in \ + *=*|--*) ;; \ + *n*) am__dry=yes; break;; \ + esac; \ + done;; \ + esac; \ + test $$am__dry = yes; \ + } +pkgdatadir = $(datadir)/@PACKAGE@ +pkgincludedir = $(includedir)/@PACKAGE@ +pkglibdir = $(libdir)/@PACKAGE@ +pkglibexecdir = $(libexecdir)/@PACKAGE@ +am__cd = CDPATH="$${ZSH_VERSION+.}$(PATH_SEPARATOR)" && cd +install_sh_DATA = $(install_sh) -c -m 644 +install_sh_PROGRAM = $(install_sh) -c +install_sh_SCRIPT = $(install_sh) -c +INSTALL_HEADER = $(INSTALL_DATA) +transform = $(program_transform_name) +NORMAL_INSTALL = : +PRE_INSTALL = : +POST_INSTALL = : +NORMAL_UNINSTALL = : +PRE_UNINSTALL = : +POST_UNINSTALL = : +build_triplet = @build@ +host_triplet = @host@ +subdir = mpf +DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in +ACLOCAL_M4 = $(top_srcdir)/aclocal.m4 +am__aclocal_m4_deps = $(top_srcdir)/acinclude.m4 \ + $(top_srcdir)/configure.ac +am__configure_deps = $(am__aclocal_m4_deps) $(CONFIGURE_DEPENDENCIES) \ + $(ACLOCAL_M4) +mkinstalldirs = $(install_sh) -d +CONFIG_HEADER = $(top_builddir)/config.h +CONFIG_CLEAN_FILES = +CONFIG_CLEAN_VPATH_FILES = +LTLIBRARIES = $(noinst_LTLIBRARIES) +libmpf_la_LIBADD = +am_libmpf_la_OBJECTS = init.lo init2.lo inits.lo set.lo set_ui.lo \ + set_si.lo set_str.lo set_d.lo set_z.lo set_q.lo iset.lo \ + iset_ui.lo iset_si.lo iset_str.lo iset_d.lo clear.lo clears.lo \ + get_str.lo dump.lo size.lo eq.lo reldiff.lo sqrt.lo random2.lo \ + inp_str.lo out_str.lo add.lo add_ui.lo sub.lo sub_ui.lo \ + ui_sub.lo mul.lo mul_ui.lo div.lo div_ui.lo cmp.lo cmp_d.lo \ + cmp_si.lo cmp_ui.lo mul_2exp.lo div_2exp.lo abs.lo neg.lo \ + get_d.lo get_d_2exp.lo set_dfl_prec.lo set_prc.lo \ + set_prc_raw.lo get_dfl_prec.lo get_prc.lo ui_div.lo sqrt_ui.lo \ + pow_ui.lo urandomb.lo swap.lo get_si.lo get_ui.lo int_p.lo \ + ceilfloor.lo trunc.lo fits_sint.lo fits_slong.lo \ + fits_sshort.lo fits_uint.lo fits_ulong.lo fits_ushort.lo +libmpf_la_OBJECTS = $(am_libmpf_la_OBJECTS) +DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir) +depcomp = +am__depfiles_maybe = +COMPILE = $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) \ + $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +LTCOMPILE = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ + --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) \ + $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) +CCLD = $(CC) +LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) \ + --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) $(AM_LDFLAGS) \ + $(LDFLAGS) -o $@ +SOURCES = $(libmpf_la_SOURCES) +DIST_SOURCES = $(libmpf_la_SOURCES) +am__can_run_installinfo = \ + case $$AM_UPDATE_INFO_DIR in \ + n|no|NO) false;; \ + *) (install-info --version) >/dev/null 2>&1;; \ + esac +ETAGS = etags +CTAGS = ctags +DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST) +ABI = @ABI@ +ACLOCAL = @ACLOCAL@ +AMTAR = @AMTAR@ +AR = @AR@ +AS = @AS@ +ASMFLAGS = @ASMFLAGS@ +AUTOCONF = @AUTOCONF@ +AUTOHEADER = @AUTOHEADER@ +AUTOMAKE = @AUTOMAKE@ +AWK = @AWK@ +CALLING_CONVENTIONS_OBJS = @CALLING_CONVENTIONS_OBJS@ +CC = @CC@ +CCAS = @CCAS@ +CC_FOR_BUILD = @CC_FOR_BUILD@ +CFLAGS = @CFLAGS@ +CPP = @CPP@ +CPPFLAGS = @CPPFLAGS@ +CPP_FOR_BUILD = @CPP_FOR_BUILD@ +CXX = @CXX@ +CXXCPP = @CXXCPP@ +CXXFLAGS = @CXXFLAGS@ +CYGPATH_W = @CYGPATH_W@ +DEFN_LONG_LONG_LIMB = @DEFN_LONG_LONG_LIMB@ +DEFS = @DEFS@ +DLLTOOL = @DLLTOOL@ +DSYMUTIL = @DSYMUTIL@ +DUMPBIN = @DUMPBIN@ +ECHO_C = @ECHO_C@ +ECHO_N = @ECHO_N@ +ECHO_T = @ECHO_T@ +EGREP = @EGREP@ +EXEEXT = @EXEEXT@ +EXEEXT_FOR_BUILD = @EXEEXT_FOR_BUILD@ +FGREP = @FGREP@ +GMP_LDFLAGS = @GMP_LDFLAGS@ +GMP_LIMB_BITS = @GMP_LIMB_BITS@ +GMP_NAIL_BITS = @GMP_NAIL_BITS@ +GREP = @GREP@ +HAVE_CLOCK_01 = @HAVE_CLOCK_01@ +HAVE_CPUTIME_01 = @HAVE_CPUTIME_01@ +HAVE_GETRUSAGE_01 = @HAVE_GETRUSAGE_01@ +HAVE_GETTIMEOFDAY_01 = @HAVE_GETTIMEOFDAY_01@ +HAVE_HOST_CPU_FAMILY_power = @HAVE_HOST_CPU_FAMILY_power@ +HAVE_HOST_CPU_FAMILY_powerpc = @HAVE_HOST_CPU_FAMILY_powerpc@ +HAVE_SIGACTION_01 = @HAVE_SIGACTION_01@ +HAVE_SIGALTSTACK_01 = @HAVE_SIGALTSTACK_01@ +HAVE_SIGSTACK_01 = @HAVE_SIGSTACK_01@ +HAVE_STACK_T_01 = @HAVE_STACK_T_01@ +HAVE_SYS_RESOURCE_H_01 = @HAVE_SYS_RESOURCE_H_01@ +INSTALL = @INSTALL@ +INSTALL_DATA = @INSTALL_DATA@ +INSTALL_PROGRAM = @INSTALL_PROGRAM@ +INSTALL_SCRIPT = @INSTALL_SCRIPT@ +INSTALL_STRIP_PROGRAM = @INSTALL_STRIP_PROGRAM@ +LD = @LD@ +LDFLAGS = @LDFLAGS@ +LEX = @LEX@ +LEXLIB = @LEXLIB@ +LEX_OUTPUT_ROOT = @LEX_OUTPUT_ROOT@ +LIBCURSES = @LIBCURSES@ +LIBGMPXX_LDFLAGS = @LIBGMPXX_LDFLAGS@ +LIBGMP_DLL = @LIBGMP_DLL@ +LIBGMP_LDFLAGS = @LIBGMP_LDFLAGS@ +LIBM = @LIBM@ +LIBM_FOR_BUILD = @LIBM_FOR_BUILD@ +LIBOBJS = @LIBOBJS@ +LIBREADLINE = @LIBREADLINE@ +LIBS = @LIBS@ +LIBTOOL = @LIBTOOL@ +LIPO = @LIPO@ +LN_S = @LN_S@ +LTLIBOBJS = @LTLIBOBJS@ +M4 = @M4@ +MAINT = @MAINT@ +MAKEINFO = @MAKEINFO@ +MANIFEST_TOOL = @MANIFEST_TOOL@ +MKDIR_P = @MKDIR_P@ +NM = @NM@ +NMEDIT = @NMEDIT@ +OBJDUMP = @OBJDUMP@ +OBJEXT = @OBJEXT@ +OTOOL = @OTOOL@ +OTOOL64 = @OTOOL64@ +PACKAGE = @PACKAGE@ +PACKAGE_BUGREPORT = @PACKAGE_BUGREPORT@ +PACKAGE_NAME = @PACKAGE_NAME@ +PACKAGE_STRING = @PACKAGE_STRING@ +PACKAGE_TARNAME = @PACKAGE_TARNAME@ +PACKAGE_URL = @PACKAGE_URL@ +PACKAGE_VERSION = @PACKAGE_VERSION@ +PATH_SEPARATOR = @PATH_SEPARATOR@ +RANLIB = @RANLIB@ +SED = @SED@ +SET_MAKE = @SET_MAKE@ +SHELL = @SHELL@ +SPEED_CYCLECOUNTER_OBJ = @SPEED_CYCLECOUNTER_OBJ@ +STRIP = @STRIP@ +TAL_OBJECT = @TAL_OBJECT@ +TUNE_LIBS = @TUNE_LIBS@ +TUNE_SQR_OBJ = @TUNE_SQR_OBJ@ +U_FOR_BUILD = @U_FOR_BUILD@ +VERSION = @VERSION@ +WITH_READLINE_01 = @WITH_READLINE_01@ +YACC = @YACC@ +YFLAGS = @YFLAGS@ +abs_builddir = @abs_builddir@ +abs_srcdir = @abs_srcdir@ +abs_top_builddir = @abs_top_builddir@ +abs_top_srcdir = @abs_top_srcdir@ +ac_ct_AR = @ac_ct_AR@ +ac_ct_CC = @ac_ct_CC@ +ac_ct_CXX = @ac_ct_CXX@ +ac_ct_DUMPBIN = @ac_ct_DUMPBIN@ +am__leading_dot = @am__leading_dot@ +am__tar = @am__tar@ +am__untar = @am__untar@ +bindir = @bindir@ +build = @build@ +build_alias = @build_alias@ +build_cpu = @build_cpu@ +build_os = @build_os@ +build_vendor = @build_vendor@ +builddir = @builddir@ +datadir = @datadir@ +datarootdir = @datarootdir@ +docdir = @docdir@ +dvidir = @dvidir@ +exec_prefix = @exec_prefix@ +gmp_srclinks = @gmp_srclinks@ +host = @host@ +host_alias = @host_alias@ +host_cpu = @host_cpu@ +host_os = @host_os@ +host_vendor = @host_vendor@ +htmldir = @htmldir@ +includedir = @includedir@ +infodir = @infodir@ +install_sh = @install_sh@ +libdir = @libdir@ +libexecdir = @libexecdir@ +localedir = @localedir@ +localstatedir = @localstatedir@ +mandir = @mandir@ +mkdir_p = @mkdir_p@ +mpn_objects = @mpn_objects@ +mpn_objs_in_libgmp = @mpn_objs_in_libgmp@ +oldincludedir = @oldincludedir@ +pdfdir = @pdfdir@ +prefix = @prefix@ +program_transform_name = @program_transform_name@ +psdir = @psdir@ +sbindir = @sbindir@ +sharedstatedir = @sharedstatedir@ +srcdir = @srcdir@ +sysconfdir = @sysconfdir@ +target_alias = @target_alias@ +top_build_prefix = @top_build_prefix@ +top_builddir = @top_builddir@ +top_srcdir = @top_srcdir@ +INCLUDES = -D__GMP_WITHIN_GMP -I$(top_srcdir) +noinst_LTLIBRARIES = libmpf.la +libmpf_la_SOURCES = \ + init.c init2.c inits.c set.c set_ui.c set_si.c set_str.c set_d.c set_z.c \ + set_q.c iset.c iset_ui.c iset_si.c iset_str.c iset_d.c clear.c clears.c \ + get_str.c dump.c size.c eq.c reldiff.c sqrt.c random2.c inp_str.c out_str.c \ + add.c add_ui.c sub.c sub_ui.c ui_sub.c mul.c mul_ui.c div.c div_ui.c \ + cmp.c cmp_d.c cmp_si.c cmp_ui.c mul_2exp.c div_2exp.c abs.c neg.c get_d.c \ + get_d_2exp.c set_dfl_prec.c set_prc.c set_prc_raw.c get_dfl_prec.c get_prc.c \ + ui_div.c sqrt_ui.c \ + pow_ui.c urandomb.c swap.c get_si.c get_ui.c int_p.c \ + ceilfloor.c trunc.c \ + fits_sint.c fits_slong.c fits_sshort.c \ + fits_uint.c fits_ulong.c fits_ushort.c \ + fits_s.h fits_u.h + +all: all-am + +.SUFFIXES: +.SUFFIXES: .c .lo .o .obj +$(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(am__configure_deps) + @for dep in $?; do \ + case '$(am__configure_deps)' in \ + *$$dep*) \ + ( cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh ) \ + && { if test -f $@; then exit 0; else break; fi; }; \ + exit 1;; \ + esac; \ + done; \ + echo ' cd $(top_srcdir) && $(AUTOMAKE) --gnu --ignore-deps mpf/Makefile'; \ + $(am__cd) $(top_srcdir) && \ + $(AUTOMAKE) --gnu --ignore-deps mpf/Makefile +.PRECIOUS: Makefile +Makefile: $(srcdir)/Makefile.in $(top_builddir)/config.status + @case '$?' in \ + *config.status*) \ + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh;; \ + *) \ + echo ' cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe)'; \ + cd $(top_builddir) && $(SHELL) ./config.status $(subdir)/$@ $(am__depfiles_maybe);; \ + esac; + +$(top_builddir)/config.status: $(top_srcdir)/configure $(CONFIG_STATUS_DEPENDENCIES) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh + +$(top_srcdir)/configure: @MAINTAINER_MODE_TRUE@ $(am__configure_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(ACLOCAL_M4): @MAINTAINER_MODE_TRUE@ $(am__aclocal_m4_deps) + cd $(top_builddir) && $(MAKE) $(AM_MAKEFLAGS) am--refresh +$(am__aclocal_m4_deps): + +clean-noinstLTLIBRARIES: + -test -z "$(noinst_LTLIBRARIES)" || rm -f $(noinst_LTLIBRARIES) + @list='$(noinst_LTLIBRARIES)'; for p in $$list; do \ + dir="`echo $$p | sed -e 's|/[^/]*$$||'`"; \ + test "$$dir" != "$$p" || dir=.; \ + echo "rm -f \"$${dir}/so_locations\""; \ + rm -f "$${dir}/so_locations"; \ + done +libmpf.la: $(libmpf_la_OBJECTS) $(libmpf_la_DEPENDENCIES) $(EXTRA_libmpf_la_DEPENDENCIES) + $(LINK) $(libmpf_la_OBJECTS) $(libmpf_la_LIBADD) $(LIBS) + +mostlyclean-compile: + -rm -f *.$(OBJEXT) + +distclean-compile: + -rm -f *.tab.c + +.c.o: + $(COMPILE) -c $< + +.c.obj: + $(COMPILE) -c `$(CYGPATH_W) '$<'` + +.c.lo: + $(LTCOMPILE) -c -o $@ $< + +mostlyclean-libtool: + -rm -f *.lo + +clean-libtool: + -rm -rf .libs _libs + +ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + mkid -fID $$unique +tags: TAGS + +TAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + set x; \ + here=`pwd`; \ + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + shift; \ + if test -z "$(ETAGS_ARGS)$$*$$unique"; then :; else \ + test -n "$$unique" || unique=$$empty_fix; \ + if test $$# -gt 0; then \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + "$$@" $$unique; \ + else \ + $(ETAGS) $(ETAGSFLAGS) $(AM_ETAGSFLAGS) $(ETAGS_ARGS) \ + $$unique; \ + fi; \ + fi +ctags: CTAGS +CTAGS: $(HEADERS) $(SOURCES) $(TAGS_DEPENDENCIES) \ + $(TAGS_FILES) $(LISP) + list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \ + unique=`for i in $$list; do \ + if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \ + done | \ + $(AWK) '{ files[$$0] = 1; nonempty = 1; } \ + END { if (nonempty) { for (i in files) print i; }; }'`; \ + test -z "$(CTAGS_ARGS)$$unique" \ + || $(CTAGS) $(CTAGSFLAGS) $(AM_CTAGSFLAGS) $(CTAGS_ARGS) \ + $$unique + +GTAGS: + here=`$(am__cd) $(top_builddir) && pwd` \ + && $(am__cd) $(top_srcdir) \ + && gtags -i $(GTAGS_ARGS) "$$here" + +distclean-tags: + -rm -f TAGS ID GTAGS GRTAGS GSYMS GPATH tags + +distdir: $(DISTFILES) + @srcdirstrip=`echo "$(srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + topsrcdirstrip=`echo "$(top_srcdir)" | sed 's/[].[^$$\\*]/\\\\&/g'`; \ + list='$(DISTFILES)'; \ + dist_files=`for file in $$list; do echo $$file; done | \ + sed -e "s|^$$srcdirstrip/||;t" \ + -e "s|^$$topsrcdirstrip/|$(top_builddir)/|;t"`; \ + case $$dist_files in \ + */*) $(MKDIR_P) `echo "$$dist_files" | \ + sed '/\//!d;s|^|$(distdir)/|;s,/[^/]*$$,,' | \ + sort -u` ;; \ + esac; \ + for file in $$dist_files; do \ + if test -f $$file || test -d $$file; then d=.; else d=$(srcdir); fi; \ + if test -d $$d/$$file; then \ + dir=`echo "/$$file" | sed -e 's,/[^/]*$$,,'`; \ + if test -d "$(distdir)/$$file"; then \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + if test -d $(srcdir)/$$file && test $$d != $(srcdir); then \ + cp -fpR $(srcdir)/$$file "$(distdir)$$dir" || exit 1; \ + find "$(distdir)/$$file" -type d ! -perm -700 -exec chmod u+rwx {} \;; \ + fi; \ + cp -fpR $$d/$$file "$(distdir)$$dir" || exit 1; \ + else \ + test -f "$(distdir)/$$file" \ + || cp -p $$d/$$file "$(distdir)/$$file" \ + || exit 1; \ + fi; \ + done +check-am: all-am +check: check-am +all-am: Makefile $(LTLIBRARIES) +installdirs: +install: install-am +install-exec: install-exec-am +install-data: install-data-am +uninstall: uninstall-am + +install-am: all-am + @$(MAKE) $(AM_MAKEFLAGS) install-exec-am install-data-am + +installcheck: installcheck-am +install-strip: + if test -z '$(STRIP)'; then \ + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + install; \ + else \ + $(MAKE) $(AM_MAKEFLAGS) INSTALL_PROGRAM="$(INSTALL_STRIP_PROGRAM)" \ + install_sh_PROGRAM="$(INSTALL_STRIP_PROGRAM)" INSTALL_STRIP_FLAG=-s \ + "INSTALL_PROGRAM_ENV=STRIPPROG='$(STRIP)'" install; \ + fi +mostlyclean-generic: + +clean-generic: + +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) + +maintainer-clean-generic: + @echo "This command is intended for maintainers to use" + @echo "it deletes files that may require special tools to rebuild." +clean: clean-am + +clean-am: clean-generic clean-libtool clean-noinstLTLIBRARIES \ + mostlyclean-am + +distclean: distclean-am + -rm -f Makefile +distclean-am: clean-am distclean-compile distclean-generic \ + distclean-tags + +dvi: dvi-am + +dvi-am: + +html: html-am + +html-am: + +info: info-am + +info-am: + +install-data-am: + +install-dvi: install-dvi-am + +install-dvi-am: + +install-exec-am: + +install-html: install-html-am + +install-html-am: + +install-info: install-info-am + +install-info-am: + +install-man: + +install-pdf: install-pdf-am + +install-pdf-am: + +install-ps: install-ps-am + +install-ps-am: + +installcheck-am: + +maintainer-clean: maintainer-clean-am + -rm -f Makefile +maintainer-clean-am: distclean-am maintainer-clean-generic + +mostlyclean: mostlyclean-am + +mostlyclean-am: mostlyclean-compile mostlyclean-generic \ + mostlyclean-libtool + +pdf: pdf-am + +pdf-am: + +ps: ps-am + +ps-am: + +uninstall-am: + +.MAKE: install-am install-strip + +.PHONY: CTAGS GTAGS all all-am check check-am clean clean-generic \ + clean-libtool clean-noinstLTLIBRARIES ctags distclean \ + distclean-compile distclean-generic distclean-libtool \ + distclean-tags distdir dvi dvi-am html html-am info info-am \ + install install-am install-data install-data-am install-dvi \ + install-dvi-am install-exec install-exec-am install-html \ + install-html-am install-info install-info-am install-man \ + install-pdf install-pdf-am install-ps install-ps-am \ + install-strip installcheck installcheck-am installdirs \ + maintainer-clean maintainer-clean-generic mostlyclean \ + mostlyclean-compile mostlyclean-generic mostlyclean-libtool \ + pdf pdf-am ps ps-am tags uninstall uninstall-am + + +# Tell versions [3.59,3.63) of GNU make to not export all variables. +# Otherwise a system limit (for SysV at least) may be exceeded. +.NOEXPORT: diff --git a/gmp/mpf/abs.c b/gmp/mpf/abs.c new file mode 100644 index 0000000000..a2bde2a4f8 --- /dev/null +++ b/gmp/mpf/abs.c @@ -0,0 +1,59 @@ +/* mpf_abs -- Compute the absolute value of a float. + +Copyright 1993-1995, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_abs (mpf_ptr r, mpf_srcptr u) +{ + mp_size_t size; + + size = ABS (u->_mp_size); + if (r != u) + { + mp_size_t prec; + mp_ptr rp, up; + + prec = r->_mp_prec + 1; /* lie not to lose precision in assignment */ + rp = r->_mp_d; + up = u->_mp_d; + + if (size > prec) + { + up += size - prec; + size = prec; + } + + MPN_COPY (rp, up, size); + r->_mp_exp = u->_mp_exp; + } + r->_mp_size = size; +} diff --git a/gmp/mpf/add.c b/gmp/mpf/add.c new file mode 100644 index 0000000000..d2a5c097c5 --- /dev/null +++ b/gmp/mpf/add.c @@ -0,0 +1,184 @@ +/* mpf_add -- Add two floats. + +Copyright 1993, 1994, 1996, 2000, 2001, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_add (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +{ + mp_srcptr up, vp; + mp_ptr rp, tp; + mp_size_t usize, vsize, rsize; + mp_size_t prec; + mp_exp_t uexp; + mp_size_t ediff; + mp_limb_t cy; + int negate; + TMP_DECL; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (usize == 0) + { + set_r_v_maybe: + if (r != v) + mpf_set (r, v); + return; + } + if (vsize == 0) + { + v = u; + goto set_r_v_maybe; + } + + /* If signs of U and V are different, perform subtraction. */ + if ((usize ^ vsize) < 0) + { + __mpf_struct v_negated; + v_negated._mp_size = -vsize; + v_negated._mp_exp = v->_mp_exp; + v_negated._mp_d = v->_mp_d; + mpf_sub (r, u, &v_negated); + return; + } + + TMP_MARK; + + /* Signs are now known to be the same. */ + negate = usize < 0; + + /* Make U be the operand with the largest exponent. */ + if (u->_mp_exp < v->_mp_exp) + { + mpf_srcptr t; + t = u; u = v; v = t; + usize = u->_mp_size; + vsize = v->_mp_size; + } + + usize = ABS (usize); + vsize = ABS (vsize); + up = u->_mp_d; + vp = v->_mp_d; + rp = r->_mp_d; + prec = r->_mp_prec; + uexp = u->_mp_exp; + ediff = u->_mp_exp - v->_mp_exp; + + /* If U extends beyond PREC, ignore the part that does. */ + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + + /* If V extends beyond PREC, ignore the part that does. + Note that this may make vsize negative. */ + if (vsize + ediff > prec) + { + vp += vsize + ediff - prec; + vsize = prec - ediff; + } + +#if 0 + /* Locate the least significant non-zero limb in (the needed parts + of) U and V, to simplify the code below. */ + while (up[0] == 0) + up++, usize--; + while (vp[0] == 0) + vp++, vsize--; +#endif + + /* Allocate temp space for the result. Allocate + just vsize + ediff later??? */ + tp = TMP_ALLOC_LIMBS (prec); + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (rp != up) + MPN_COPY_INCR (rp, up, usize); + rsize = usize; + } + else + { + /* uuuu | uuuu | uuuu | uuuu | uuuu */ + /* vvvvvvv | vv | vvvvv | v | vv */ + + if (usize > ediff) + { + /* U and V partially overlaps. */ + if (vsize + ediff <= usize) + { + /* uuuu */ + /* v */ + mp_size_t size; + size = usize - ediff - vsize; + MPN_COPY (tp, up, size); + cy = mpn_add (tp + size, up + size, usize - size, vp, vsize); + rsize = usize; + } + else + { + /* uuuu */ + /* vvvvv */ + mp_size_t size; + size = vsize + ediff - usize; + MPN_COPY (tp, vp, size); + cy = mpn_add (tp + size, up, usize, vp + size, usize - ediff); + rsize = vsize + ediff; + } + } + else + { + /* uuuu */ + /* vv */ + mp_size_t size; + size = vsize + ediff - usize; + MPN_COPY (tp, vp, vsize); + MPN_ZERO (tp + vsize, ediff - usize); + MPN_COPY (tp + size, up, usize); + cy = 0; + rsize = size + usize; + } + + MPN_COPY (rp, tp, rsize); + rp[rsize] = cy; + rsize += cy; + uexp += cy; + } + + r->_mp_size = negate ? -rsize : rsize; + r->_mp_exp = uexp; + TMP_FREE; +} diff --git a/gmp/mpf/add_ui.c b/gmp/mpf/add_ui.c new file mode 100644 index 0000000000..b1e57d04c1 --- /dev/null +++ b/gmp/mpf/add_ui.c @@ -0,0 +1,153 @@ +/* mpf_add_ui -- Add a float and an unsigned integer. + +Copyright 1993, 1994, 1996, 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_add_ui (mpf_ptr sum, mpf_srcptr u, unsigned long int v) +{ + mp_srcptr up = u->_mp_d; + mp_ptr sump = sum->_mp_d; + mp_size_t usize, sumsize; + mp_size_t prec = sum->_mp_prec; + mp_exp_t uexp = u->_mp_exp; + + usize = u->_mp_size; + if (usize <= 0) + { + if (usize == 0) + { + mpf_set_ui (sum, v); + return; + } + else + { + __mpf_struct u_negated; + u_negated._mp_size = -usize; + u_negated._mp_exp = u->_mp_exp; + u_negated._mp_d = u->_mp_d; + mpf_sub_ui (sum, &u_negated, v); + sum->_mp_size = -(sum->_mp_size); + return; + } + } + + if (v == 0) + { + sum_is_u: + if (u != sum) + { + sumsize = MIN (usize, prec + 1); + MPN_COPY (sum->_mp_d, up + usize - sumsize, sumsize); + sum->_mp_size = sumsize; + sum->_mp_exp = u->_mp_exp; + } + return; + } + + if (uexp > 0) + { + /* U >= 1. */ + if (uexp > prec) + { + /* U >> V, V is not part of final result. */ + goto sum_is_u; + } + else + { + /* U's "limb point" is somewhere between the first limb + and the PREC:th limb. + Both U and V are part of the final result. */ + if (uexp > usize) + { + /* uuuuuu0000. */ + /* + v. */ + /* We begin with moving U to the top of SUM, to handle + samevar(U,SUM). */ + MPN_COPY_DECR (sump + uexp - usize, up, usize); + sump[0] = v; + MPN_ZERO (sump + 1, uexp - usize - 1); +#if 0 /* What is this??? */ + if (sum == u) + MPN_COPY (sum->_mp_d, sump, uexp); +#endif + sum->_mp_size = uexp; + sum->_mp_exp = uexp; + } + else + { + /* uuuuuu.uuuu */ + /* + v. */ + mp_limb_t cy_limb; + if (usize > prec) + { + /* Ignore excess limbs in U. */ + up += usize - prec; + usize -= usize - prec; /* Eq. usize = prec */ + } + if (sump != up) + MPN_COPY_INCR (sump, up, usize - uexp); + cy_limb = mpn_add_1 (sump + usize - uexp, up + usize - uexp, + uexp, (mp_limb_t) v); + sump[usize] = cy_limb; + sum->_mp_size = usize + cy_limb; + sum->_mp_exp = uexp + cy_limb; + } + } + } + else + { + /* U < 1, so V > U for sure. */ + /* v. */ + /* .0000uuuu */ + if ((-uexp) >= prec) + { + sump[0] = v; + sum->_mp_size = 1; + sum->_mp_exp = 1; + } + else + { + if (usize + (-uexp) + 1 > prec) + { + /* Ignore excess limbs in U. */ + up += usize + (-uexp) + 1 - prec; + usize -= usize + (-uexp) + 1 - prec; + } + if (sump != up) + MPN_COPY_INCR (sump, up, usize); + MPN_ZERO (sump + usize, -uexp); + sump[usize + (-uexp)] = v; + sum->_mp_size = usize + (-uexp) + 1; + sum->_mp_exp = 1; + } + } +} diff --git a/gmp/mpf/ceilfloor.c b/gmp/mpf/ceilfloor.c new file mode 100644 index 0000000000..302e2b8ae5 --- /dev/null +++ b/gmp/mpf/ceilfloor.c @@ -0,0 +1,126 @@ +/* mpf_ceil, mpf_floor -- round an mpf to an integer. + +Copyright 2001, 2004, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* dir==1 for ceil, dir==-1 for floor + + Notice the use of prec+1 ensures mpf_ceil and mpf_floor are equivalent to + mpf_set if u is already an integer. */ + +static void __gmpf_ceil_or_floor (REGPARM_2_1 (mpf_ptr, mpf_srcptr, int)) REGPARM_ATTR (1); +#define mpf_ceil_or_floor(r,u,dir) __gmpf_ceil_or_floor (REGPARM_2_1 (r, u, dir)) + +REGPARM_ATTR (1) static void +mpf_ceil_or_floor (mpf_ptr r, mpf_srcptr u, int dir) +{ + mp_ptr rp, up, p; + mp_size_t size, asize, prec; + mp_exp_t exp; + + size = SIZ(u); + if (size == 0) + { + zero: + SIZ(r) = 0; + EXP(r) = 0; + return; + } + + rp = PTR(r); + exp = EXP(u); + if (exp <= 0) + { + /* u is only a fraction */ + if ((size ^ dir) < 0) + goto zero; + rp[0] = 1; + EXP(r) = 1; + SIZ(r) = dir; + return; + } + EXP(r) = exp; + + up = PTR(u); + asize = ABS (size); + up += asize; + + /* skip fraction part of u */ + asize = MIN (asize, exp); + + /* don't lose precision in the copy */ + prec = PREC (r) + 1; + + /* skip excess over target precision */ + asize = MIN (asize, prec); + + up -= asize; + + if ((size ^ dir) >= 0) + { + /* rounding direction matches sign, must increment if ignored part is + non-zero */ + for (p = PTR(u); p != up; p++) + { + if (*p != 0) + { + if (mpn_add_1 (rp, up, asize, CNST_LIMB(1))) + { + /* was all 0xFF..FFs, which have become zeros, giving just + a carry */ + rp[0] = 1; + asize = 1; + EXP(r)++; + } + SIZ(r) = (size >= 0 ? asize : -asize); + return; + } + } + } + + SIZ(r) = (size >= 0 ? asize : -asize); + if (rp != up) + MPN_COPY_INCR (rp, up, asize); +} + + +void +mpf_ceil (mpf_ptr r, mpf_srcptr u) +{ + mpf_ceil_or_floor (r, u, 1); +} + +void +mpf_floor (mpf_ptr r, mpf_srcptr u) +{ + mpf_ceil_or_floor (r, u, -1); +} diff --git a/gmp/mpf/clear.c b/gmp/mpf/clear.c new file mode 100644 index 0000000000..2df0de579b --- /dev/null +++ b/gmp/mpf/clear.c @@ -0,0 +1,39 @@ +/* mpf_clear -- de-allocate the space occupied by the dynamic digit space of + an integer. + +Copyright 1993-1995, 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_clear (mpf_ptr m) +{ + (*__gmp_free_func) (m->_mp_d, (size_t) (m->_mp_prec + 1) * GMP_LIMB_BYTES); +} diff --git a/gmp/mpf/clears.c b/gmp/mpf/clears.c new file mode 100644 index 0000000000..addbe8faf3 --- /dev/null +++ b/gmp/mpf/clears.c @@ -0,0 +1,49 @@ +/* mpf_clears() -- Clear multiple mpf_t variables. + +Copyright 2009, 2014 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdarg.h> +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_clears (mpf_ptr x, ...) +{ + va_list ap; + + va_start (ap, x); + + while (x != NULL) + { + (*__gmp_free_func) (x->_mp_d, (size_t) (x->_mp_prec + 1) * GMP_LIMB_BYTES); + x = va_arg (ap, mpf_ptr); + } + va_end (ap); +} diff --git a/gmp/mpf/cmp.c b/gmp/mpf/cmp.c new file mode 100644 index 0000000000..ab22c3f89c --- /dev/null +++ b/gmp/mpf/cmp.c @@ -0,0 +1,117 @@ +/* mpf_cmp -- Compare two floats. + +Copyright 1993, 1994, 1996, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +mpf_cmp (mpf_srcptr u, mpf_srcptr v) __GMP_NOTHROW +{ + mp_srcptr up, vp; + mp_size_t usize, vsize; + mp_exp_t uexp, vexp; + int cmp; + int usign; + + uexp = u->_mp_exp; + vexp = v->_mp_exp; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* 1. Are the signs different? */ + if ((usize ^ vsize) >= 0) + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + /* vsize >= 0 */ + return -(vsize != 0); + if (vsize == 0) + /* usize >= 0 */ + return usize != 0; + /* Fall out. */ + } + else + { + /* Either U or V is negative, but not both. */ + return usize >= 0 ? 1 : -1; + } + + /* U and V have the same sign and are both non-zero. */ + + usign = usize >= 0 ? 1 : -1; + + /* 2. Are the exponents different? */ + if (uexp > vexp) + return usign; + if (uexp < vexp) + return -usign; + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + +#define STRICT_MPF_NORMALIZATION 0 +#if ! STRICT_MPF_NORMALIZATION + /* Ignore zeroes at the low end of U and V. */ + while (up[0] == 0) + { + up++; + usize--; + } + while (vp[0] == 0) + { + vp++; + vsize--; + } +#endif + + if (usize > vsize) + { + cmp = mpn_cmp (up + usize - vsize, vp, vsize); + if (cmp == 0) + return usign; + } + else if (vsize > usize) + { + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp == 0) + return -usign; + } + else + { + cmp = mpn_cmp (up, vp, usize); + if (cmp == 0) + return 0; + } + return cmp > 0 ? usign : -usign; +} diff --git a/gmp/mpf/cmp_d.c b/gmp/mpf/cmp_d.c new file mode 100644 index 0000000000..52893a781e --- /dev/null +++ b/gmp/mpf/cmp_d.c @@ -0,0 +1,60 @@ +/* mpf_cmp_d -- compare mpf and double. + +Copyright 2001, 2003 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "config.h" + +#if HAVE_FLOAT_H +#include <float.h> /* for DBL_MAX */ +#endif + +#include "gmp.h" +#include "gmp-impl.h" + +int +mpf_cmp_d (mpf_srcptr f, double d) +{ + mp_limb_t darray[LIMBS_PER_DOUBLE]; + mpf_t df; + + /* d=NaN has no sensible return value, so raise an exception. + d=Inf or -Inf is always bigger than z. */ + DOUBLE_NAN_INF_ACTION (d, + __gmp_invalid_operation (), + return (d < 0.0 ? 1 : -1)); + + if (d == 0.0) + return SIZ(f); + + PTR(df) = darray; + SIZ(df) = (d >= 0.0 ? LIMBS_PER_DOUBLE : -LIMBS_PER_DOUBLE); + EXP(df) = __gmp_extract_double (darray, ABS(d)); + + return mpf_cmp (f, df); +} diff --git a/gmp/mpf/cmp_si.c b/gmp/mpf/cmp_si.c new file mode 100644 index 0000000000..eaa8b87da9 --- /dev/null +++ b/gmp/mpf/cmp_si.c @@ -0,0 +1,118 @@ +/* mpf_cmp_si -- Compare a float with a signed integer. + +Copyright 1993-1995, 1999-2002, 2004, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +mpf_cmp_si (mpf_srcptr u, long int vval) __GMP_NOTHROW +{ + mp_srcptr up; + mp_size_t usize; + mp_exp_t uexp; + mp_limb_t ulimb; + int usign; + unsigned long abs_vval; + + uexp = u->_mp_exp; + usize = u->_mp_size; + + /* 1. Are the signs different? */ + if ((usize < 0) == (vval < 0)) /* don't use xor, type size may differ */ + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + /* vval >= 0 */ + return -(vval != 0); + if (vval == 0) + /* usize >= 0 */ + return usize != 0; + /* Fall out. */ + } + else + { + /* Either U or V is negative, but not both. */ + return usize >= 0 ? 1 : -1; + } + + /* U and V have the same sign and are both non-zero. */ + + usign = usize >= 0 ? 1 : -1; + usize = ABS (usize); + abs_vval = ABS_CAST (unsigned long, vval); + + /* 2. Are the exponents different (V's exponent == 1)? */ +#if GMP_NAIL_BITS != 0 + if (uexp > 1 + (abs_vval > GMP_NUMB_MAX)) + return usign; + if (uexp < 1 + (abs_vval > GMP_NUMB_MAX)) + return -usign; +#else + if (uexp > 1) + return usign; + if (uexp < 1) + return -usign; +#endif + + up = u->_mp_d; + + ulimb = up[usize - 1]; +#if GMP_NAIL_BITS != 0 + if (usize >= 2 && uexp == 2) + { + if ((ulimb >> GMP_NAIL_BITS) != 0) + return usign; + ulimb = (ulimb << GMP_NUMB_BITS) | up[usize - 2]; + usize--; + } +#endif + usize--; + + /* 3. Compare the most significant mantissa limb with V. */ + if (ulimb > abs_vval) + return usign; + else if (ulimb < abs_vval) + return -usign; + + /* Ignore zeroes at the low end of U. */ + while (*up == 0) + { + up++; + usize--; + } + + /* 4. Now, if the number of limbs are different, we have a difference + since we have made sure the trailing limbs are not zero. */ + if (usize > 0) + return usign; + + /* Wow, we got zero even if we tried hard to avoid it. */ + return 0; +} diff --git a/gmp/mpf/cmp_ui.c b/gmp/mpf/cmp_ui.c new file mode 100644 index 0000000000..ccb76c6ce0 --- /dev/null +++ b/gmp/mpf/cmp_ui.c @@ -0,0 +1,100 @@ +/* mpf_cmp_ui -- Compare a float with an unsigned integer. + +Copyright 1993-1995, 1999, 2001, 2002 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +mpf_cmp_ui (mpf_srcptr u, unsigned long int vval) __GMP_NOTHROW +{ + mp_srcptr up; + mp_size_t usize; + mp_exp_t uexp; + mp_limb_t ulimb; + + uexp = u->_mp_exp; + usize = u->_mp_size; + + /* 1. Is U negative? */ + if (usize < 0) + return -1; + /* We rely on usize being non-negative in the code that follows. */ + + if (vval == 0) + return usize != 0; + + /* 2. Are the exponents different (V's exponent == 1)? */ +#if GMP_NAIL_BITS != 0 + if (uexp > 1 + (vval > GMP_NUMB_MAX)) + return 1; + if (uexp < 1 + (vval > GMP_NUMB_MAX)) + return -1; +#else + if (uexp > 1) + return 1; + if (uexp < 1) + return -1; +#endif + + up = u->_mp_d; + + ulimb = up[usize - 1]; +#if GMP_NAIL_BITS != 0 + if (usize >= 2 && uexp == 2) + { + if ((ulimb >> GMP_NAIL_BITS) != 0) + return 1; + ulimb = (ulimb << GMP_NUMB_BITS) | up[usize - 2]; + usize--; + } +#endif + usize--; + + /* 3. Compare the most significant mantissa limb with V. */ + if (ulimb > vval) + return 1; + else if (ulimb < vval) + return -1; + + /* Ignore zeroes at the low end of U. */ + while (*up == 0) + { + up++; + usize--; + } + + /* 4. Now, if the number of limbs are different, we have a difference + since we have made sure the trailing limbs are not zero. */ + if (usize > 0) + return 1; + + /* Wow, we got zero even if we tried hard to avoid it. */ + return 0; +} diff --git a/gmp/mpf/div.c b/gmp/mpf/div.c new file mode 100644 index 0000000000..af38cb8698 --- /dev/null +++ b/gmp/mpf/div.c @@ -0,0 +1,138 @@ +/* mpf_div -- Divide two floats. + +Copyright 1993, 1994, 1996, 2000-2002, 2004, 2005, 2010, 2012 Free Software +Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Not done: + + No attempt is made to identify an overlap u==v. The result will be + correct (1.0), but a full actual division is done whereas of course + x/x==1 needs no work. Such a call is not a sensible thing to make, and + it's left to an application to notice and optimize if it might arise + somehow through pointer aliasing or whatever. + + Enhancements: + + The high quotient limb is non-zero when high{up,vsize} >= {vp,vsize}. We + could make that comparison and use qsize==prec instead of qsize==prec+1, + to save one limb in the division. + + If r==u but the size is enough bigger than prec that there won't be an + overlap between quotient and dividend in mpn_div_q, then we can avoid + copying up,usize. This would only arise from a prec reduced with + mpf_set_prec_raw and will be pretty unusual, but might be worthwhile if + it could be worked into the copy_u decision cleanly. */ + +void +mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +{ + mp_srcptr up, vp; + mp_ptr rp, tp, new_vp; + mp_size_t usize, vsize, rsize, prospective_rsize, tsize, zeros; + mp_size_t sign_quotient, prec, high_zero, chop; + mp_exp_t rexp; + int copy_u; + TMP_DECL; + + usize = SIZ(u); + vsize = SIZ(v); + + if (UNLIKELY (vsize == 0)) + DIVIDE_BY_ZERO; + + if (usize == 0) + { + SIZ(r) = 0; + EXP(r) = 0; + return; + } + + sign_quotient = usize ^ vsize; + usize = ABS (usize); + vsize = ABS (vsize); + prec = PREC(r); + + TMP_MARK; + rexp = EXP(u) - EXP(v) + 1; + + rp = PTR(r); + up = PTR(u); + vp = PTR(v); + + prospective_rsize = usize - vsize + 1; /* quot from using given u,v sizes */ + rsize = prec + 1; /* desired quot */ + + zeros = rsize - prospective_rsize; /* padding u to give rsize */ + copy_u = (zeros > 0 || rp == up); /* copy u if overlap or padding */ + + chop = MAX (-zeros, 0); /* negative zeros means shorten u */ + up += chop; + usize -= chop; + zeros += chop; /* now zeros >= 0 */ + + tsize = usize + zeros; /* size for possible copy of u */ + + /* copy and possibly extend u if necessary */ + if (copy_u) + { + tp = TMP_ALLOC_LIMBS (tsize + 1); /* +1 for mpn_div_q's scratch needs */ + MPN_ZERO (tp, zeros); + MPN_COPY (tp+zeros, up, usize); + up = tp; + usize = tsize; + } + else + { + tp = TMP_ALLOC_LIMBS (usize + 1); + } + + /* ensure divisor doesn't overlap quotient */ + if (rp == vp) + { + new_vp = TMP_ALLOC_LIMBS (vsize); + MPN_COPY (new_vp, vp, vsize); + vp = new_vp; + } + + ASSERT (usize-vsize+1 == rsize); + mpn_div_q (rp, up, usize, vp, vsize, tp); + + /* strip possible zero high limb */ + high_zero = (rp[rsize-1] == 0); + rsize -= high_zero; + rexp -= high_zero; + + SIZ(r) = sign_quotient >= 0 ? rsize : -rsize; + EXP(r) = rexp; + TMP_FREE; +} diff --git a/gmp/mpf/div_2exp.c b/gmp/mpf/div_2exp.c new file mode 100644 index 0000000000..fef8152050 --- /dev/null +++ b/gmp/mpf/div_2exp.c @@ -0,0 +1,139 @@ +/* mpf_div_2exp -- Divide a float by 2^n. + +Copyright 1993, 1994, 1996, 2000-2002, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Multiples of GMP_NUMB_BITS in exp simply mean an amount subtracted from + EXP(u) to set EXP(r). The remainder exp%GMP_NUMB_BITS is then a right + shift for the limb data. + + If exp%GMP_NUMB_BITS == 0 then there's no shifting, we effectively just + do an mpz_set with changed EXP(r). Like mpz_set we take prec+1 limbs in + this case. Although just prec would suffice, it's nice to have + mpf_div_2exp with exp==0 come out the same as mpz_set. + + When shifting we take up to prec many limbs from the input. Our shift is + cy = mpn_rshift (PTR(r)+1, PTR(u)+k, ...), where k is the number of low + limbs dropped from u, and the carry out is stored to PTR(r)[0]. We don't + try to work extra bits from PTR(u)[k-1] (when k>=1 makes it available) + into that low carry limb. Just prec limbs (with the high non-zero) from + the input is enough bits for the application requested precision, no need + to do extra work. + + If r==u the shift will have overlapping operands. When k>=1 (ie. when + usize > prec), the overlap is in the style supported by rshift (ie. dst + <= src). + + But when r==u and k==0 (ie. usize <= prec), we would have an invalid + overlap (mpn_rshift (rp+1, rp, ...)). In this case we must instead use + mpn_lshift (PTR(r), PTR(u), size, NUMB-shift). An lshift by NUMB-shift + bits gives identical data of course, it's just its overlap restrictions + which differ. + + In both shift cases, the resulting data is abs_usize+1 limbs. "adj" is + used to add +1 to that size if the high is non-zero (it may of course + have become zero by the shifting). EXP(u) is the exponent just above + those abs_usize+1 limbs, so it gets -1+adj, which means -1 if the high is + zero, or no change if the high is non-zero. + + Enhancements: + + The way mpn_lshift is used means successive mpf_div_2exp calls on the + same operand will accumulate low zero limbs, until prec+1 limbs is + reached. This is wasteful for subsequent operations. When abs_usize <= + prec, we should test the low exp%GMP_NUMB_BITS many bits of PTR(u)[0], + ie. those which would be shifted out by an mpn_rshift. If they're zero + then use that mpn_rshift. */ + +void +mpf_div_2exp (mpf_ptr r, mpf_srcptr u, mp_bitcnt_t exp) +{ + mp_srcptr up; + mp_ptr rp = r->_mp_d; + mp_size_t usize; + mp_size_t abs_usize; + mp_size_t prec = r->_mp_prec; + mp_exp_t uexp = u->_mp_exp; + + usize = u->_mp_size; + + if (UNLIKELY (usize == 0)) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + abs_usize = ABS (usize); + up = u->_mp_d; + + if (exp % GMP_NUMB_BITS == 0) + { + prec++; /* retain more precision here as we don't need + to account for carry-out here */ + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + } + if (rp != up) + MPN_COPY_INCR (rp, up, abs_usize); + r->_mp_exp = uexp - exp / GMP_NUMB_BITS; + } + else + { + mp_limb_t cy_limb; + mp_size_t adj; + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + /* Use mpn_rshift since mpn_lshift operates downwards, and we + therefore would clobber part of U before using that part, in case + R is the same variable as U. */ + cy_limb = mpn_rshift (rp + 1, up, abs_usize, exp % GMP_NUMB_BITS); + rp[0] = cy_limb; + adj = rp[abs_usize] != 0; + } + else + { + cy_limb = mpn_lshift (rp, up, abs_usize, + GMP_NUMB_BITS - exp % GMP_NUMB_BITS); + rp[abs_usize] = cy_limb; + adj = cy_limb != 0; + } + + abs_usize += adj; + r->_mp_exp = uexp - exp / GMP_NUMB_BITS - 1 + adj; + } + r->_mp_size = usize >= 0 ? abs_usize : -abs_usize; +} diff --git a/gmp/mpf/div_ui.c b/gmp/mpf/div_ui.c new file mode 100644 index 0000000000..9be7e680be --- /dev/null +++ b/gmp/mpf/div_ui.c @@ -0,0 +1,111 @@ +/* mpf_div_ui -- Divide a float with an unsigned integer. + +Copyright 1993, 1994, 1996, 2000-2002, 2004, 2005, 2012 Free Software +Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +mpf_div_ui (mpf_ptr r, mpf_srcptr u, unsigned long int v) +{ + mp_srcptr up; + mp_ptr rp, tp, rtp; + mp_size_t usize; + mp_size_t rsize, tsize; + mp_size_t sign_quotient; + mp_size_t prec; + mp_limb_t q_limb; + mp_exp_t rexp; + TMP_DECL; + +#if BITS_PER_ULONG > GMP_NUMB_BITS /* avoid warnings about shift amount */ + if (v > GMP_NUMB_MAX) + { + mpf_t vf; + mp_limb_t vl[2]; + SIZ(vf) = 2; + EXP(vf) = 2; + PTR(vf) = vl; + vl[0] = v & GMP_NUMB_MASK; + vl[1] = v >> GMP_NUMB_BITS; + mpf_div (r, u, vf); + return; + } +#endif + + if (UNLIKELY (v == 0)) + DIVIDE_BY_ZERO; + + usize = u->_mp_size; + + if (usize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + sign_quotient = usize; + usize = ABS (usize); + prec = r->_mp_prec; + + TMP_MARK; + + rp = r->_mp_d; + up = u->_mp_d; + + tsize = 1 + prec; + tp = TMP_ALLOC_LIMBS (tsize + 1); + + if (usize > tsize) + { + up += usize - tsize; + usize = tsize; + rtp = tp; + } + else + { + MPN_ZERO (tp, tsize - usize); + rtp = tp + (tsize - usize); + } + + /* Move the dividend to the remainder. */ + MPN_COPY (rtp, up, usize); + + mpn_divmod_1 (rp, tp, tsize, (mp_limb_t) v); + q_limb = rp[tsize - 1]; + + rsize = tsize - (q_limb == 0); + rexp = u->_mp_exp - (q_limb == 0); + r->_mp_size = sign_quotient >= 0 ? rsize : -rsize; + r->_mp_exp = rexp; + TMP_FREE; +} diff --git a/gmp/mpf/dump.c b/gmp/mpf/dump.c new file mode 100644 index 0000000000..af67105923 --- /dev/null +++ b/gmp/mpf/dump.c @@ -0,0 +1,53 @@ +/* mpf_dump -- Dump a float to stdout. + + THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE. IT IS NOT SAFE TO + CALL THIS FUNCTION DIRECTLY. IN FACT, IT IS ALMOST GUARANTEED THAT THIS + FUNCTION WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE. + + +Copyright 1993-1995, 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> +#include <string.h> /* for strlen */ +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_dump (mpf_srcptr u) +{ + mp_exp_t exp; + char *str; + + str = mpf_get_str (0, &exp, 10, 0, u); + if (str[0] == '-') + printf ("-0.%se%ld\n", str + 1, exp); + else + printf ("0.%se%ld\n", str, exp); + (*__gmp_free_func) (str, strlen (str) + 1); +} diff --git a/gmp/mpf/eq.c b/gmp/mpf/eq.c new file mode 100644 index 0000000000..30c6befd06 --- /dev/null +++ b/gmp/mpf/eq.c @@ -0,0 +1,150 @@ +/* mpf_eq -- Compare two floats up to a specified bit #. + +Copyright 1993, 1995, 1996, 2001, 2002, 2008, 2009, 2012 Free Software +Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +int +mpf_eq (mpf_srcptr u, mpf_srcptr v, mp_bitcnt_t n_bits) +{ + mp_srcptr up, vp, p; + mp_size_t usize, vsize, minsize, maxsize, n_limbs, i, size; + mp_exp_t uexp, vexp; + mp_limb_t diff; + int cnt; + + uexp = u->_mp_exp; + vexp = v->_mp_exp; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* 1. Are the signs different? */ + if ((usize ^ vsize) >= 0) + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + return vsize == 0; + if (vsize == 0) + return 0; + + /* Fall out. */ + } + else + { + /* Either U or V is negative, but not both. */ + return 0; + } + + /* U and V have the same sign and are both non-zero. */ + + /* 2. Are the exponents different? */ + if (uexp != vexp) + return 0; + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + + up += usize; /* point just above most significant limb */ + vp += vsize; /* point just above most significant limb */ + + count_leading_zeros (cnt, up[-1]); + if ((vp[-1] >> (GMP_LIMB_BITS - 1 - cnt)) != 1) + return 0; /* msb positions different */ + + n_bits += cnt - GMP_NAIL_BITS; + n_limbs = (n_bits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; + + usize = MIN (usize, n_limbs); + vsize = MIN (vsize, n_limbs); + +#if 0 + /* Ignore zeros at the low end of U and V. */ + while (up[0] == 0) + up++, usize--; + while (vp[0] == 0) + vp++, vsize--; +#endif + + minsize = MIN (usize, vsize); + maxsize = usize + vsize - minsize; + + up -= minsize; /* point at most significant common limb */ + vp -= minsize; /* point at most significant common limb */ + + /* Compare the most significant part which has explicit limbs for U and V. */ + for (i = minsize - 1; i > 0; i--) + { + if (up[i] != vp[i]) + return 0; + } + + n_bits -= (maxsize - 1) * GMP_NUMB_BITS; + + size = maxsize - minsize; + if (size != 0) + { + if (up[0] != vp[0]) + return 0; + + /* Now either U or V has its limbs consumed, i.e, continues with an + infinite number of implicit zero limbs. Check that the other operand + has just zeros in the corresponding, relevant part. */ + + if (usize > vsize) + p = up - size; + else + p = vp - size; + + for (i = size - 1; i > 0; i--) + { + if (p[i] != 0) + return 0; + } + + diff = p[0]; + } + else + { + /* Both U or V has its limbs consumed. */ + + diff = up[0] ^ vp[0]; + } + + if (n_bits < GMP_NUMB_BITS) + diff >>= GMP_NUMB_BITS - n_bits; + + return diff == 0; +} diff --git a/gmp/mpf/fits_s.h b/gmp/mpf/fits_s.h new file mode 100644 index 0000000000..ec2635f656 --- /dev/null +++ b/gmp/mpf/fits_s.h @@ -0,0 +1,75 @@ +/* mpf_fits_s*_p -- test whether an mpf fits a C signed type. + +Copyright 2001, 2002 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Notice this is equivalent to mpz_set_f + mpz_fits_s*_p. */ + +int +FUNCTION (mpf_srcptr f) __GMP_NOTHROW +{ + mp_size_t fs, fn; + mp_srcptr fp; + mp_exp_t exp; + mp_limb_t fl; + + fs = SIZ(f); + if (fs == 0) + return 1; /* zero fits */ + + exp = EXP(f); + if (exp < 1) + return 1; /* -1 < f < 1 truncates to zero, so fits */ + + fp = PTR(f); + fn = ABS (fs); + + if (exp == 1) + { + fl = fp[fn-1]; + } +#if GMP_NAIL_BITS != 0 + else if (exp == 2 && MAXIMUM > GMP_NUMB_MAX) + { + fl = fp[fn-1]; + if ((fl >> GMP_NAIL_BITS) != 0) + return 0; + fl = (fl << GMP_NUMB_BITS); + if (fn >= 2) + fl |= fp[fn-2]; + } +#endif + else + return 0; + + return fl <= (fs >= 0 ? (mp_limb_t) MAXIMUM : - (mp_limb_t) MINIMUM); +} diff --git a/gmp/mpf/fits_sint.c b/gmp/mpf/fits_sint.c new file mode 100644 index 0000000000..26ace07c38 --- /dev/null +++ b/gmp/mpf/fits_sint.c @@ -0,0 +1,36 @@ +/* mpf_fits_sint_p -- test whether an mpf fits an int. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_sint_p +#define MAXIMUM INT_MAX +#define MINIMUM INT_MIN + +#include "fits_s.h" diff --git a/gmp/mpf/fits_slong.c b/gmp/mpf/fits_slong.c new file mode 100644 index 0000000000..25db68c47d --- /dev/null +++ b/gmp/mpf/fits_slong.c @@ -0,0 +1,36 @@ +/* mpf_fits_slong_p -- test whether an mpf fits a long. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_slong_p +#define MAXIMUM LONG_MAX +#define MINIMUM LONG_MIN + +#include "fits_s.h" diff --git a/gmp/mpf/fits_sshort.c b/gmp/mpf/fits_sshort.c new file mode 100644 index 0000000000..3bfc5a4a34 --- /dev/null +++ b/gmp/mpf/fits_sshort.c @@ -0,0 +1,36 @@ +/* mpf_fits_sshort_p -- test whether an mpf fits a short. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_sshort_p +#define MAXIMUM SHRT_MAX +#define MINIMUM SHRT_MIN + +#include "fits_s.h" diff --git a/gmp/mpf/fits_u.h b/gmp/mpf/fits_u.h new file mode 100644 index 0000000000..65ac60e09a --- /dev/null +++ b/gmp/mpf/fits_u.h @@ -0,0 +1,74 @@ +/* mpf_fits_u*_p -- test whether an mpf fits a C unsigned type. + +Copyright 2001, 2002, 2013 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Notice this is equivalent to mpz_set_f + mpz_fits_u*_p. */ + +int +FUNCTION (mpf_srcptr f) __GMP_NOTHROW +{ + mp_size_t fn; + mp_srcptr fp; + mp_exp_t exp; + mp_limb_t fl; + + exp = EXP(f); + if (exp < 1) + return 1; /* -1 < f < 1 truncates to zero, so fits */ + + fn = SIZ(f); + if (fn <= 0) + return fn == 0; /* zero fits, negatives don't */ + + fp = PTR(f); + + if (exp == 1) + { + fl = fp[fn-1]; + } +#if GMP_NAIL_BITS != 0 + else if (exp == 2 && MAXIMUM > GMP_NUMB_MAX) + { + fl = fp[fn-1]; + if ((fl >> GMP_NAIL_BITS) != 0) + return 0; + fl = (fl << GMP_NUMB_BITS); + if (fn >= 2) + fl |= fp[fn-2]; + } +#endif + else + return 0; + + return fl <= MAXIMUM; +} diff --git a/gmp/mpf/fits_uint.c b/gmp/mpf/fits_uint.c new file mode 100644 index 0000000000..4b107b04d8 --- /dev/null +++ b/gmp/mpf/fits_uint.c @@ -0,0 +1,35 @@ +/* mpf_fits_uint_p -- test whether an mpf fits an unsigned int. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_uint_p +#define MAXIMUM UINT_MAX + +#include "fits_u.h" diff --git a/gmp/mpf/fits_ulong.c b/gmp/mpf/fits_ulong.c new file mode 100644 index 0000000000..1db688ce4c --- /dev/null +++ b/gmp/mpf/fits_ulong.c @@ -0,0 +1,35 @@ +/* mpf_fits_ulong_p -- test whether an mpf fits an unsigned long. + +Copyright 2001, 2002 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_ulong_p +#define MAXIMUM ULONG_MAX + +#include "fits_u.h" diff --git a/gmp/mpf/fits_ushort.c b/gmp/mpf/fits_ushort.c new file mode 100644 index 0000000000..76a3fd9d52 --- /dev/null +++ b/gmp/mpf/fits_ushort.c @@ -0,0 +1,35 @@ +/* mpf_fits_ushort_p -- test whether an mpf fits an unsigned short. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + + +#define FUNCTION mpf_fits_ushort_p +#define MAXIMUM USHRT_MAX + +#include "fits_u.h" diff --git a/gmp/mpf/get_d.c b/gmp/mpf/get_d.c new file mode 100644 index 0000000000..8f6f9bbea1 --- /dev/null +++ b/gmp/mpf/get_d.c @@ -0,0 +1,47 @@ +/* double mpf_get_d (mpf_t src) -- return SRC truncated to a double. + +Copyright 1996, 2001-2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +double +mpf_get_d (mpf_srcptr src) +{ + mp_size_t size, abs_size; + long exp; + + size = SIZ (src); + if (UNLIKELY (size == 0)) + return 0.0; + + abs_size = ABS (size); + exp = (EXP (src) - abs_size) * GMP_NUMB_BITS; + return mpn_get_d (PTR (src), abs_size, size, exp); +} diff --git a/gmp/mpf/get_d_2exp.c b/gmp/mpf/get_d_2exp.c new file mode 100644 index 0000000000..17ce229f25 --- /dev/null +++ b/gmp/mpf/get_d_2exp.c @@ -0,0 +1,61 @@ +/* double mpf_get_d_2exp (signed long int *exp, mpf_t src). + +Copyright 2001-2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +double +mpf_get_d_2exp (signed long int *exp2, mpf_srcptr src) +{ + mp_size_t size, abs_size; + mp_srcptr ptr; + int cnt; + long exp; + + size = SIZ(src); + if (UNLIKELY (size == 0)) + { + *exp2 = 0; + return 0.0; + } + + ptr = PTR(src); + abs_size = ABS (size); + count_leading_zeros (cnt, ptr[abs_size - 1]); + cnt -= GMP_NAIL_BITS; + + exp = EXP(src) * GMP_NUMB_BITS - cnt; + *exp2 = exp; + + return mpn_get_d (ptr, abs_size, (mp_size_t) 0, + (long) - (abs_size * GMP_NUMB_BITS - cnt)); +} diff --git a/gmp/mpf/get_dfl_prec.c b/gmp/mpf/get_dfl_prec.c new file mode 100644 index 0000000000..9a773d83ac --- /dev/null +++ b/gmp/mpf/get_dfl_prec.c @@ -0,0 +1,39 @@ +/* mpf_get_default_prec -- return default precision in bits. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +mp_bitcnt_t +mpf_get_default_prec (void) __GMP_NOTHROW +{ + return __GMPF_PREC_TO_BITS (__gmp_default_fp_limb_precision); +} diff --git a/gmp/mpf/get_prc.c b/gmp/mpf/get_prc.c new file mode 100644 index 0000000000..3b3283a48c --- /dev/null +++ b/gmp/mpf/get_prc.c @@ -0,0 +1,38 @@ +/* mpf_get_prec(x) -- Return the precision in bits of x. + +Copyright 1996, 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +mp_bitcnt_t +mpf_get_prec (mpf_srcptr x) __GMP_NOTHROW +{ + return __GMPF_PREC_TO_BITS (x->_mp_prec); +} diff --git a/gmp/mpf/get_si.c b/gmp/mpf/get_si.c new file mode 100644 index 0000000000..5b63dbd425 --- /dev/null +++ b/gmp/mpf/get_si.c @@ -0,0 +1,87 @@ +/* mpf_get_si -- mpf to long conversion + +Copyright 2001, 2002, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Any fraction bits are truncated, meaning simply discarded. + + For values bigger than a long, the low bits are returned, like + mpz_get_si, but this isn't documented. + + Notice this is equivalent to mpz_set_f + mpz_get_si. + + + Implementation: + + fl is established in basically the same way as for mpf_get_ui, see that + code for explanations of the conditions. + + However unlike mpf_get_ui we need an explicit return 0 for exp<=0. When + f is a negative fraction (ie. size<0 and exp<=0) we can't let fl==0 go + through to the zany final "~ ((fl - 1) & LONG_MAX)", that would give + -0x80000000 instead of the desired 0. */ + +long +mpf_get_si (mpf_srcptr f) __GMP_NOTHROW +{ + mp_exp_t exp; + mp_size_t size, abs_size; + mp_srcptr fp; + mp_limb_t fl; + + exp = EXP (f); + size = SIZ (f); + fp = PTR (f); + + /* fraction alone truncates to zero + this also covers zero, since we have exp==0 for zero */ + if (exp <= 0) + return 0L; + + /* there are some limbs above the radix point */ + + fl = 0; + abs_size = ABS (size); + if (abs_size >= exp) + fl = fp[abs_size-exp]; + +#if BITS_PER_ULONG > GMP_NUMB_BITS + if (exp > 1 && abs_size+1 >= exp) + fl |= fp[abs_size - exp + 1] << GMP_NUMB_BITS; +#endif + + if (size > 0) + return fl & LONG_MAX; + else + /* this form necessary to correctly handle -0x80..00 */ + return -1 - (long) ((fl - 1) & LONG_MAX); +} diff --git a/gmp/mpf/get_str.c b/gmp/mpf/get_str.c new file mode 100644 index 0000000000..98af03272e --- /dev/null +++ b/gmp/mpf/get_str.c @@ -0,0 +1,330 @@ +/* mpf_get_str (digit_ptr, exp, base, n_digits, a) -- Convert the floating + point number A to a base BASE number and store N_DIGITS raw digits at + DIGIT_PTR, and the base BASE exponent in the word pointed to by EXP. For + example, the number 3.1416 would be returned as "31416" in DIGIT_PTR and + 1 in EXP. + +Copyright 1993-1997, 2000-2003, 2005, 2006, 2011 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdlib.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" /* for count_leading_zeros */ + +/* Could use some more work. + + 1. Allocation is excessive. Try to combine areas. Perhaps use result + string area for temp limb space? + 2. We generate up to two limbs of extra digits. This is because we don't + check the exact number of bits in the input operand, and from that + compute an accurate exponent (variable e in the code). It would be + cleaner and probably somewhat faster to change this. +*/ + +/* Compute base^exp and return the most significant prec limbs in rp[]. + Put the count of omitted low limbs in *ign. + Return the actual size (which might be less than prec). + Allocation of rp[] and the temporary tp[] should be 2*prec+2 limbs. */ +static mp_size_t +mpn_pow_1_highpart (mp_ptr rp, mp_size_t *ignp, + mp_limb_t base, unsigned long exp, + mp_size_t prec, mp_ptr tp) +{ + mp_size_t ign; /* counts number of ignored low limbs in r */ + mp_size_t off; /* keeps track of offset where value starts */ + mp_ptr passed_rp = rp; + mp_size_t rn; + int cnt; + int i; + + if (exp == 0) + { + rp[0] = 1; + *ignp = 0; + return 1; + } + + rp[0] = base; + rn = 1; + off = 0; + ign = 0; + count_leading_zeros (cnt, exp); + for (i = GMP_LIMB_BITS - cnt - 2; i >= 0; i--) + { + mpn_sqr (tp, rp + off, rn); + rn = 2 * rn; + rn -= tp[rn - 1] == 0; + ign <<= 1; + + off = 0; + if (rn > prec) + { + ign += rn - prec; + off = rn - prec; + rn = prec; + } + MP_PTR_SWAP (rp, tp); + + if (((exp >> i) & 1) != 0) + { + mp_limb_t cy; + cy = mpn_mul_1 (rp, rp + off, rn, base); + rp[rn] = cy; + rn += cy != 0; + off = 0; + } + } + + if (rn > prec) + { + ASSERT (rn == prec + 1); + + ign += rn - prec; + rp += rn - prec; + rn = prec; + } + + /* With somewhat less than 50% probability, we can skip this copy. */ + if (passed_rp != rp + off) + MPN_COPY_INCR (passed_rp, rp + off, rn); + *ignp = ign; + return rn; +} + +char * +mpf_get_str (char *dbuf, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u) +{ + mp_exp_t ue; + mp_size_t n_limbs_needed; + size_t max_digits; + mp_ptr up, pp, tp; + mp_size_t un, pn, tn; + unsigned char *tstr; + mp_exp_t exp_in_base; + size_t n_digits_computed; + mp_size_t i; + const char *num_to_text; + size_t alloc_size = 0; + char *dp; + TMP_DECL; + + up = PTR(u); + un = ABSIZ(u); + ue = EXP(u); + + if (base >= 0) + { + num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz"; + if (base <= 1) + base = 10; + else if (base > 36) + { + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz"; + if (base > 62) + return NULL; + } + } + else + { + base = -base; + if (base <= 1) + base = 10; + else if (base > 36) + return NULL; + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + + MPF_SIGNIFICANT_DIGITS (max_digits, base, PREC(u)); + if (n_digits == 0 || n_digits > max_digits) + n_digits = max_digits; + + if (dbuf == 0) + { + /* We didn't get a string from the user. Allocate one (and return + a pointer to it) with space for `-' and terminating null. */ + alloc_size = n_digits + 2; + dbuf = (char *) (*__gmp_allocate_func) (n_digits + 2); + } + + if (un == 0) + { + *exp = 0; + *dbuf = 0; + n_digits = 0; + goto done; + } + + TMP_MARK; + + /* Allocate temporary digit space. We can't put digits directly in the user + area, since we generate more digits than requested. (We allocate + 2 * GMP_LIMB_BITS extra bytes because of the digit block nature of the + conversion.) */ + tstr = (unsigned char *) TMP_ALLOC (n_digits + 2 * GMP_LIMB_BITS + 3); + + LIMBS_PER_DIGIT_IN_BASE (n_limbs_needed, n_digits, base); + + if (ue <= n_limbs_needed) + { + /* We need to multiply number by base^n to get an n_digits integer part. */ + mp_size_t n_more_limbs_needed, ign, off; + unsigned long e; + + n_more_limbs_needed = n_limbs_needed - ue; + DIGITS_IN_BASE_PER_LIMB (e, n_more_limbs_needed, base); + + if (un > n_limbs_needed) + { + up += un - n_limbs_needed; + un = n_limbs_needed; + } + pp = TMP_ALLOC_LIMBS (2 * n_limbs_needed + 2); + tp = TMP_ALLOC_LIMBS (2 * n_limbs_needed + 2); + + pn = mpn_pow_1_highpart (pp, &ign, (mp_limb_t) base, e, n_limbs_needed, tp); + if (un > pn) + mpn_mul (tp, up, un, pp, pn); /* FIXME: mpn_mul_highpart */ + else + mpn_mul (tp, pp, pn, up, un); /* FIXME: mpn_mul_highpart */ + tn = un + pn; + tn -= tp[tn - 1] == 0; + off = un - ue - ign; + if (off < 0) + { + MPN_COPY_DECR (tp - off, tp, tn); + MPN_ZERO (tp, -off); + tn -= off; + off = 0; + } + n_digits_computed = mpn_get_str (tstr, base, tp + off, tn - off); + + exp_in_base = n_digits_computed - e; + } + else + { + /* We need to divide number by base^n to get an n_digits integer part. */ + mp_size_t n_less_limbs_needed, ign, off, xn; + unsigned long e; + mp_ptr dummyp, xp; + + n_less_limbs_needed = ue - n_limbs_needed; + DIGITS_IN_BASE_PER_LIMB (e, n_less_limbs_needed, base); + + if (un > n_limbs_needed) + { + up += un - n_limbs_needed; + un = n_limbs_needed; + } + pp = TMP_ALLOC_LIMBS (2 * n_limbs_needed + 2); + tp = TMP_ALLOC_LIMBS (2 * n_limbs_needed + 2); + + pn = mpn_pow_1_highpart (pp, &ign, (mp_limb_t) base, e, n_limbs_needed, tp); + + xn = n_limbs_needed + (n_less_limbs_needed-ign); + xp = TMP_ALLOC_LIMBS (xn); + off = xn - un; + MPN_ZERO (xp, off); + MPN_COPY (xp + off, up, un); + + dummyp = TMP_ALLOC_LIMBS (pn); + mpn_tdiv_qr (tp, dummyp, (mp_size_t) 0, xp, xn, pp, pn); + tn = xn - pn + 1; + tn -= tp[tn - 1] == 0; + n_digits_computed = mpn_get_str (tstr, base, tp, tn); + + exp_in_base = n_digits_computed + e; + } + + /* We should normally have computed too many digits. Round the result + at the point indicated by n_digits. */ + if (n_digits_computed > n_digits) + { + size_t i; + /* Round the result. */ + if (tstr[n_digits] * 2 >= base) + { + n_digits_computed = n_digits; + for (i = n_digits - 1;; i--) + { + unsigned int x; + x = ++(tstr[i]); + if (x != base) + break; + n_digits_computed--; + if (i == 0) + { + /* We had something like `bbbbbbb...bd', where 2*d >= base + and `b' denotes digit with significance base - 1. + This rounds up to `1', increasing the exponent. */ + tstr[0] = 1; + n_digits_computed = 1; + exp_in_base++; + break; + } + } + } + } + + /* We might have fewer digits than requested as a result of rounding above, + (i.e. 0.999999 => 1.0) or because we have a number that simply doesn't + need many digits in this base (e.g., 0.125 in base 10). */ + if (n_digits > n_digits_computed) + n_digits = n_digits_computed; + + /* Remove trailing 0. There can be many zeros. */ + while (n_digits != 0 && tstr[n_digits - 1] == 0) + n_digits--; + + dp = dbuf + (SIZ(u) < 0); + + /* Translate to ASCII and copy to result string. */ + for (i = 0; i < n_digits; i++) + dp[i] = num_to_text[tstr[i]]; + dp[n_digits] = 0; + + *exp = exp_in_base; + + if (SIZ(u) < 0) + { + dbuf[0] = '-'; + n_digits++; + } + + TMP_FREE; + + done: + /* If the string was alloced then resize it down to the actual space + required. */ + if (alloc_size != 0) + { + __GMP_REALLOCATE_FUNC_MAYBE_TYPE (dbuf, alloc_size, n_digits + 1, char); + } + + return dbuf; +} diff --git a/gmp/mpf/get_ui.c b/gmp/mpf/get_ui.c new file mode 100644 index 0000000000..eb9b30e69c --- /dev/null +++ b/gmp/mpf/get_ui.c @@ -0,0 +1,102 @@ +/* mpf_get_ui -- mpf to ulong conversion + +Copyright 2001, 2002, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Any fraction bits are truncated, meaning simply discarded. + + For values bigger than a ulong, the low bits are returned (the low + absolute value bits actually), like mpz_get_ui, but this isn't + documented. + + Notice this is equivalent to mpz_set_f + mpz_get_ui. + + + Implementation: + + The limb just above the radix point for us to extract is ptr[size-exp]. + + We need to check that the size-exp index falls in our available data + range, 0 to size-1 inclusive. We test this without risk of an overflow + involving exp by requiring size>=exp (giving size-exp >= 0) and exp>0 + (giving size-exp <= size-1). + + Notice if size==0 there's no fetch, since of course size>=exp and exp>0 + can only be true if size>0. So there's no special handling for size==0, + it comes out as 0 the same as any other time we have no data at our + target index. + + For nails, the second limb above the radix point is also required, this + is ptr[size-exp+1]. + + Again we need to check that size-exp+1 falls in our data range, 0 to + size-1 inclusive. We test without risk of overflow by requiring + size+1>=exp (giving size-exp+1 >= 0) and exp>1 (giving size-exp+1 <= + size-1). + + And again if size==0 these second fetch conditions are not satisfied + either since size+1>=exp and exp>1 are only true if size>0. + + The code is arranged with exp>0 wrapping the exp>1 test since exp>1 is + mis-compiled by alpha gcc prior to version 3.4. It re-writes it as + exp-1>0, which is incorrect when exp==MP_EXP_T_MIN. By having exp>0 + tested first we ensure MP_EXP_T_MIN doesn't reach exp>1. */ + +unsigned long +mpf_get_ui (mpf_srcptr f) __GMP_NOTHROW +{ + mp_size_t size; + mp_exp_t exp; + mp_srcptr fp; + mp_limb_t fl; + + exp = EXP (f); + size = SIZ (f); + fp = PTR (f); + + fl = 0; + if (exp > 0) + { + /* there are some limbs above the radix point */ + + size = ABS (size); + if (size >= exp) + fl = fp[size-exp]; + +#if BITS_PER_ULONG > GMP_NUMB_BITS + if (exp > 1 && size+1 >= exp) + fl += (fp[size-exp+1] << GMP_NUMB_BITS); +#endif + } + + return (unsigned long) fl; +} diff --git a/gmp/mpf/init.c b/gmp/mpf/init.c new file mode 100644 index 0000000000..d8590f27a7 --- /dev/null +++ b/gmp/mpf/init.c @@ -0,0 +1,42 @@ +/* mpf_init() -- Make a new multiple precision number with value 0. + +Copyright 1993-1995, 2000, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init (mpf_ptr r) +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_size = 0; + r->_mp_exp = 0; + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); +} diff --git a/gmp/mpf/init2.c b/gmp/mpf/init2.c new file mode 100644 index 0000000000..a7891e33dd --- /dev/null +++ b/gmp/mpf/init2.c @@ -0,0 +1,44 @@ +/* mpf_init2() -- Make a new multiple precision number with value 0. + +Copyright 1993-1995, 2000, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init2 (mpf_ptr r, mp_bitcnt_t prec_in_bits) +{ + mp_size_t prec; + + prec = __GMPF_BITS_TO_PREC (prec_in_bits); + r->_mp_size = 0; + r->_mp_exp = 0; + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); +} diff --git a/gmp/mpf/inits.c b/gmp/mpf/inits.c new file mode 100644 index 0000000000..fb14c6b0fd --- /dev/null +++ b/gmp/mpf/inits.c @@ -0,0 +1,49 @@ +/* mpf_inits() -- Initialize multiple mpf_t variables and set them to 0. + +Copyright 2009 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdarg.h> +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_inits (mpf_ptr x, ...) +{ + va_list ap; + + va_start (ap, x); + + while (x != NULL) + { + mpf_init (x); + x = va_arg (ap, mpf_ptr); + } + va_end (ap); +} diff --git a/gmp/mpf/inp_str.c b/gmp/mpf/inp_str.c new file mode 100644 index 0000000000..45cc34cebb --- /dev/null +++ b/gmp/mpf/inp_str.c @@ -0,0 +1,93 @@ +/* mpf_inp_str(dest_float, stream, base) -- Input a number in base + BASE from stdio stream STREAM and store the result in DEST_FLOAT. + +Copyright 1996, 2000-2002, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> +#include <ctype.h> +#include "gmp.h" +#include "gmp-impl.h" + +size_t +mpf_inp_str (mpf_ptr rop, FILE *stream, int base) +{ + char *str; + size_t alloc_size, str_size; + int c; + int res; + size_t nread; + + if (stream == 0) + stream = stdin; + + alloc_size = 100; + str = (char *) (*__gmp_allocate_func) (alloc_size); + str_size = 0; + nread = 0; + + /* Skip whitespace. */ + do + { + c = getc (stream); + nread++; + } + while (isspace (c)); + + for (;;) + { + if (str_size >= alloc_size) + { + size_t old_alloc_size = alloc_size; + alloc_size = alloc_size * 3 / 2; + str = (char *) (*__gmp_reallocate_func) (str, old_alloc_size, alloc_size); + } + if (c == EOF || isspace (c)) + break; + str[str_size++] = c; + c = getc (stream); + } + ungetc (c, stream); + nread--; + + if (str_size >= alloc_size) + { + size_t old_alloc_size = alloc_size; + alloc_size = alloc_size * 3 / 2; + str = (char *) (*__gmp_reallocate_func) (str, old_alloc_size, alloc_size); + } + str[str_size] = 0; + + res = mpf_set_str (rop, str, base); + (*__gmp_free_func) (str, alloc_size); + + if (res == -1) + return 0; /* error */ + + return str_size + nread; +} diff --git a/gmp/mpf/int_p.c b/gmp/mpf/int_p.c new file mode 100644 index 0000000000..91e62266be --- /dev/null +++ b/gmp/mpf/int_p.c @@ -0,0 +1,59 @@ +/* mpf_integer_p -- test whether an mpf is an integer */ + +/* +Copyright 2001, 2002 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +int +mpf_integer_p (mpf_srcptr f) __GMP_NOTHROW +{ + mp_srcptr ptr; + mp_exp_t exp; + mp_size_t size, frac, i; + + size = SIZ (f); + if (size == 0) + return 1; /* zero is an integer */ + + exp = EXP (f); + if (exp <= 0) + return 0; /* has only fraction limbs */ + + /* any fraction limbs must be zero */ + frac = ABS (size) - exp; + ptr = PTR (f); + for (i = 0; i < frac; i++) + if (ptr[i] != 0) + return 0; + + return 1; +} diff --git a/gmp/mpf/iset.c b/gmp/mpf/iset.c new file mode 100644 index 0000000000..c8c35e50e9 --- /dev/null +++ b/gmp/mpf/iset.c @@ -0,0 +1,62 @@ +/* mpf_init_set -- Initialize a float and assign it from another float. + +Copyright 1993-1995, 2000, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init_set (mpf_ptr r, mpf_srcptr s) +{ + mp_ptr rp, sp; + mp_size_t ssize, size; + mp_size_t prec; + + prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); + r->_mp_prec = prec; + + prec++; /* lie not to lose precision in assignment */ + ssize = s->_mp_size; + size = ABS (ssize); + + rp = r->_mp_d; + sp = s->_mp_d; + + if (size > prec) + { + sp += size - prec; + size = prec; + } + + r->_mp_exp = s->_mp_exp; + r->_mp_size = ssize >= 0 ? size : -size; + + MPN_COPY (rp, sp, size); +} diff --git a/gmp/mpf/iset_d.c b/gmp/mpf/iset_d.c new file mode 100644 index 0000000000..d128db9b4d --- /dev/null +++ b/gmp/mpf/iset_d.c @@ -0,0 +1,42 @@ +/* mpf_init_set_d -- Initialize a float and assign it from a double. + +Copyright 1993-1995, 2000, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init_set_d (mpf_ptr r, double val) +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); + + mpf_set_d (r, val); +} diff --git a/gmp/mpf/iset_si.c b/gmp/mpf/iset_si.c new file mode 100644 index 0000000000..f7e9005086 --- /dev/null +++ b/gmp/mpf/iset_si.c @@ -0,0 +1,58 @@ +/* mpf_init_set_si() -- Initialize a float and assign it from a signed int. + +Copyright 1993-1995, 2000, 2001, 2003, 2004, 2012 Free Software Foundation, +Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init_set_si (mpf_ptr r, long int val) +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + mp_size_t size; + mp_limb_t vl; + + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); + + vl = (mp_limb_t) ABS_CAST (unsigned long int, val); + + r->_mp_d[0] = vl & GMP_NUMB_MASK; + size = vl != 0; + +#if BITS_PER_ULONG > GMP_NUMB_BITS + vl >>= GMP_NUMB_BITS; + r->_mp_d[1] = vl; + size += (vl != 0); +#endif + + r->_mp_exp = size; + r->_mp_size = val >= 0 ? size : -size; +} diff --git a/gmp/mpf/iset_str.c b/gmp/mpf/iset_str.c new file mode 100644 index 0000000000..a181f80c93 --- /dev/null +++ b/gmp/mpf/iset_str.c @@ -0,0 +1,44 @@ +/* mpf_init_set_str -- Initialize a float and assign it from a string. + +Copyright 1995, 1996, 2000, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +mpf_init_set_str (mpf_ptr r, const char *s, int base) +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_size = 0; + r->_mp_exp = 0; + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); + + return mpf_set_str (r, s, base); +} diff --git a/gmp/mpf/iset_ui.c b/gmp/mpf/iset_ui.c new file mode 100644 index 0000000000..f047982116 --- /dev/null +++ b/gmp/mpf/iset_ui.c @@ -0,0 +1,53 @@ +/* mpf_init_set_ui() -- Initialize a float and assign it from an unsigned int. + +Copyright 1993-1995, 2000, 2001, 2003, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_init_set_ui (mpf_ptr r, unsigned long int val) +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + mp_size_t size; + + r->_mp_prec = prec; + r->_mp_d = (mp_ptr) (*__gmp_allocate_func) ((size_t) (prec + 1) * GMP_LIMB_BYTES); + r->_mp_d[0] = val & GMP_NUMB_MASK; + size = (val != 0); + +#if BITS_PER_ULONG > GMP_NUMB_BITS + val >>= GMP_NUMB_BITS; + r->_mp_d[1] = val; + size += (val != 0); +#endif + + r->_mp_size = size; + r->_mp_exp = size; +} diff --git a/gmp/mpf/mul.c b/gmp/mpf/mul.c new file mode 100644 index 0000000000..41d1db7c14 --- /dev/null +++ b/gmp/mpf/mul.c @@ -0,0 +1,96 @@ +/* mpf_mul -- Multiply two floats. + +Copyright 1993, 1994, 1996, 2001, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_mul (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +{ + mp_srcptr up, vp; + mp_size_t usize, vsize; + mp_size_t sign_product; + mp_size_t prec = r->_mp_prec; + TMP_DECL; + + TMP_MARK; + usize = u->_mp_size; + vsize = v->_mp_size; + sign_product = usize ^ vsize; + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + + if (usize == 0 || vsize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; /* ??? */ + } + else + { + mp_size_t rsize; + mp_limb_t cy_limb; + mp_ptr rp, tp; + mp_size_t adj; + + rsize = usize + vsize; + tp = TMP_ALLOC_LIMBS (rsize); + cy_limb = (usize >= vsize + ? mpn_mul (tp, up, usize, vp, vsize) + : mpn_mul (tp, vp, vsize, up, usize)); + + adj = cy_limb == 0; + rsize -= adj; + prec++; + if (rsize > prec) + { + tp += rsize - prec; + rsize = prec; + } + rp = r->_mp_d; + MPN_COPY (rp, tp, rsize); + r->_mp_exp = u->_mp_exp + v->_mp_exp - adj; + r->_mp_size = sign_product >= 0 ? rsize : -rsize; + } + TMP_FREE; +} diff --git a/gmp/mpf/mul_2exp.c b/gmp/mpf/mul_2exp.c new file mode 100644 index 0000000000..83df2176d3 --- /dev/null +++ b/gmp/mpf/mul_2exp.c @@ -0,0 +1,133 @@ +/* mpf_mul_2exp -- Multiply a float by 2^n. + +Copyright 1993, 1994, 1996, 2000-2002, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Multiples of GMP_NUMB_BITS in exp simply mean an amount added to EXP(u) + to set EXP(r). The remainder exp%GMP_NUMB_BITS is then a left shift for + the limb data. + + If exp%GMP_NUMB_BITS == 0 then there's no shifting, we effectively just + do an mpz_set with changed EXP(r). Like mpz_set we take prec+1 limbs in + this case. Although just prec would suffice, it's nice to have + mpf_mul_2exp with exp==0 come out the same as mpz_set. + + When shifting we take up to prec many limbs from the input. Our shift is + cy = mpn_lshift (PTR(r), PTR(u)+k, size, ...), where k is the number of + low limbs dropped from u, and the carry out is stored to PTR(r)[size]. + + It may be noted that the low limb PTR(r)[0] doesn't incorporate bits from + PTR(u)[k-1] (when k>=1 makes that limb available). Taking just prec + limbs from the input (with the high non-zero) is enough bits for the + application requested precision, there's no need for extra work. + + If r==u the shift will have overlapping operands. When k==0 (ie. when + usize <= prec), the overlap is supported by lshift (ie. dst == src). + + But when r==u and k>=1 (ie. usize > prec), we would have an invalid + overlap (ie. mpn_lshift (rp, rp+k, ...)). In this case we must instead + use mpn_rshift (PTR(r)+1, PTR(u)+k, size, NUMB-shift) with the carry out + stored to PTR(r)[0]. An rshift by NUMB-shift bits like this gives + identical data, it's just its overlap restrictions which differ. + + Enhancements: + + The way mpn_lshift is used means successive mpf_mul_2exp calls on the + same operand will accumulate low zero limbs, until prec+1 limbs is + reached. This is wasteful for subsequent operations. When abs_usize <= + prec, we should test the low exp%GMP_NUMB_BITS many bits of PTR(u)[0], + ie. those which would be shifted out by an mpn_rshift. If they're zero + then use that mpn_rshift. */ + +void +mpf_mul_2exp (mpf_ptr r, mpf_srcptr u, mp_bitcnt_t exp) +{ + mp_srcptr up; + mp_ptr rp = r->_mp_d; + mp_size_t usize; + mp_size_t abs_usize; + mp_size_t prec = r->_mp_prec; + mp_exp_t uexp = u->_mp_exp; + + usize = u->_mp_size; + + if (UNLIKELY (usize == 0)) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + abs_usize = ABS (usize); + up = u->_mp_d; + + if (exp % GMP_NUMB_BITS == 0) + { + prec++; /* retain more precision here as we don't need + to account for carry-out here */ + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + } + if (rp != up) + MPN_COPY_INCR (rp, up, abs_usize); + r->_mp_exp = uexp + exp / GMP_NUMB_BITS; + } + else + { + mp_limb_t cy_limb; + mp_size_t adj; + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + /* Use mpn_rshift since mpn_lshift operates downwards, and we + therefore would clobber part of U before using that part, in case + R is the same variable as U. */ + cy_limb = mpn_rshift (rp + 1, up, abs_usize, + GMP_NUMB_BITS - exp % GMP_NUMB_BITS); + rp[0] = cy_limb; + adj = rp[abs_usize] != 0; + } + else + { + cy_limb = mpn_lshift (rp, up, abs_usize, exp % GMP_NUMB_BITS); + rp[abs_usize] = cy_limb; + adj = cy_limb != 0; + } + + abs_usize += adj; + r->_mp_exp = uexp + exp / GMP_NUMB_BITS + adj; + } + r->_mp_size = usize >= 0 ? abs_usize : -abs_usize; +} diff --git a/gmp/mpf/mul_ui.c b/gmp/mpf/mul_ui.c new file mode 100644 index 0000000000..031b2fe06f --- /dev/null +++ b/gmp/mpf/mul_ui.c @@ -0,0 +1,182 @@ +/* mpf_mul_ui -- Multiply a float and an unsigned integer. + +Copyright 1993, 1994, 1996, 2001, 2003, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* The core operation is a multiply of PREC(r) limbs from u by v, producing + either PREC(r) or PREC(r)+1 result limbs. If u is shorter than PREC(r), + then we take only as much as it has. If u is longer we incorporate a + carry from the lower limbs. + + If u has just 1 extra limb, then the carry to add is high(up[0]*v). That + is of course what mpn_mul_1 would do if it was called with PREC(r)+1 + limbs of input. + + If u has more than 1 extra limb, then there can be a further carry bit + out of lower uncalculated limbs (the way the low of one product adds to + the high of the product below it). This is of course what an mpn_mul_1 + would do if it was called with the full u operand. But we instead work + downwards explicitly, until a carry occurs or until a value other than + GMP_NUMB_MAX occurs (that being the only value a carry bit can propagate + across). + + The carry determination normally requires two umul_ppmm's, only rarely + will GMP_NUMB_MAX occur and require further products. + + The carry limb is conveniently added into the mul_1 using mpn_mul_1c when + that function exists, otherwise a subsequent mpn_add_1 is needed. + + Clearly when mpn_mul_1c is used the carry must be calculated first. But + this is also the case when add_1 is used, since if r==u and ABSIZ(r) > + PREC(r) then the mpn_mul_1 overwrites the low part of the input. + + A reuse r==u with size > prec can occur from a size PREC(r)+1 in the + usual way, or it can occur from an mpf_set_prec_raw leaving a bigger + sized value. In both cases we can end up calling mpn_mul_1 with + overlapping src and dst regions, but this will be with dst < src and such + an overlap is permitted. + + Not done: + + No attempt is made to determine in advance whether the result will be + PREC(r) or PREC(r)+1 limbs. If it's going to be PREC(r)+1 then we could + take one less limb from u and generate just PREC(r), that of course + satisfying application requested precision. But any test counting bits + or forming the high product would almost certainly take longer than the + incremental cost of an extra limb in mpn_mul_1. + + Enhancements: + + Repeated mpf_mul_ui's with an even v will accumulate low zero bits on the + result, leaving low zero limbs after a while, which it might be nice to + strip to save work in subsequent operations. Calculating the low limb + explicitly would let us direct mpn_mul_1 to put the balance at rp when + the low is zero (instead of normally rp+1). But it's not clear whether + this would be worthwhile. Explicit code for the low limb will probably + be slower than having it done in mpn_mul_1, so we need to consider how + often a zero will be stripped and how much that's likely to save + later. */ + +void +mpf_mul_ui (mpf_ptr r, mpf_srcptr u, unsigned long int v) +{ + mp_srcptr up; + mp_size_t usize; + mp_size_t size; + mp_size_t prec, excess; + mp_limb_t cy_limb, vl, cbit, cin; + mp_ptr rp; + + usize = u->_mp_size; + if (UNLIKELY (v == 0) || UNLIKELY (usize == 0)) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + +#if BITS_PER_ULONG > GMP_NUMB_BITS /* avoid warnings about shift amount */ + if (v > GMP_NUMB_MAX) + { + mpf_t vf; + mp_limb_t vp[2]; + vp[0] = v & GMP_NUMB_MASK; + vp[1] = v >> GMP_NUMB_BITS; + PTR(vf) = vp; + SIZ(vf) = 2; + ASSERT_CODE (PREC(vf) = 2); + EXP(vf) = 2; + mpf_mul (r, u, vf); + return; + } +#endif + + size = ABS (usize); + prec = r->_mp_prec; + up = u->_mp_d; + vl = v; + excess = size - prec; + cin = 0; + + if (excess > 0) + { + /* up is bigger than desired rp, shorten it to prec limbs and + determine a carry-in */ + + mp_limb_t vl_shifted = vl << GMP_NAIL_BITS; + mp_limb_t hi, lo, next_lo, sum; + mp_size_t i; + + /* high limb of top product */ + i = excess - 1; + umul_ppmm (cin, lo, up[i], vl_shifted); + + /* and carry bit out of products below that, if any */ + for (;;) + { + i--; + if (i < 0) + break; + + umul_ppmm (hi, next_lo, up[i], vl_shifted); + lo >>= GMP_NAIL_BITS; + ADDC_LIMB (cbit, sum, hi, lo); + cin += cbit; + lo = next_lo; + + /* Continue only if the sum is GMP_NUMB_MAX. GMP_NUMB_MAX is the + only value a carry from below can propagate across. If we've + just seen the carry out (ie. cbit!=0) then sum!=GMP_NUMB_MAX, + so this test stops us for that case too. */ + if (LIKELY (sum != GMP_NUMB_MAX)) + break; + } + + up += excess; + size = prec; + } + + rp = r->_mp_d; +#if HAVE_NATIVE_mpn_mul_1c + cy_limb = mpn_mul_1c (rp, up, size, vl, cin); +#else + cy_limb = mpn_mul_1 (rp, up, size, vl); + __GMPN_ADD_1 (cbit, rp, rp, size, cin); + cy_limb += cbit; +#endif + rp[size] = cy_limb; + cy_limb = cy_limb != 0; + r->_mp_exp = u->_mp_exp + cy_limb; + size += cy_limb; + r->_mp_size = usize >= 0 ? size : -size; +} diff --git a/gmp/mpf/neg.c b/gmp/mpf/neg.c new file mode 100644 index 0000000000..553018f306 --- /dev/null +++ b/gmp/mpf/neg.c @@ -0,0 +1,62 @@ +/* mpf_neg -- Negate a float. + +Copyright 1993-1995, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_neg (mpf_ptr r, mpf_srcptr u) +{ + mp_size_t size; + + size = -u->_mp_size; + if (r != u) + { + mp_size_t prec; + mp_size_t asize; + mp_ptr rp, up; + + prec = r->_mp_prec + 1; /* lie not to lose precision in assignment */ + asize = ABS (size); + rp = r->_mp_d; + up = u->_mp_d; + + if (asize > prec) + { + up += asize - prec; + asize = prec; + } + + MPN_COPY (rp, up, asize); + r->_mp_exp = u->_mp_exp; + size = size >= 0 ? asize : -asize; + } + r->_mp_size = size; +} diff --git a/gmp/mpf/out_str.c b/gmp/mpf/out_str.c new file mode 100644 index 0000000000..200da74806 --- /dev/null +++ b/gmp/mpf/out_str.c @@ -0,0 +1,117 @@ +/* mpf_out_str (stream, base, n_digits, op) -- Print N_DIGITS digits from + the float OP to STREAM in base BASE. Return the number of characters + written, or 0 if an error occurred. + +Copyright 1996, 1997, 2001, 2002, 2005, 2011 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#define _GNU_SOURCE /* for DECIMAL_POINT in langinfo.h */ + +#include "config.h" + +#include <stdio.h> +#include <string.h> + +#if HAVE_LANGINFO_H +#include <langinfo.h> /* for nl_langinfo */ +#endif + +#if HAVE_LOCALE_H +#include <locale.h> /* for localeconv */ +#endif + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +size_t +mpf_out_str (FILE *stream, int base, size_t n_digits, mpf_srcptr op) +{ + char *str; + mp_exp_t exp; + size_t written; + TMP_DECL; + + TMP_MARK; + + if (base == 0) + base = 10; + if (n_digits == 0) + MPF_SIGNIFICANT_DIGITS (n_digits, base, op->_mp_prec); + + if (stream == 0) + stream = stdout; + + /* Consider these changes: + * Don't allocate memory here for huge n_digits; pass NULL to mpf_get_str. + * Make mpf_get_str allocate extra space when passed NULL, to avoid + allocating two huge string buffers. + * Implement more/other allocation reductions tricks. */ + + str = (char *) TMP_ALLOC (n_digits + 2); /* extra for minus sign and \0 */ + + mpf_get_str (str, &exp, base, n_digits, op); + n_digits = strlen (str); + + written = 0; + + /* Write sign */ + if (str[0] == '-') + { + str++; + fputc ('-', stream); + written = 1; + n_digits--; + } + + { + const char *point = GMP_DECIMAL_POINT; + size_t pointlen = strlen (point); + putc ('0', stream); + fwrite (point, 1, pointlen, stream); + written += pointlen + 1; + } + + /* Write mantissa */ + { + size_t fwret; + fwret = fwrite (str, 1, n_digits, stream); + written += fwret; + } + + /* Write exponent */ + { + int fpret; + fpret = fprintf (stream, (base <= 10 ? "e%ld" : "@%ld"), exp); + written += fpret; + } + + TMP_FREE; + return ferror (stream) ? 0 : written; +} diff --git a/gmp/mpf/pow_ui.c b/gmp/mpf/pow_ui.c new file mode 100644 index 0000000000..a5af431059 --- /dev/null +++ b/gmp/mpf/pow_ui.c @@ -0,0 +1,54 @@ +/* mpf_pow_ui -- Compute b^e. + +Copyright 1998, 1999, 2001, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_pow_ui (mpf_ptr r, mpf_srcptr b, unsigned long int e) +{ + mpf_t b2; + + mpf_init2 (b2, mpf_get_prec (r)); + mpf_set (b2, b); + + if ((e & 1) != 0) + mpf_set (r, b); + else + mpf_set_ui (r, 1); + while (e >>= 1) + { + mpf_mul (b2, b2, b2); + if ((e & 1) != 0) + mpf_mul (r, r, b2); + } + + mpf_clear (b2); +} diff --git a/gmp/mpf/random2.c b/gmp/mpf/random2.c new file mode 100644 index 0000000000..4d7f37e976 --- /dev/null +++ b/gmp/mpf/random2.c @@ -0,0 +1,67 @@ +/* mpf_random2 -- Generate a positive random mpf_t of specified size, with + long runs of consecutive ones and zeros in the binary representation. + Intended for testing of other MP routines. + +Copyright 1995, 1996, 2001-2003 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +void +mpf_random2 (mpf_ptr x, mp_size_t xs, mp_exp_t exp) +{ + mp_size_t xn; + mp_size_t prec; + mp_limb_t elimb; + + xn = ABS (xs); + prec = PREC(x); + + if (xn == 0) + { + EXP(x) = 0; + SIZ(x) = 0; + return; + } + + if (xn > prec + 1) + xn = prec + 1; + + /* General random mantissa. */ + mpn_random2 (PTR(x), xn); + + /* Generate random exponent. */ + _gmp_rand (&elimb, RANDS, GMP_NUMB_BITS); + exp = ABS (exp); + exp = elimb % (2 * exp + 1) - exp; + + EXP(x) = exp; + SIZ(x) = xs < 0 ? -xn : xn; +} diff --git a/gmp/mpf/reldiff.c b/gmp/mpf/reldiff.c new file mode 100644 index 0000000000..f49da08166 --- /dev/null +++ b/gmp/mpf/reldiff.c @@ -0,0 +1,65 @@ +/* mpf_reldiff -- Generate the relative difference of two floats. + +Copyright 1996, 2001, 2004, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* The precision we use for d = x-y is based on what mpf_div will want from + the dividend. It calls mpn_div_q to produce a quotient of rprec+1 limbs. + So rprec+1 == dsize - xsize + 1, hence dprec = rprec+xsize. */ + +void +mpf_reldiff (mpf_t rdiff, mpf_srcptr x, mpf_srcptr y) +{ + if (UNLIKELY (SIZ(x) == 0)) + { + mpf_set_ui (rdiff, (unsigned long int) (mpf_sgn (y) != 0)); + } + else + { + mp_size_t dprec; + mpf_t d; + TMP_DECL; + + TMP_MARK; + dprec = PREC(rdiff) + ABSIZ(x); + ASSERT (PREC(rdiff)+1 == dprec - ABSIZ(x) + 1); + + PREC(d) = dprec; + PTR(d) = TMP_ALLOC_LIMBS (dprec + 1); + + mpf_sub (d, x, y); + SIZ(d) = ABSIZ(d); + mpf_div (rdiff, d, x); + + TMP_FREE; + } +} diff --git a/gmp/mpf/set.c b/gmp/mpf/set.c new file mode 100644 index 0000000000..ec8161d779 --- /dev/null +++ b/gmp/mpf/set.c @@ -0,0 +1,56 @@ +/* mpf_set -- Assign a float from another float. + +Copyright 1993-1995, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set (mpf_ptr r, mpf_srcptr u) +{ + mp_ptr rp, up; + mp_size_t size, asize; + mp_size_t prec; + + prec = r->_mp_prec + 1; /* lie not to lose precision in assignment */ + size = u->_mp_size; + asize = ABS (size); + rp = r->_mp_d; + up = u->_mp_d; + + if (asize > prec) + { + up += asize - prec; + asize = prec; + } + + r->_mp_exp = u->_mp_exp; + r->_mp_size = size >= 0 ? asize : -asize; + MPN_COPY_INCR (rp, up, asize); +} diff --git a/gmp/mpf/set_d.c b/gmp/mpf/set_d.c new file mode 100644 index 0000000000..100194d116 --- /dev/null +++ b/gmp/mpf/set_d.c @@ -0,0 +1,60 @@ +/* mpf_set_d -- Assign a float from a double. + +Copyright 1993-1996, 2001, 2003, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "config.h" + +#if HAVE_FLOAT_H +#include <float.h> /* for DBL_MAX */ +#endif + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set_d (mpf_ptr r, double d) +{ + int negative; + + DOUBLE_NAN_INF_ACTION (d, + __gmp_invalid_operation (), + __gmp_invalid_operation ()); + + if (UNLIKELY (d == 0)) + { + SIZ(r) = 0; + EXP(r) = 0; + return; + } + negative = d < 0; + d = ABS (d); + + SIZ(r) = negative ? -LIMBS_PER_DOUBLE : LIMBS_PER_DOUBLE; + EXP(r) = __gmp_extract_double (PTR(r), d); +} diff --git a/gmp/mpf/set_dfl_prec.c b/gmp/mpf/set_dfl_prec.c new file mode 100644 index 0000000000..04c9a55135 --- /dev/null +++ b/gmp/mpf/set_dfl_prec.c @@ -0,0 +1,40 @@ +/* mpf_set_default_prec -- + +Copyright 1993-1995, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +mp_size_t __gmp_default_fp_limb_precision = __GMPF_BITS_TO_PREC (53); + +void +mpf_set_default_prec (mp_bitcnt_t prec_in_bits) __GMP_NOTHROW +{ + __gmp_default_fp_limb_precision = __GMPF_BITS_TO_PREC (prec_in_bits); +} diff --git a/gmp/mpf/set_prc.c b/gmp/mpf/set_prc.c new file mode 100644 index 0000000000..30ba06c6e6 --- /dev/null +++ b/gmp/mpf/set_prc.c @@ -0,0 +1,69 @@ +/* mpf_set_prec(x) -- Change the precision of x. + +Copyright 1993-1995, 2000, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* A full new_prec+1 limbs are always retained, even though just new_prec + would satisfy the requested precision. If size==new_prec+1 then + certainly new_prec+1 should be kept since no copying is needed in that + case. If just new_prec was kept for size>new_prec+1 it'd be a bit + inconsistent. */ + +void +mpf_set_prec (mpf_ptr x, mp_bitcnt_t new_prec_in_bits) +{ + mp_size_t old_prec, new_prec, new_prec_plus1; + mp_size_t size, sign; + mp_ptr xp; + + new_prec = __GMPF_BITS_TO_PREC (new_prec_in_bits); + old_prec = PREC(x); + + /* do nothing if already the right precision */ + if (new_prec == old_prec) + return; + + PREC(x) = new_prec; + new_prec_plus1 = new_prec + 1; + + /* retain most significant limbs */ + sign = SIZ(x); + size = ABS (sign); + xp = PTR(x); + if (size > new_prec_plus1) + { + SIZ(x) = (sign >= 0 ? new_prec_plus1 : -new_prec_plus1); + MPN_COPY_INCR (xp, xp + size - new_prec_plus1, new_prec_plus1); + } + + PTR(x) = __GMP_REALLOCATE_FUNC_LIMBS (xp, old_prec+1, new_prec_plus1); +} diff --git a/gmp/mpf/set_prc_raw.c b/gmp/mpf/set_prc_raw.c new file mode 100644 index 0000000000..7799442299 --- /dev/null +++ b/gmp/mpf/set_prc_raw.c @@ -0,0 +1,40 @@ +/* mpf_set_prec_raw(x,bits) -- Change the precision of x without changing + allocation. For proper operation, the original precision need to be reset + sooner or later. + +Copyright 1996, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set_prec_raw (mpf_ptr x, mp_bitcnt_t prec_in_bits) __GMP_NOTHROW +{ + x->_mp_prec = __GMPF_BITS_TO_PREC (prec_in_bits); +} diff --git a/gmp/mpf/set_q.c b/gmp/mpf/set_q.c new file mode 100644 index 0000000000..c5739b2abe --- /dev/null +++ b/gmp/mpf/set_q.c @@ -0,0 +1,155 @@ +/* mpf_set_q (mpf_t rop, mpq_t op) -- Convert the rational op to the float rop. + +Copyright 1996, 1999, 2001, 2002, 2004, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +/* As usual the aim is to produce PREC(r) limbs, with the high non-zero. + The basic mpn_tdiv_qr produces a quotient of nsize-dsize+1 limbs, with + either the high or second highest limb non-zero. We arrange for + nsize-dsize+1 to equal prec+1, hence giving either prec or prec+1 result + limbs at PTR(r). + + nsize-dsize+1 == prec+1 is achieved by adjusting num(q), either dropping + low limbs if it's too big, or padding with low zeros if it's too small. + The full given den(q) is always used. + + We cannot truncate den(q), because even when it's much bigger than prec + the last limbs can still influence the final quotient. Often they don't, + but we leave optimization of that to a prospective quotient-only mpn + division. + + Not done: + + If den(q) is a power of 2 then we may end up with low zero limbs on the + result. But nothing is done about this, since it should be unlikely on + random data, and can be left to an application to call mpf_div_2exp if it + might occur with any frequency. + + Enhancements: + + The high quotient limb is non-zero when high{np,dsize} >= {dp,dsize}. We + could make that comparison and use qsize==prec instead of qsize==prec+1, + to save one limb in the division. + + Future: + + If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of + padding n with zeros in temporary space. + + If/when a quotient-only division exists it can be used here immediately. + remp is only to satisfy mpn_tdiv_qr, the remainder is not used. */ + +void +mpf_set_q (mpf_t r, mpq_srcptr q) +{ + mp_srcptr np, dp; + mp_size_t prec, nsize, dsize, qsize, prospective_qsize, tsize, zeros; + mp_size_t sign_quotient, high_zero; + mp_ptr qp, tp, remp; + mp_exp_t exp; + TMP_DECL; + + ASSERT (SIZ(&q->_mp_den) > 0); /* canonical q */ + + nsize = SIZ (&q->_mp_num); + dsize = SIZ (&q->_mp_den); + + if (UNLIKELY (nsize == 0)) + { + SIZ (r) = 0; + EXP (r) = 0; + return; + } + + TMP_MARK; + + prec = PREC (r); + qp = PTR (r); + + sign_quotient = nsize; + nsize = ABS (nsize); + np = PTR (&q->_mp_num); + dp = PTR (&q->_mp_den); + + prospective_qsize = nsize - dsize + 1; /* q from using given n,d sizes */ + exp = prospective_qsize; /* ie. number of integer limbs */ + qsize = prec + 1; /* desired q */ + + zeros = qsize - prospective_qsize; /* n zeros to get desired qsize */ + tsize = nsize + zeros; /* possible copy of n */ + + if (WANT_TMP_DEBUG) + { + /* separate alloc blocks, for malloc debugging */ + remp = TMP_ALLOC_LIMBS (dsize); + tp = NULL; + if (zeros > 0) + tp = TMP_ALLOC_LIMBS (tsize); + } + else + { + /* one alloc with a conditionalized size, for efficiency */ + mp_size_t size = dsize + (zeros > 0 ? tsize : 0); + remp = TMP_ALLOC_LIMBS (size); + tp = remp + dsize; + } + + if (zeros > 0) + { + /* pad n with zeros into temporary space */ + MPN_ZERO (tp, zeros); + MPN_COPY (tp+zeros, np, nsize); + np = tp; + nsize = tsize; + } + else + { + /* shorten n to get desired qsize */ + nsize += zeros; + np -= zeros; + } + + ASSERT (nsize-dsize+1 == qsize); + mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize); + + /* strip possible zero high limb */ + high_zero = (qp[qsize-1] == 0); + qsize -= high_zero; + exp -= high_zero; + + EXP (r) = exp; + SIZ (r) = sign_quotient >= 0 ? qsize : -qsize; + + TMP_FREE; +} diff --git a/gmp/mpf/set_si.c b/gmp/mpf/set_si.c new file mode 100644 index 0000000000..9c47c7511f --- /dev/null +++ b/gmp/mpf/set_si.c @@ -0,0 +1,53 @@ +/* mpf_set_si() -- Assign a float from a signed int. + +Copyright 1993-1995, 2000-2002, 2004, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set_si (mpf_ptr dest, long val) +{ + mp_size_t size; + mp_limb_t vl; + + vl = (mp_limb_t) ABS_CAST (unsigned long int, val); + + dest->_mp_d[0] = vl & GMP_NUMB_MASK; + size = vl != 0; + +#if BITS_PER_ULONG > GMP_NUMB_BITS + vl >>= GMP_NUMB_BITS; + dest->_mp_d[1] = vl; + size += (vl != 0); +#endif + + dest->_mp_exp = size; + dest->_mp_size = val >= 0 ? size : -size; +} diff --git a/gmp/mpf/set_str.c b/gmp/mpf/set_str.c new file mode 100644 index 0000000000..9053accade --- /dev/null +++ b/gmp/mpf/set_str.c @@ -0,0 +1,402 @@ +/* mpf_set_str (dest, string, base) -- Convert the string STRING + in base BASE to a float in dest. If BASE is zero, the leading characters + of STRING is used to figure out the base. + +Copyright 1993-1997, 2000-2003, 2005, 2007, 2008, 2011, 2013 Free Software +Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +/* + This still needs work, as suggested by some FIXME comments. + 1. Don't depend on superfluous mantissa digits. + 2. Allocate temp space more cleverly. + 3. Use mpn_div_q instead of mpn_lshift+mpn_divrem. +*/ + +#define _GNU_SOURCE /* for DECIMAL_POINT in langinfo.h */ + +#include "config.h" + +#include <stdlib.h> +#include <string.h> +#include <ctype.h> + +#if HAVE_LANGINFO_H +#include <langinfo.h> /* for nl_langinfo */ +#endif + +#if HAVE_LOCALE_H +#include <locale.h> /* for localeconv */ +#endif + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +#define digit_value_tab __gmp_digit_value_tab + +/* Compute base^exp and return the most significant prec limbs in rp[]. + Put the count of omitted low limbs in *ign. + Return the actual size (which might be less than prec). */ +static mp_size_t +mpn_pow_1_highpart (mp_ptr rp, mp_size_t *ignp, + mp_limb_t base, mp_exp_t exp, + mp_size_t prec, mp_ptr tp) +{ + mp_size_t ign; /* counts number of ignored low limbs in r */ + mp_size_t off; /* keeps track of offset where value starts */ + mp_ptr passed_rp = rp; + mp_size_t rn; + int cnt; + int i; + + rp[0] = base; + rn = 1; + off = 0; + ign = 0; + count_leading_zeros (cnt, exp); + for (i = GMP_LIMB_BITS - cnt - 2; i >= 0; i--) + { + mpn_sqr (tp, rp + off, rn); + rn = 2 * rn; + rn -= tp[rn - 1] == 0; + ign <<= 1; + + off = 0; + if (rn > prec) + { + ign += rn - prec; + off = rn - prec; + rn = prec; + } + MP_PTR_SWAP (rp, tp); + + if (((exp >> i) & 1) != 0) + { + mp_limb_t cy; + cy = mpn_mul_1 (rp, rp + off, rn, base); + rp[rn] = cy; + rn += cy != 0; + off = 0; + } + } + + if (rn > prec) + { + ign += rn - prec; + rp += rn - prec; + rn = prec; + } + + MPN_COPY_INCR (passed_rp, rp + off, rn); + *ignp = ign; + return rn; +} + +int +mpf_set_str (mpf_ptr x, const char *str, int base) +{ + size_t str_size; + char *s, *begs; + size_t i, j; + int c; + int negative; + char *dotpos = 0; + const char *expptr; + int exp_base; + const char *point = GMP_DECIMAL_POINT; + size_t pointlen = strlen (point); + const unsigned char *digit_value; + TMP_DECL; + + c = (unsigned char) *str; + + /* Skip whitespace. */ + while (isspace (c)) + c = (unsigned char) *++str; + + negative = 0; + if (c == '-') + { + negative = 1; + c = (unsigned char) *++str; + } + + /* Default base to decimal. */ + if (base == 0) + base = 10; + + exp_base = base; + + if (base < 0) + { + exp_base = 10; + base = -base; + } + + digit_value = digit_value_tab; + if (base > 36) + { + /* For bases > 36, use the collating sequence + 0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz. */ + digit_value += 208; + if (base > 62) + return -1; /* too large base */ + } + + /* Require at least one digit, possibly after an initial decimal point. */ + if (digit_value[c] >= base) + { + /* not a digit, must be a decimal point */ + for (i = 0; i < pointlen; i++) + if (str[i] != point[i]) + return -1; + if (digit_value[(unsigned char) str[pointlen]] >= base) + return -1; + } + + /* Locate exponent part of the input. Look from the right of the string, + since the exponent is usually a lot shorter than the mantissa. */ + expptr = NULL; + str_size = strlen (str); + for (i = str_size - 1; i > 0; i--) + { + c = (unsigned char) str[i]; + if (c == '@' || (base <= 10 && (c == 'e' || c == 'E'))) + { + expptr = str + i + 1; + str_size = i; + break; + } + } + + TMP_MARK; + s = begs = (char *) TMP_ALLOC (str_size + 1); + + /* Loop through mantissa, converting it from ASCII to raw byte values. */ + for (i = 0; i < str_size; i++) + { + c = (unsigned char) *str; + if (!isspace (c)) + { + int dig; + + for (j = 0; j < pointlen; j++) + if (str[j] != point[j]) + goto not_point; + if (1) + { + if (dotpos != 0) + { + /* already saw a decimal point, another is invalid */ + TMP_FREE; + return -1; + } + dotpos = s; + str += pointlen - 1; + i += pointlen - 1; + } + else + { + not_point: + dig = digit_value[c]; + if (dig >= base) + { + TMP_FREE; + return -1; + } + *s++ = dig; + } + } + c = (unsigned char) *++str; + } + + str_size = s - begs; + + { + long exp_in_base; + mp_size_t ra, ma, rn, mn; + int cnt; + mp_ptr mp, tp, rp; + mp_exp_t exp_in_limbs; + mp_size_t prec = PREC(x) + 1; + int divflag; + mp_size_t madj, radj; + +#if 0 + size_t n_chars_needed; + + /* This breaks things like 0.000...0001. To safely ignore superfluous + digits, we need to skip over leading zeros. */ + /* Just consider the relevant leading digits of the mantissa. */ + LIMBS_PER_DIGIT_IN_BASE (n_chars_needed, prec, base); + if (str_size > n_chars_needed) + str_size = n_chars_needed; +#endif + + LIMBS_PER_DIGIT_IN_BASE (ma, str_size, base); + mp = TMP_ALLOC_LIMBS (ma); + mn = mpn_set_str (mp, (unsigned char *) begs, str_size, base); + + if (mn == 0) + { + SIZ(x) = 0; + EXP(x) = 0; + TMP_FREE; + return 0; + } + + madj = 0; + /* Ignore excess limbs in MP,MSIZE. */ + if (mn > prec) + { + madj = mn - prec; + mp += mn - prec; + mn = prec; + } + + if (expptr != 0) + { + /* Scan and convert the exponent, in base exp_base. */ + long dig, minus, plusminus; + c = (unsigned char) *expptr; + minus = -(long) (c == '-'); + plusminus = minus | -(long) (c == '+'); + expptr -= plusminus; /* conditional increment */ + c = (unsigned char) *expptr++; + dig = digit_value[c]; + if (dig >= exp_base) + { + TMP_FREE; + return -1; + } + exp_in_base = dig; + c = (unsigned char) *expptr++; + dig = digit_value[c]; + while (dig < exp_base) + { + exp_in_base = exp_in_base * exp_base; + exp_in_base += dig; + c = (unsigned char) *expptr++; + dig = digit_value[c]; + } + exp_in_base = (exp_in_base ^ minus) - minus; /* conditional negation */ + } + else + exp_in_base = 0; + if (dotpos != 0) + exp_in_base -= s - dotpos; + divflag = exp_in_base < 0; + exp_in_base = ABS (exp_in_base); + + if (exp_in_base == 0) + { + MPN_COPY (PTR(x), mp, mn); + SIZ(x) = negative ? -mn : mn; + EXP(x) = mn + madj; + TMP_FREE; + return 0; + } + + ra = 2 * (prec + 1); + rp = TMP_ALLOC_LIMBS (ra); + tp = TMP_ALLOC_LIMBS (ra); + rn = mpn_pow_1_highpart (rp, &radj, (mp_limb_t) base, exp_in_base, prec, tp); + + if (divflag) + { +#if 0 + /* FIXME: Should use mpn_div_q here. */ + ... + mpn_div_q (tp, mp, mn, rp, rn, scratch); + ... +#else + mp_ptr qp; + mp_limb_t qlimb; + if (mn < rn) + { + /* Pad out MP,MSIZE for current divrem semantics. */ + mp_ptr tmp = TMP_ALLOC_LIMBS (rn + 1); + MPN_ZERO (tmp, rn - mn); + MPN_COPY (tmp + rn - mn, mp, mn); + mp = tmp; + madj -= rn - mn; + mn = rn; + } + if ((rp[rn - 1] & GMP_NUMB_HIGHBIT) == 0) + { + mp_limb_t cy; + count_leading_zeros (cnt, rp[rn - 1]); + cnt -= GMP_NAIL_BITS; + mpn_lshift (rp, rp, rn, cnt); + cy = mpn_lshift (mp, mp, mn, cnt); + if (cy) + mp[mn++] = cy; + } + + qp = TMP_ALLOC_LIMBS (prec + 1); + qlimb = mpn_divrem (qp, prec - (mn - rn), mp, mn, rp, rn); + tp = qp; + exp_in_limbs = qlimb + (mn - rn) + (madj - radj); + rn = prec; + if (qlimb != 0) + { + tp[prec] = qlimb; + /* Skip the least significant limb not to overrun the destination + variable. */ + tp++; + } +#endif + } + else + { + tp = TMP_ALLOC_LIMBS (rn + mn); + if (rn > mn) + mpn_mul (tp, rp, rn, mp, mn); + else + mpn_mul (tp, mp, mn, rp, rn); + rn += mn; + rn -= tp[rn - 1] == 0; + exp_in_limbs = rn + madj + radj; + + if (rn > prec) + { + tp += rn - prec; + rn = prec; + exp_in_limbs += 0; + } + } + + MPN_COPY (PTR(x), tp, rn); + SIZ(x) = negative ? -rn : rn; + EXP(x) = exp_in_limbs; + TMP_FREE; + return 0; + } +} diff --git a/gmp/mpf/set_ui.c b/gmp/mpf/set_ui.c new file mode 100644 index 0000000000..617bce13c1 --- /dev/null +++ b/gmp/mpf/set_ui.c @@ -0,0 +1,49 @@ +/* mpf_set_ui() -- Assign a float from an unsigned int. + +Copyright 1993-1995, 2001, 2002, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set_ui (mpf_ptr f, unsigned long val) +{ + mp_size_t size; + + f->_mp_d[0] = val & GMP_NUMB_MASK; + size = val != 0; + +#if BITS_PER_ULONG > GMP_NUMB_BITS + val >>= GMP_NUMB_BITS; + f->_mp_d[1] = val; + size += (val != 0); +#endif + + f->_mp_exp = f->_mp_size = size; +} diff --git a/gmp/mpf/set_z.c b/gmp/mpf/set_z.c new file mode 100644 index 0000000000..fe91904563 --- /dev/null +++ b/gmp/mpf/set_z.c @@ -0,0 +1,57 @@ +/* mpf_set_z -- Assign a float from an integer. + +Copyright 1996, 2001, 2004 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_set_z (mpf_ptr r, mpz_srcptr u) +{ + mp_ptr rp, up; + mp_size_t size, asize; + mp_size_t prec; + + prec = PREC (r) + 1; + size = SIZ (u); + asize = ABS (size); + rp = PTR (r); + up = PTR (u); + + EXP (r) = asize; + + if (asize > prec) + { + up += asize - prec; + asize = prec; + } + + SIZ (r) = size >= 0 ? asize : -asize; + MPN_COPY (rp, up, asize); +} diff --git a/gmp/mpf/size.c b/gmp/mpf/size.c new file mode 100644 index 0000000000..c6b22b64cc --- /dev/null +++ b/gmp/mpf/size.c @@ -0,0 +1,39 @@ +/* mpf_size(x) -- return the number of limbs currently used by the + value of the float X. + +Copyright 1993-1995, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +size_t +mpf_size (mpf_srcptr f) __GMP_NOTHROW +{ + return __GMP_ABS (f->_mp_size); +} diff --git a/gmp/mpf/sqrt.c b/gmp/mpf/sqrt.c new file mode 100644 index 0000000000..44502244e2 --- /dev/null +++ b/gmp/mpf/sqrt.c @@ -0,0 +1,113 @@ +/* mpf_sqrt -- Compute the square root of a float. + +Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005, 2012 Free Software +Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" + + +/* As usual, the aim is to produce PREC(r) limbs of result, with the high + limb non-zero. This is accomplished by applying mpn_sqrtrem to either + 2*prec or 2*prec-1 limbs, both such sizes resulting in prec limbs. + + The choice between 2*prec or 2*prec-1 limbs is based on the input + exponent. With b=2^GMP_NUMB_BITS the limb base then we can think of + effectively taking out a factor b^(2k), for suitable k, to get to an + integer input of the desired size ready for mpn_sqrtrem. It must be an + even power taken out, ie. an even number of limbs, so the square root + gives factor b^k and the radix point is still on a limb boundary. So if + EXP(r) is even we'll get an even number of input limbs 2*prec, or if + EXP(r) is odd we get an odd number 2*prec-1. + + Further limbs below the 2*prec or 2*prec-1 used don't affect the result + and are simply truncated. This can be seen by considering an integer x, + with s=floor(sqrt(x)). s is the unique integer satisfying s^2 <= x < + (s+1)^2. Notice that adding a fraction part to x (ie. some further bits) + doesn't change the inequality, s remains the unique solution. Working + suitable factors of 2 into this argument lets it apply to an intended + precision at any position for any x, not just the integer binary point. + + If the input is smaller than 2*prec or 2*prec-1, then we just pad with + zeros, that of course being our usual interpretation of short inputs. + The effect is to extend the root beyond the size of the input (for + instance into fractional limbs if u is an integer). */ + +void +mpf_sqrt (mpf_ptr r, mpf_srcptr u) +{ + mp_size_t usize; + mp_ptr up, tp; + mp_size_t prec, tsize; + mp_exp_t uexp, expodd; + TMP_DECL; + + usize = u->_mp_size; + if (UNLIKELY (usize <= 0)) + { + if (usize < 0) + SQRT_OF_NEGATIVE; + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK; + + uexp = u->_mp_exp; + prec = r->_mp_prec; + up = u->_mp_d; + + expodd = (uexp & 1); + tsize = 2 * prec - expodd; + r->_mp_size = prec; + r->_mp_exp = (uexp + expodd) / 2; /* ceil(uexp/2) */ + + /* root size is ceil(tsize/2), this will be our desired "prec" limbs */ + ASSERT ((tsize + 1) / 2 == prec); + + tp = TMP_ALLOC_LIMBS (tsize); + + if (usize > tsize) + { + up += usize - tsize; + usize = tsize; + MPN_COPY (tp, up, tsize); + } + else + { + MPN_ZERO (tp, tsize - usize); + MPN_COPY (tp + (tsize - usize), up, usize); + } + + mpn_sqrtrem (r->_mp_d, NULL, tp, tsize); + + TMP_FREE; +} diff --git a/gmp/mpf/sqrt_ui.c b/gmp/mpf/sqrt_ui.c new file mode 100644 index 0000000000..82dec7bcb3 --- /dev/null +++ b/gmp/mpf/sqrt_ui.c @@ -0,0 +1,109 @@ +/* mpf_sqrt_ui -- Compute the square root of an unsigned integer. + +Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005 Free Software Foundation, +Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" + + +/* As usual the aim is to produce PREC(r) limbs of result with the high limb + non-zero. That high limb will end up floor(sqrt(u)), and limbs below are + produced by padding the input with zeros, two for each desired result + limb, being 2*(prec-1) for a total 2*prec-1 limbs passed to mpn_sqrtrem. + The way mpn_sqrtrem calculates floor(sqrt(x)) ensures the root is correct + to the intended accuracy, ie. truncated to prec limbs. + + With nails, u might be two limbs, in which case a total 2*prec limbs is + passed to mpn_sqrtrem (still giving a prec limb result). If uhigh is + zero we adjust back to 2*prec-1, since mpn_sqrtrem requires the high + non-zero. 2*prec limbs are always allocated, even when uhigh is zero, so + the store of uhigh can be done without a conditional. + + u==0 is a special case so the rest of the code can assume the result is + non-zero (ie. will have a non-zero high limb on the result). + + Not done: + + No attempt is made to identify perfect squares. It's considered this can + be left to an application if it might occur with any frequency. As it + stands, mpn_sqrtrem does its normal amount of work on a perfect square + followed by zero limbs, though of course only an mpn_sqrtrem1 would be + actually needed. We also end up leaving our mpf result with lots of low + trailing zeros, slowing down subsequent operations. + + We're not aware of any optimizations that can be made using the fact the + input has lots of trailing zeros (apart from the perfect square + case). */ + + +/* 1 if we (might) need two limbs for u */ +#define U2 (GMP_NUMB_BITS < BITS_PER_ULONG) + +void +mpf_sqrt_ui (mpf_ptr r, unsigned long int u) +{ + mp_size_t rsize, zeros; + mp_ptr tp; + mp_size_t prec; + TMP_DECL; + + if (UNLIKELY (u == 0)) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK; + + prec = r->_mp_prec; + zeros = 2 * prec - 2; + rsize = zeros + 1 + U2; + + tp = TMP_ALLOC_LIMBS (rsize); + + MPN_ZERO (tp, zeros); + tp[zeros] = u & GMP_NUMB_MASK; + +#if U2 + { + mp_limb_t uhigh = u >> GMP_NUMB_BITS; + tp[zeros + 1] = uhigh; + rsize -= (uhigh == 0); + } +#endif + + mpn_sqrtrem (r->_mp_d, NULL, tp, rsize); + + r->_mp_size = prec; + r->_mp_exp = 1; + TMP_FREE; +} diff --git a/gmp/mpf/sub.c b/gmp/mpf/sub.c new file mode 100644 index 0000000000..3aaf192790 --- /dev/null +++ b/gmp/mpf/sub.c @@ -0,0 +1,419 @@ +/* mpf_sub -- Subtract two floats. + +Copyright 1993-1996, 1999-2002, 2004, 2005, 2011 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +{ + mp_srcptr up, vp; + mp_ptr rp, tp; + mp_size_t usize, vsize, rsize; + mp_size_t prec; + mp_exp_t exp; + mp_size_t ediff; + int negate; + TMP_DECL; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (usize == 0) + { + mpf_neg (r, v); + return; + } + if (vsize == 0) + { + if (r != u) + mpf_set (r, u); + return; + } + + /* If signs of U and V are different, perform addition. */ + if ((usize ^ vsize) < 0) + { + __mpf_struct v_negated; + v_negated._mp_size = -vsize; + v_negated._mp_exp = v->_mp_exp; + v_negated._mp_d = v->_mp_d; + mpf_add (r, u, &v_negated); + return; + } + + TMP_MARK; + + /* Signs are now known to be the same. */ + negate = usize < 0; + + /* Make U be the operand with the largest exponent. */ + if (u->_mp_exp < v->_mp_exp) + { + mpf_srcptr t; + t = u; u = v; v = t; + negate ^= 1; + usize = u->_mp_size; + vsize = v->_mp_size; + } + + usize = ABS (usize); + vsize = ABS (vsize); + up = u->_mp_d; + vp = v->_mp_d; + rp = r->_mp_d; + prec = r->_mp_prec + 1; + exp = u->_mp_exp; + ediff = u->_mp_exp - v->_mp_exp; + + /* If ediff is 0 or 1, we might have a situation where the operands are + extremely close. We need to scan the operands from the most significant + end ignore the initial parts that are equal. */ + if (ediff <= 1) + { + if (ediff == 0) + { + /* Skip leading limbs in U and V that are equal. */ + if (up[usize - 1] == vp[vsize - 1]) + { + /* This loop normally exits immediately. Optimize for that. */ + do + { + usize--; + vsize--; + exp--; + + if (usize == 0) + { + /* u cancels high limbs of v, result is rest of v */ + negate ^= 1; + cancellation: + /* strip high zeros before truncating to prec */ + while (vsize != 0 && vp[vsize - 1] == 0) + { + vsize--; + exp--; + } + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + MPN_COPY_INCR (rp, vp, vsize); + rsize = vsize; + goto done; + } + if (vsize == 0) + { + vp = up; + vsize = usize; + goto cancellation; + } + } + while (up[usize - 1] == vp[vsize - 1]); + } + + if (up[usize - 1] < vp[vsize - 1]) + { + /* For simplicity, swap U and V. Note that since the loop above + wouldn't have exited unless up[usize - 1] and vp[vsize - 1] + were non-equal, this if-statement catches all cases where U + is smaller than V. */ + MPN_SRCPTR_SWAP (up,usize, vp,vsize); + negate ^= 1; + /* negating ediff not necessary since it is 0. */ + } + + /* Check for + x+1 00000000 ... + x ffffffff ... */ + if (up[usize - 1] != vp[vsize - 1] + 1) + goto general_case; + usize--; + vsize--; + exp--; + } + else /* ediff == 1 */ + { + /* Check for + 1 00000000 ... + 0 ffffffff ... */ + + if (up[usize - 1] != 1 || vp[vsize - 1] != GMP_NUMB_MAX + || (usize >= 2 && up[usize - 2] != 0)) + goto general_case; + + usize--; + exp--; + } + + /* Skip sequences of 00000000/ffffffff */ + while (vsize != 0 && usize != 0 && up[usize - 1] == 0 + && vp[vsize - 1] == GMP_NUMB_MAX) + { + usize--; + vsize--; + exp--; + } + + if (usize == 0) + { + while (vsize != 0 && vp[vsize - 1] == GMP_NUMB_MAX) + { + vsize--; + exp--; + } + } + + if (usize > prec - 1) + { + up += usize - (prec - 1); + usize = prec - 1; + } + if (vsize > prec - 1) + { + vp += vsize - (prec - 1); + vsize = prec - 1; + } + + tp = TMP_ALLOC_LIMBS (prec); + { + mp_limb_t cy_limb; + if (vsize == 0) + { + mp_size_t size, i; + size = usize; + for (i = 0; i < size; i++) + tp[i] = up[i]; + tp[size] = 1; + rsize = size + 1; + exp++; + goto normalize; + } + if (usize == 0) + { + mp_size_t size, i; + size = vsize; + for (i = 0; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1); + rsize = vsize; + if (cy_limb == 0) + { + tp[rsize] = 1; + rsize++; + exp++; + } + goto normalize; + } + if (usize >= vsize) + { + /* uuuu */ + /* vv */ + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else /* (usize < vsize) */ + { + /* uuuu */ + /* vvvvvvv */ + mp_size_t size, i; + size = vsize - usize; + for (i = 0; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + cy_limb = mpn_sub_n (tp + size, up, vp + size, usize); + cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1); + rsize = vsize; + } + if (cy_limb == 0) + { + tp[rsize] = 1; + rsize++; + exp++; + } + goto normalize; + } + } + +general_case: + /* If U extends beyond PREC, ignore the part that does. */ + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + + /* If V extends beyond PREC, ignore the part that does. + Note that this may make vsize negative. */ + if (vsize + ediff > prec) + { + vp += vsize + ediff - prec; + vsize = prec - ediff; + } + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (rp != up) + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + /* Allocate temp space for the result. Allocate + just vsize + ediff later??? */ + tp = TMP_ALLOC_LIMBS (prec); + + /* Locate the least significant non-zero limb in (the needed + parts of) U and V, to simplify the code below. */ + for (;;) + { + if (vsize == 0) + { + MPN_COPY (rp, up, usize); + rsize = usize; + goto done; + } + if (vp[0] != 0) + break; + vp++, vsize--; + } + for (;;) + { + if (usize == 0) + { + MPN_COPY (rp, vp, vsize); + rsize = vsize; + negate ^= 1; + goto done; + } + if (up[0] != 0) + break; + up++, usize--; + } + + /* uuuu | uuuu | uuuu | uuuu | uuuu */ + /* vvvvvvv | vv | vvvvv | v | vv */ + + if (usize > ediff) + { + /* U and V partially overlaps. */ + if (ediff == 0) + { + /* Have to compare the leading limbs of u and v + to determine whether to compute u - v or v - u. */ + if (usize >= vsize) + { + /* uuuu */ + /* vv */ + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else /* (usize < vsize) */ + { + /* uuuu */ + /* vvvvvvv */ + mp_size_t size, i; + size = vsize - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + mpn_sub_n (tp + size, up, vp + size, usize); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize; + } + } + else + { + if (vsize + ediff <= usize) + { + /* uuuu */ + /* v */ + mp_size_t size; + size = usize - ediff - vsize; + MPN_COPY (tp, up, size); + mpn_sub (tp + size, up + size, usize - size, vp, vsize); + rsize = usize; + } + else + { + /* uuuu */ + /* vvvvv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + mpn_sub (tp + size, up, usize, vp + size, usize - ediff); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize + ediff; + } + } + } + else + { + /* uuuu */ + /* vv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < vsize; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + for (i = vsize; i < size; i++) + tp[i] = GMP_NUMB_MAX; + mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1); + rsize = size + usize; + } + + normalize: + /* Full normalize. Optimize later. */ + while (rsize != 0 && tp[rsize - 1] == 0) + { + rsize--; + exp--; + } + MPN_COPY (rp, tp, rsize); + } + + done: + r->_mp_size = negate ? -rsize : rsize; + if (rsize == 0) + exp = 0; + r->_mp_exp = exp; + TMP_FREE; +} diff --git a/gmp/mpf/sub_ui.c b/gmp/mpf/sub_ui.c new file mode 100644 index 0000000000..cf9b88eb00 --- /dev/null +++ b/gmp/mpf/sub_ui.c @@ -0,0 +1,51 @@ +/* mpf_sub_ui -- Subtract an unsigned integer from a float. + +Copyright 1993, 1994, 1996, 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_sub_ui (mpf_ptr sum, mpf_srcptr u, unsigned long int v) +{ + __mpf_struct vv; + mp_limb_t vl; + + if (v == 0) + { + mpf_set (sum, u); + return; + } + + vl = v; + vv._mp_size = 1; + vv._mp_d = &vl; + vv._mp_exp = 1; + mpf_sub (sum, u, &vv); +} diff --git a/gmp/mpf/swap.c b/gmp/mpf/swap.c new file mode 100644 index 0000000000..a370652876 --- /dev/null +++ b/gmp/mpf/swap.c @@ -0,0 +1,57 @@ +/* mpf_swap (U, V) -- Swap U and V. + +Copyright 1997, 1998, 2000, 2001, 2013 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_swap (mpf_ptr u, mpf_ptr v) __GMP_NOTHROW +{ + mp_ptr tptr; + mp_size_t tprec; + mp_size_t tsiz; + mp_exp_t texp; + + tprec = PREC(u); + PREC(u) = PREC(v); + PREC(v) = tprec; + + tsiz = SIZ(u); + SIZ(u) = SIZ(v); + SIZ(v) = tsiz; + + texp = EXP(u); + EXP(u) = EXP(v); + EXP(v) = texp; + + tptr = PTR(u); + PTR(u) = PTR(v); + PTR(v) = tptr; +} diff --git a/gmp/mpf/trunc.c b/gmp/mpf/trunc.c new file mode 100644 index 0000000000..5f94f7aec8 --- /dev/null +++ b/gmp/mpf/trunc.c @@ -0,0 +1,75 @@ +/* mpf_trunc -- truncate an mpf to an integer. + +Copyright 2001 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + + +/* Notice the use of prec+1 ensures mpf_trunc is equivalent to mpf_set if u + is already an integer. */ + +void +mpf_trunc (mpf_ptr r, mpf_srcptr u) +{ + mp_ptr rp; + mp_srcptr up; + mp_size_t size, asize, prec; + mp_exp_t exp; + + exp = EXP(u); + size = SIZ(u); + if (size == 0 || exp <= 0) + { + /* u is only a fraction */ + SIZ(r) = 0; + EXP(r) = 0; + return; + } + + up = PTR(u); + EXP(r) = exp; + asize = ABS (size); + up += asize; + + /* skip fraction part of u */ + asize = MIN (asize, exp); + + /* don't lose precision in the copy */ + prec = PREC(r) + 1; + + /* skip excess over target precision */ + asize = MIN (asize, prec); + + up -= asize; + rp = PTR(r); + SIZ(r) = (size >= 0 ? asize : -asize); + if (rp != up) + MPN_COPY_INCR (rp, up, asize); +} diff --git a/gmp/mpf/ui_div.c b/gmp/mpf/ui_div.c new file mode 100644 index 0000000000..ceb881ebcc --- /dev/null +++ b/gmp/mpf/ui_div.c @@ -0,0 +1,128 @@ +/* mpf_ui_div -- Divide an unsigned integer with a float. + +Copyright 1993-1996, 2000-2002, 2004, 2005, 2012 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include <stdio.h> /* for NULL */ +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + + +void +mpf_ui_div (mpf_ptr r, unsigned long int u, mpf_srcptr v) +{ + mp_srcptr vp; + mp_ptr rp, tp, remp, new_vp; + mp_size_t vsize; + mp_size_t rsize, prospective_rsize, zeros, tsize, high_zero; + mp_size_t sign_quotient; + mp_size_t prec; + mp_exp_t rexp; + TMP_DECL; + + vsize = v->_mp_size; + sign_quotient = vsize; + + if (UNLIKELY (vsize == 0)) + DIVIDE_BY_ZERO; + + if (UNLIKELY (u == 0)) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + vsize = ABS (vsize); + prec = r->_mp_prec; + + TMP_MARK; + rexp = 1 - v->_mp_exp + 1; + + rp = r->_mp_d; + vp = v->_mp_d; + + prospective_rsize = 1 - vsize + 1; /* quot from using given u,v sizes */ + rsize = prec + 1; /* desired quot size */ + + zeros = rsize - prospective_rsize; /* padding u to give rsize */ + tsize = 1 + zeros; /* u with zeros */ + + if (WANT_TMP_DEBUG) + { + /* separate alloc blocks, for malloc debugging */ + remp = TMP_ALLOC_LIMBS (vsize); + tp = TMP_ALLOC_LIMBS (tsize); + new_vp = NULL; + if (rp == vp) + new_vp = TMP_ALLOC_LIMBS (vsize); + } + else + { + /* one alloc with calculated size, for efficiency */ + mp_size_t size = vsize + tsize + (rp == vp ? vsize : 0); + remp = TMP_ALLOC_LIMBS (size); + tp = remp + vsize; + new_vp = tp + tsize; + } + + /* ensure divisor doesn't overlap quotient */ + if (rp == vp) + { + MPN_COPY (new_vp, vp, vsize); + vp = new_vp; + } + + MPN_ZERO (tp, tsize-1); + + tp[tsize - 1] = u & GMP_NUMB_MASK; +#if BITS_PER_ULONG > GMP_NUMB_BITS + if (u > GMP_NUMB_MAX) + { + /* tsize-vsize+1 == rsize, so tsize >= rsize. rsize == prec+1 >= 2, + so tsize >= 2, hence there's room for 2-limb u with nails */ + ASSERT (tsize >= 2); + tp[tsize - 1] = u >> GMP_NUMB_BITS; + tp[tsize - 2] = u & GMP_NUMB_MASK; + rexp++; + } +#endif + + ASSERT (tsize-vsize+1 == rsize); + mpn_tdiv_qr (rp, remp, (mp_size_t) 0, tp, tsize, vp, vsize); + + /* strip possible zero high limb */ + high_zero = (rp[rsize-1] == 0); + rsize -= high_zero; + rexp -= high_zero; + + r->_mp_size = sign_quotient >= 0 ? rsize : -rsize; + r->_mp_exp = rexp; + TMP_FREE; +} diff --git a/gmp/mpf/ui_sub.c b/gmp/mpf/ui_sub.c new file mode 100644 index 0000000000..b7a536eb8b --- /dev/null +++ b/gmp/mpf/ui_sub.c @@ -0,0 +1,336 @@ +/* mpf_ui_sub -- Subtract a float from an unsigned long int. + +Copyright 1993-1996, 2001, 2002, 2005 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v) +{ + mp_srcptr up, vp; + mp_ptr rp, tp; + mp_size_t usize, vsize, rsize; + mp_size_t prec; + mp_exp_t uexp; + mp_size_t ediff; + int negate; + mp_limb_t ulimb; + TMP_DECL; + + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (u == 0) + { + mpf_neg (r, v); + return; + } + if (vsize == 0) + { + mpf_set_ui (r, u); + return; + } + + /* If signs of U and V are different, perform addition. */ + if (vsize < 0) + { + __mpf_struct v_negated; + v_negated._mp_size = -vsize; + v_negated._mp_exp = v->_mp_exp; + v_negated._mp_d = v->_mp_d; + mpf_add_ui (r, &v_negated, u); + return; + } + + TMP_MARK; + + /* Signs are now known to be the same. */ + + ulimb = u; + /* Make U be the operand with the largest exponent. */ + if (1 < v->_mp_exp) + { + negate = 1; + usize = ABS (vsize); + vsize = 1; + up = v->_mp_d; + vp = &ulimb; + rp = r->_mp_d; + prec = r->_mp_prec + 1; + uexp = v->_mp_exp; + ediff = uexp - 1; + } + else + { + negate = 0; + usize = 1; + vsize = ABS (vsize); + up = &ulimb; + vp = v->_mp_d; + rp = r->_mp_d; + prec = r->_mp_prec; + uexp = 1; + ediff = 1 - v->_mp_exp; + } + + /* Ignore leading limbs in U and V that are equal. Doing + this helps increase the precision of the result. */ + if (ediff == 0) + { + /* This loop normally exits immediately. Optimize for that. */ + for (;;) + { + usize--; + vsize--; + if (up[usize] != vp[vsize]) + break; + uexp--; + if (usize == 0) + goto Lu0; + if (vsize == 0) + goto Lv0; + } + usize++; + vsize++; + /* Note that either operand (but not both operands) might now have + leading zero limbs. It matters only that U is unnormalized if + vsize is now zero, and vice versa. And it is only in that case + that we have to adjust uexp. */ + if (vsize == 0) + Lv0: + while (usize != 0 && up[usize - 1] == 0) + usize--, uexp--; + if (usize == 0) + Lu0: + while (vsize != 0 && vp[vsize - 1] == 0) + vsize--, uexp--; + } + + /* If U extends beyond PREC, ignore the part that does. */ + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + + /* If V extends beyond PREC, ignore the part that does. + Note that this may make vsize negative. */ + if (vsize + ediff > prec) + { + vp += vsize + ediff - prec; + vsize = prec - ediff; + } + + /* Allocate temp space for the result. Allocate + just vsize + ediff later??? */ + tp = TMP_ALLOC_LIMBS (prec); + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (tp != up) + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + /* Locate the least significant non-zero limb in (the needed + parts of) U and V, to simplify the code below. */ + for (;;) + { + if (vsize == 0) + { + MPN_COPY (rp, up, usize); + rsize = usize; + goto done; + } + if (vp[0] != 0) + break; + vp++, vsize--; + } + for (;;) + { + if (usize == 0) + { + MPN_COPY (rp, vp, vsize); + rsize = vsize; + negate ^= 1; + goto done; + } + if (up[0] != 0) + break; + up++, usize--; + } + + /* uuuu | uuuu | uuuu | uuuu | uuuu */ + /* vvvvvvv | vv | vvvvv | v | vv */ + + if (usize > ediff) + { + /* U and V partially overlaps. */ + if (ediff == 0) + { + /* Have to compare the leading limbs of u and v + to determine whether to compute u - v or v - u. */ + if (usize > vsize) + { + /* uuuu */ + /* vv */ + int cmp; + cmp = mpn_cmp (up + usize - vsize, vp, vsize); + if (cmp >= 0) + { + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else + { + /* vv */ /* Swap U and V. */ + /* uuuu */ + mp_size_t size, i; + size = usize - vsize; + tp[0] = -up[0] & GMP_NUMB_MASK; + for (i = 1; i < size; i++) + tp[i] = ~up[i] & GMP_NUMB_MASK; + mpn_sub_n (tp + size, vp, up + size, vsize); + mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1); + negate ^= 1; + rsize = usize; + } + } + else if (usize < vsize) + { + /* uuuu */ + /* vvvvvvv */ + int cmp; + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp > 0) + { + mp_size_t size, i; + size = vsize - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + mpn_sub_n (tp + size, up, vp + size, usize); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize; + } + else + { + /* vvvvvvv */ /* Swap U and V. */ + /* uuuu */ + /* This is the only place we can get 0.0. */ + mp_size_t size; + size = vsize - usize; + MPN_COPY (tp, vp, size); + mpn_sub_n (tp + size, vp + size, up, usize); + negate ^= 1; + rsize = vsize; + } + } + else + { + /* uuuu */ + /* vvvv */ + int cmp; + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp > 0) + { + mpn_sub_n (tp, up, vp, usize); + rsize = usize; + } + else + { + mpn_sub_n (tp, vp, up, usize); + negate ^= 1; + rsize = usize; + /* can give zero */ + } + } + } + else + { + if (vsize + ediff <= usize) + { + /* uuuu */ + /* v */ + mp_size_t size; + size = usize - ediff - vsize; + MPN_COPY (tp, up, size); + mpn_sub (tp + size, up + size, usize - size, vp, vsize); + rsize = usize; + } + else + { + /* uuuu */ + /* vvvvv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < size; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + mpn_sub (tp + size, up, usize, vp + size, usize - ediff); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize + ediff; + } + } + } + else + { + /* uuuu */ + /* vv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0] & GMP_NUMB_MASK; + for (i = 1; i < vsize; i++) + tp[i] = ~vp[i] & GMP_NUMB_MASK; + for (i = vsize; i < size; i++) + tp[i] = GMP_NUMB_MAX; + mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1); + rsize = size + usize; + } + + /* Full normalize. Optimize later. */ + while (rsize != 0 && tp[rsize - 1] == 0) + { + rsize--; + uexp--; + } + MPN_COPY (rp, tp, rsize); + } + + done: + r->_mp_size = negate ? -rsize : rsize; + r->_mp_exp = uexp; + TMP_FREE; +} diff --git a/gmp/mpf/urandomb.c b/gmp/mpf/urandomb.c new file mode 100644 index 0000000000..72271e8762 --- /dev/null +++ b/gmp/mpf/urandomb.c @@ -0,0 +1,69 @@ +/* mpf_urandomb (rop, state, nbits) -- Generate a uniform pseudorandom + real number between 0 (inclusive) and 1 (exclusive) of size NBITS, + using STATE as the random state previously initialized by a call to + gmp_randinit(). + +Copyright 1999-2002 Free Software Foundation, Inc. + +This file is part of the GNU MP Library. + +The GNU MP Library is free software; you can redistribute it and/or modify +it under the terms of either: + + * the GNU Lesser General Public License as published by the Free + Software Foundation; either version 3 of the License, or (at your + option) any later version. + +or + + * the GNU General Public License as published by the Free Software + Foundation; either version 2 of the License, or (at your option) any + later version. + +or both in parallel, as here. + +The GNU MP 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 General Public License +for more details. + +You should have received copies of the GNU General Public License and the +GNU Lesser General Public License along with the GNU MP Library. If not, +see https://www.gnu.org/licenses/. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +mpf_urandomb (mpf_t rop, gmp_randstate_t rstate, mp_bitcnt_t nbits) +{ + mp_ptr rp; + mp_size_t nlimbs; + mp_exp_t exp; + mp_size_t prec; + + rp = PTR (rop); + nlimbs = BITS_TO_LIMBS (nbits); + prec = PREC (rop); + + if (nlimbs > prec + 1 || nlimbs == 0) + { + nlimbs = prec + 1; + nbits = nlimbs * GMP_NUMB_BITS; + } + + _gmp_rand (rp, rstate, nbits); + + /* If nbits isn't a multiple of GMP_NUMB_BITS, shift up. */ + if (nbits % GMP_NUMB_BITS != 0) + mpn_lshift (rp, rp, nlimbs, GMP_NUMB_BITS - nbits % GMP_NUMB_BITS); + + exp = 0; + while (nlimbs != 0 && rp[nlimbs - 1] == 0) + { + nlimbs--; + exp--; + } + EXP (rop) = exp; + SIZ (rop) = nlimbs; +} |