Revert "Remove gmp from AOS"
This reverts commit f37c97684f0910a3f241394549392f00145ab0f7.
We need gmp for SymEngine for symbolicmanipultion in C++
Change-Id: Ia13216d1715cf96944f7b4f422b7a799f921d4a4
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/third_party/gmp/mpf/Makefile.am b/third_party/gmp/mpf/Makefile.am
new file mode 100644
index 0000000..ff7e48f
--- /dev/null
+++ b/third_party/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/.
+
+
+AM_CPPFLAGS = -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_z.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/third_party/gmp/mpf/Makefile.in b/third_party/gmp/mpf/Makefile.in
new file mode 100644
index 0000000..8acf389
--- /dev/null
+++ b/third_party/gmp/mpf/Makefile.in
@@ -0,0 +1,662 @@
+# Makefile.in generated by automake 1.15 from Makefile.am.
+# @configure_input@
+
+# Copyright (C) 1994-2014 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__is_gnu_make = { \
+ if test -z '$(MAKELEVEL)'; then \
+ false; \
+ elif test -n '$(MAKE_HOST)'; then \
+ true; \
+ elif test -n '$(MAKE_VERSION)' && test -n '$(CURDIR)'; then \
+ true; \
+ else \
+ false; \
+ fi; \
+}
+am__make_running_with_option = \
+ case $${target_option-} in \
+ ?) ;; \
+ *) echo "am__make_running_with_option: internal error: invalid" \
+ "target option '$${target_option-}' specified" >&2; \
+ exit 1;; \
+ esac; \
+ has_opt=no; \
+ sane_makeflags=$$MAKEFLAGS; \
+ if $(am__is_gnu_make); then \
+ sane_makeflags=$$MFLAGS; \
+ else \
+ case $$MAKEFLAGS in \
+ *\\[\ \ ]*) \
+ bs=\\; \
+ sane_makeflags=`printf '%s\n' "$$MAKEFLAGS" \
+ | sed "s/$$bs$$bs[$$bs $$bs ]*//g"`;; \
+ esac; \
+ fi; \
+ skip_next=no; \
+ strip_trailopt () \
+ { \
+ flg=`printf '%s\n' "$$flg" | sed "s/$$1.*$$//"`; \
+ }; \
+ for flg in $$sane_makeflags; do \
+ test $$skip_next = yes && { skip_next=no; continue; }; \
+ case $$flg in \
+ *=*|--*) continue;; \
+ -*I) strip_trailopt 'I'; skip_next=yes;; \
+ -*I?*) strip_trailopt 'I';; \
+ -*O) strip_trailopt 'O'; skip_next=yes;; \
+ -*O?*) strip_trailopt 'O';; \
+ -*l) strip_trailopt 'l'; skip_next=yes;; \
+ -*l?*) strip_trailopt 'l';; \
+ -[dEDm]) skip_next=yes;; \
+ -[JT]) skip_next=yes;; \
+ esac; \
+ case $$flg in \
+ *$$target_option*) has_opt=yes; break;; \
+ esac; \
+ done; \
+ test $$has_opt = yes
+am__make_dryrun = (target_option=n; $(am__make_running_with_option))
+am__make_keepgoing = (target_option=k; $(am__make_running_with_option))
+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
+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)
+DIST_COMMON = $(srcdir)/Makefile.am $(am__DIST_COMMON)
+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_z.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)
+AM_V_lt = $(am__v_lt_@AM_V@)
+am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@)
+am__v_lt_0 = --silent
+am__v_lt_1 =
+AM_V_P = $(am__v_P_@AM_V@)
+am__v_P_ = $(am__v_P_@AM_DEFAULT_V@)
+am__v_P_0 = false
+am__v_P_1 = :
+AM_V_GEN = $(am__v_GEN_@AM_V@)
+am__v_GEN_ = $(am__v_GEN_@AM_DEFAULT_V@)
+am__v_GEN_0 = @echo " GEN " $@;
+am__v_GEN_1 =
+AM_V_at = $(am__v_at_@AM_V@)
+am__v_at_ = $(am__v_at_@AM_DEFAULT_V@)
+am__v_at_0 = @
+am__v_at_1 =
+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) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \
+ $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) \
+ $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) \
+ $(AM_CFLAGS) $(CFLAGS)
+AM_V_CC = $(am__v_CC_@AM_V@)
+am__v_CC_ = $(am__v_CC_@AM_DEFAULT_V@)
+am__v_CC_0 = @echo " CC " $@;
+am__v_CC_1 =
+CCLD = $(CC)
+LINK = $(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) \
+ $(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
+ $(AM_LDFLAGS) $(LDFLAGS) -o $@
+AM_V_CCLD = $(am__v_CCLD_@AM_V@)
+am__v_CCLD_ = $(am__v_CCLD_@AM_DEFAULT_V@)
+am__v_CCLD_0 = @echo " CCLD " $@;
+am__v_CCLD_1 =
+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
+am__tagged_files = $(HEADERS) $(SOURCES) $(TAGS_FILES) $(LISP)
+# Read a list of newline-separated strings from the standard input,
+# and print each of them once, without duplicates. Input order is
+# *not* preserved.
+am__uniquify_input = $(AWK) '\
+ BEGIN { nonempty = 0; } \
+ { items[$$0] = 1; nonempty = 1; } \
+ END { if (nonempty) { for (i in items) print i; }; } \
+'
+# Make sure the list of sources is unique. This is necessary because,
+# e.g., the same source file might be shared among _SOURCES variables
+# for different programs/libraries.
+am__define_uniq_tagged_files = \
+ list='$(am__tagged_files)'; \
+ unique=`for i in $$list; do \
+ if test -f "$$i"; then echo $$i; else echo $(srcdir)/$$i; fi; \
+ done | $(am__uniquify_input)`
+ETAGS = etags
+CTAGS = ctags
+am__DIST_COMMON = $(srcdir)/Makefile.in
+DISTFILES = $(DIST_COMMON) $(DIST_SOURCES) $(TEXINFOS) $(EXTRA_DIST)
+ABI = @ABI@
+ACLOCAL = @ACLOCAL@
+AMTAR = @AMTAR@
+AM_DEFAULT_VERBOSITY = @AM_DEFAULT_VERBOSITY@
+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@
+LT_SYS_LIBRARY_PATH = @LT_SYS_LIBRARY_PATH@
+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@
+AM_CPPFLAGS = -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_z.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
+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)'; \
+ locs=`for p in $$list; do echo $$p; done | \
+ sed 's|^[^/]*$$|.|; s|/[^/]*$$||; s|$$|/so_locations|' | \
+ sort -u`; \
+ test -z "$$locs" || { \
+ echo rm -f $${locs}; \
+ rm -f $${locs}; \
+ }
+
+libmpf.la: $(libmpf_la_OBJECTS) $(libmpf_la_DEPENDENCIES) $(EXTRA_libmpf_la_DEPENDENCIES)
+ $(AM_V_CCLD)$(LINK) $(libmpf_la_OBJECTS) $(libmpf_la_LIBADD) $(LIBS)
+
+mostlyclean-compile:
+ -rm -f *.$(OBJEXT)
+
+distclean-compile:
+ -rm -f *.tab.c
+
+.c.o:
+ $(AM_V_CC)$(COMPILE) -c -o $@ $<
+
+.c.obj:
+ $(AM_V_CC)$(COMPILE) -c -o $@ `$(CYGPATH_W) '$<'`
+
+.c.lo:
+ $(AM_V_CC)$(LTCOMPILE) -c -o $@ $<
+
+mostlyclean-libtool:
+ -rm -f *.lo
+
+clean-libtool:
+ -rm -rf .libs _libs
+
+ID: $(am__tagged_files)
+ $(am__define_uniq_tagged_files); mkid -fID $$unique
+tags: tags-am
+TAGS: tags
+
+tags-am: $(TAGS_DEPENDENCIES) $(am__tagged_files)
+ set x; \
+ here=`pwd`; \
+ $(am__define_uniq_tagged_files); \
+ 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-am
+
+CTAGS: ctags
+ctags-am: $(TAGS_DEPENDENCIES) $(am__tagged_files)
+ $(am__define_uniq_tagged_files); \
+ 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"
+cscopelist: cscopelist-am
+
+cscopelist-am: $(am__tagged_files)
+ list='$(am__tagged_files)'; \
+ case "$(srcdir)" in \
+ [\\/]* | ?:[\\/]*) sdir="$(srcdir)" ;; \
+ *) sdir=$(subdir)/$(srcdir) ;; \
+ esac; \
+ for i in $$list; do \
+ if test -f "$$i"; then \
+ echo "$(subdir)/$$i"; \
+ else \
+ echo "$$sdir/$$i"; \
+ fi; \
+ done >> $(top_builddir)/cscope.files
+
+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 TAGS all all-am check check-am clean clean-generic \
+ clean-libtool clean-noinstLTLIBRARIES cscopelist-am ctags \
+ ctags-am 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 tags-am uninstall uninstall-am
+
+.PRECIOUS: Makefile
+
+
+# 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/third_party/gmp/mpf/abs.c b/third_party/gmp/mpf/abs.c
new file mode 100644
index 0000000..1642a46
--- /dev/null
+++ b/third_party/gmp/mpf/abs.c
@@ -0,0 +1,58 @@
+/* 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-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/third_party/gmp/mpf/add.c b/third_party/gmp/mpf/add.c
new file mode 100644
index 0000000..77a9e47
--- /dev/null
+++ b/third_party/gmp/mpf/add.c
@@ -0,0 +1,183 @@
+/* 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-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/third_party/gmp/mpf/add_ui.c b/third_party/gmp/mpf/add_ui.c
new file mode 100644
index 0000000..1e2a94b
--- /dev/null
+++ b/third_party/gmp/mpf/add_ui.c
@@ -0,0 +1,152 @@
+/* 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-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/third_party/gmp/mpf/ceilfloor.c b/third_party/gmp/mpf/ceilfloor.c
new file mode 100644
index 0000000..9bb6638
--- /dev/null
+++ b/third_party/gmp/mpf/ceilfloor.c
@@ -0,0 +1,125 @@
+/* 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-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/third_party/gmp/mpf/clear.c b/third_party/gmp/mpf/clear.c
new file mode 100644
index 0000000..0939e03
--- /dev/null
+++ b/third_party/gmp/mpf/clear.c
@@ -0,0 +1,38 @@
+/* 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-impl.h"
+
+void
+mpf_clear (mpf_ptr x)
+{
+ __GMP_FREE_FUNC_LIMBS (PTR(x), PREC(x) + 1);
+}
diff --git a/third_party/gmp/mpf/clears.c b/third_party/gmp/mpf/clears.c
new file mode 100644
index 0000000..115fa19
--- /dev/null
+++ b/third_party/gmp/mpf/clears.c
@@ -0,0 +1,49 @@
+/* mpf_clears() -- Clear multiple mpf_t variables.
+
+Copyright 2009, 2014, 2015 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 "gmp-impl.h"
+
+void
+mpf_clears (mpf_ptr x, ...)
+{
+ va_list ap;
+
+ va_start (ap, x);
+
+ do
+ {
+ __GMP_FREE_FUNC_LIMBS (PTR(x), PREC(x) + 1);
+ x = va_arg (ap, mpf_ptr);
+ }
+ while (x != NULL);
+
+ va_end (ap);
+}
diff --git a/third_party/gmp/mpf/cmp.c b/third_party/gmp/mpf/cmp.c
new file mode 100644
index 0000000..3518b51
--- /dev/null
+++ b/third_party/gmp/mpf/cmp.c
@@ -0,0 +1,113 @@
+/* mpf_cmp -- Compare two floats.
+
+Copyright 1993, 1994, 1996, 2001, 2015 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-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;
+
+ usize = SIZ(u);
+ vsize = SIZ(v);
+ usign = usize >= 0 ? 1 : -1;
+
+ /* 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 usign;
+ }
+
+ /* U and V have the same sign and are both non-zero. */
+
+ uexp = EXP(u);
+ vexp = EXP(v);
+
+ /* 2. Are the exponents different? */
+ if (uexp > vexp)
+ return usign;
+ if (uexp < vexp)
+ return -usign;
+
+ usize = ABS (usize);
+ vsize = ABS (vsize);
+
+ up = PTR (u);
+ vp = PTR (v);
+
+#define STRICT_MPF_NORMALIZATION 0
+#if ! STRICT_MPF_NORMALIZATION
+ /* Ignore zeroes at the low end of U and V. */
+ do {
+ mp_limb_t tl;
+ tl = up[0];
+ MPN_STRIP_LOW_ZEROS_NOT_ZERO (up, usize, tl);
+ tl = vp[0];
+ MPN_STRIP_LOW_ZEROS_NOT_ZERO (vp, vsize, tl);
+ } while (0);
+#endif
+
+ if (usize > vsize)
+ {
+ cmp = mpn_cmp (up + usize - vsize, vp, vsize);
+ /* if (cmp == 0) */
+ /* return usign; */
+ ++cmp;
+ }
+ 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/third_party/gmp/mpf/cmp_d.c b/third_party/gmp/mpf/cmp_d.c
new file mode 100644
index 0000000..3fa099b
--- /dev/null
+++ b/third_party/gmp/mpf/cmp_d.c
@@ -0,0 +1,59 @@
+/* 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-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/third_party/gmp/mpf/cmp_si.c b/third_party/gmp/mpf/cmp_si.c
new file mode 100644
index 0000000..d8d9880
--- /dev/null
+++ b/third_party/gmp/mpf/cmp_si.c
@@ -0,0 +1,109 @@
+/* mpf_cmp_si -- Compare a float with a signed integer.
+
+Copyright 1993-1995, 1999-2002, 2004, 2012, 2015 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-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;
+
+ usize = SIZ (u);
+
+ /* 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. */
+
+ /* 2. Are the exponents different (V's exponent == 1)? */
+ uexp = EXP (u);
+ usign = usize >= 0 ? 1 : -1;
+ usize = ABS (usize);
+ abs_vval = ABS_CAST (unsigned long, vval);
+
+#if GMP_NAIL_BITS != 0
+ if (uexp != 1 + (abs_vval > GMP_NUMB_MAX))
+ return (uexp < 1 + (abs_vval > GMP_NUMB_MAX)) ? -usign : usign;
+#else
+ if (uexp != 1)
+ return (uexp < 1) ? -usign : usign;
+#endif
+
+ up = PTR (u);
+
+ ASSERT (usize > 0);
+ ulimb = up[--usize];
+#if GMP_NAIL_BITS != 0
+ if (uexp == 2)
+ {
+ if ((ulimb >> GMP_NAIL_BITS) != 0)
+ return usign;
+ ulimb = (ulimb << GMP_NUMB_BITS);
+ if (usize != 0) ulimb |= up[--usize];
+ }
+#endif
+
+ /* 3. Compare the most significant mantissa limb with V. */
+ if (ulimb != abs_vval)
+ return (ulimb < abs_vval) ? -usign : usign;
+
+ /* Ignore zeroes at the low end of U. */
+ for (; *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/third_party/gmp/mpf/cmp_ui.c b/third_party/gmp/mpf/cmp_ui.c
new file mode 100644
index 0000000..a9a6036
--- /dev/null
+++ b/third_party/gmp/mpf/cmp_ui.c
@@ -0,0 +1,87 @@
+/* mpf_cmp_ui -- Compare a float with an unsigned integer.
+
+Copyright 1993-1995, 1999, 2001, 2002, 2015 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-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;
+
+ usize = SIZ (u);
+
+ /* 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)? */
+ uexp = EXP (u);
+
+#if GMP_NAIL_BITS != 0
+ if (uexp != 1 + (vval > GMP_NUMB_MAX))
+ return (uexp < 1 + (vval > GMP_NUMB_MAX)) ? -1 : 1;
+#else
+ if (uexp != 1)
+ return (uexp < 1) ? -1 : 1;
+#endif
+
+ up = PTR (u);
+
+ ASSERT (usize > 0);
+ ulimb = up[--usize];
+#if GMP_NAIL_BITS != 0
+ if (uexp == 2)
+ {
+ if ((ulimb >> GMP_NAIL_BITS) != 0)
+ return 1;
+ ulimb = (ulimb << GMP_NUMB_BITS);
+ if (usize != 0) ulimb |= up[--usize];
+ }
+#endif
+
+ /* 3. Compare the most significant mantissa limb with V. */
+ if (ulimb != vval)
+ return (ulimb < vval) ? -1 : 1;
+
+ /* Ignore zeroes at the low end of U. */
+ for (; *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. */
+ return (usize > 0);
+}
diff --git a/third_party/gmp/mpf/cmp_z.c b/third_party/gmp/mpf/cmp_z.c
new file mode 100644
index 0000000..279980f
--- /dev/null
+++ b/third_party/gmp/mpf/cmp_z.c
@@ -0,0 +1,45 @@
+/* mpf_cmp_z -- Compare a float with an integer.
+
+Copyright 2015 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-impl.h"
+
+int
+mpf_cmp_z (mpf_srcptr u, mpz_srcptr v) __GMP_NOTHROW
+{
+ mpf_t vf;
+ mp_size_t size;
+
+ SIZ (vf) = size = SIZ (v);
+ EXP (vf) = size = ABS (size);
+ /* PREC (vf) = size; */
+ PTR (vf) = PTR (v);
+
+ return mpf_cmp (u, vf);
+}
diff --git a/third_party/gmp/mpf/div.c b/third_party/gmp/mpf/div.c
new file mode 100644
index 0000000..d13af75
--- /dev/null
+++ b/third_party/gmp/mpf/div.c
@@ -0,0 +1,137 @@
+/* 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-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/third_party/gmp/mpf/div_2exp.c b/third_party/gmp/mpf/div_2exp.c
new file mode 100644
index 0000000..ad552c1
--- /dev/null
+++ b/third_party/gmp/mpf/div_2exp.c
@@ -0,0 +1,138 @@
+/* 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-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/third_party/gmp/mpf/div_ui.c b/third_party/gmp/mpf/div_ui.c
new file mode 100644
index 0000000..e1b0112
--- /dev/null
+++ b/third_party/gmp/mpf/div_ui.c
@@ -0,0 +1,110 @@
+/* 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-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/third_party/gmp/mpf/dump.c b/third_party/gmp/mpf/dump.c
new file mode 100644
index 0000000..cd37dab
--- /dev/null
+++ b/third_party/gmp/mpf/dump.c
@@ -0,0 +1,52 @@
+/* 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-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/third_party/gmp/mpf/eq.c b/third_party/gmp/mpf/eq.c
new file mode 100644
index 0000000..cddb9d5
--- /dev/null
+++ b/third_party/gmp/mpf/eq.c
@@ -0,0 +1,149 @@
+/* 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-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/third_party/gmp/mpf/fits_s.h b/third_party/gmp/mpf/fits_s.h
new file mode 100644
index 0000000..80e74be
--- /dev/null
+++ b/third_party/gmp/mpf/fits_s.h
@@ -0,0 +1,71 @@
+/* mpf_fits_s*_p -- test whether an mpf fits a C signed type.
+
+Copyright 2001, 2002, 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 "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;
+
+ exp = EXP(f);
+ if (exp < 1)
+ return 1; /* -1 < f < 1 truncates to zero, so fits */
+
+ fs = SIZ (f);
+ 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 : NEG_CAST (mp_limb_t, MINIMUM));
+}
diff --git a/third_party/gmp/mpf/fits_sint.c b/third_party/gmp/mpf/fits_sint.c
new file mode 100644
index 0000000..26ace07
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/fits_slong.c b/third_party/gmp/mpf/fits_slong.c
new file mode 100644
index 0000000..25db68c
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/fits_sshort.c b/third_party/gmp/mpf/fits_sshort.c
new file mode 100644
index 0000000..3bfc5a4
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/fits_u.h b/third_party/gmp/mpf/fits_u.h
new file mode 100644
index 0000000..bd7ca78
--- /dev/null
+++ b/third_party/gmp/mpf/fits_u.h
@@ -0,0 +1,73 @@
+/* mpf_fits_u*_p -- test whether an mpf fits a C unsigned type.
+
+Copyright 2001, 2002, 2013, 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 "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) /* zero catched by exp == 0 */
+ return 0; /* negatives don't fit */
+
+ 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/third_party/gmp/mpf/fits_uint.c b/third_party/gmp/mpf/fits_uint.c
new file mode 100644
index 0000000..4b107b0
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/fits_ulong.c b/third_party/gmp/mpf/fits_ulong.c
new file mode 100644
index 0000000..1db688c
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/fits_ushort.c b/third_party/gmp/mpf/fits_ushort.c
new file mode 100644
index 0000000..76a3fd9
--- /dev/null
+++ b/third_party/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/third_party/gmp/mpf/get_d.c b/third_party/gmp/mpf/get_d.c
new file mode 100644
index 0000000..34826fb
--- /dev/null
+++ b/third_party/gmp/mpf/get_d.c
@@ -0,0 +1,46 @@
+/* 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-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/third_party/gmp/mpf/get_d_2exp.c b/third_party/gmp/mpf/get_d_2exp.c
new file mode 100644
index 0000000..440a753
--- /dev/null
+++ b/third_party/gmp/mpf/get_d_2exp.c
@@ -0,0 +1,56 @@
+/* double mpf_get_d_2exp (signed long int *exp, mpf_t src).
+
+Copyright 2001-2004, 2017 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-impl.h"
+#include "longlong.h"
+
+
+double
+mpf_get_d_2exp (signed long int *expptr, mpf_srcptr src)
+{
+ mp_size_t size, abs_size;
+ mp_srcptr ptr;
+ int cnt;
+
+ size = SIZ(src);
+ if (UNLIKELY (size == 0))
+ {
+ *expptr = 0;
+ return 0.0;
+ }
+
+ ptr = PTR(src);
+ abs_size = ABS (size);
+ count_leading_zeros (cnt, ptr[abs_size - 1]);
+ cnt -= GMP_NAIL_BITS;
+
+ *expptr = EXP(src) * GMP_NUMB_BITS - cnt;
+ return mpn_get_d (ptr, abs_size, size, -(abs_size * GMP_NUMB_BITS - cnt));
+}
diff --git a/third_party/gmp/mpf/get_dfl_prec.c b/third_party/gmp/mpf/get_dfl_prec.c
new file mode 100644
index 0000000..13fc514
--- /dev/null
+++ b/third_party/gmp/mpf/get_dfl_prec.c
@@ -0,0 +1,38 @@
+/* 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-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/third_party/gmp/mpf/get_prc.c b/third_party/gmp/mpf/get_prc.c
new file mode 100644
index 0000000..8dee99e
--- /dev/null
+++ b/third_party/gmp/mpf/get_prc.c
@@ -0,0 +1,37 @@
+/* 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-impl.h"
+
+mp_bitcnt_t
+mpf_get_prec (mpf_srcptr x) __GMP_NOTHROW
+{
+ return __GMPF_PREC_TO_BITS (x->_mp_prec);
+}
diff --git a/third_party/gmp/mpf/get_si.c b/third_party/gmp/mpf/get_si.c
new file mode 100644
index 0000000..6ac4d44
--- /dev/null
+++ b/third_party/gmp/mpf/get_si.c
@@ -0,0 +1,86 @@
+/* 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-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/third_party/gmp/mpf/get_str.c b/third_party/gmp/mpf/get_str.c
new file mode 100644
index 0000000..946c4ae
--- /dev/null
+++ b/third_party/gmp/mpf/get_str.c
@@ -0,0 +1,320 @@
+/* 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, 2015, 2017 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-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);
+
+ num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
+ if (base > 1)
+ {
+ if (base <= 36)
+ num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz";
+ else if (UNLIKELY (base > 62))
+ return NULL;
+ }
+ else if (base > -2)
+ {
+ base = 10;
+ }
+ else
+ {
+ base = -base;
+ if (UNLIKELY (base > 36))
+ return NULL;
+ }
+
+ 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 = __GMP_ALLOCATE_FUNC_TYPE (n_digits + 2, char);
+ }
+
+ 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 (un > n_limbs_needed)
+ {
+ up += un - n_limbs_needed;
+ un = n_limbs_needed;
+ }
+
+ TMP_ALLOC_LIMBS_2 (pp, 2 * n_limbs_needed + 4,
+ tp, 2 * n_limbs_needed + 4);
+
+ 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);
+
+ pn = mpn_pow_1_highpart (pp, &ign, (mp_limb_t) base, e, n_limbs_needed + 1, 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);
+
+ pn = mpn_pow_1_highpart (pp, &ign, (mp_limb_t) base, e, n_limbs_needed + 1, 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/third_party/gmp/mpf/get_ui.c b/third_party/gmp/mpf/get_ui.c
new file mode 100644
index 0000000..e7b9333
--- /dev/null
+++ b/third_party/gmp/mpf/get_ui.c
@@ -0,0 +1,101 @@
+/* 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-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/third_party/gmp/mpf/init.c b/third_party/gmp/mpf/init.c
new file mode 100644
index 0000000..26ab262
--- /dev/null
+++ b/third_party/gmp/mpf/init.c
@@ -0,0 +1,41 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+}
diff --git a/third_party/gmp/mpf/init2.c b/third_party/gmp/mpf/init2.c
new file mode 100644
index 0000000..b90a08a
--- /dev/null
+++ b/third_party/gmp/mpf/init2.c
@@ -0,0 +1,43 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+}
diff --git a/third_party/gmp/mpf/inits.c b/third_party/gmp/mpf/inits.c
new file mode 100644
index 0000000..b6d054f
--- /dev/null
+++ b/third_party/gmp/mpf/inits.c
@@ -0,0 +1,49 @@
+/* mpf_inits() -- Initialize multiple mpf_t variables and set them to 0.
+
+Copyright 2009, 2015 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 "gmp-impl.h"
+
+void
+mpf_inits (mpf_ptr x, ...)
+{
+ va_list ap;
+
+ va_start (ap, x);
+
+ do
+ {
+ mpf_init (x);
+ x = va_arg (ap, mpf_ptr);
+ }
+ while (x != NULL);
+
+ va_end (ap);
+}
diff --git a/third_party/gmp/mpf/inp_str.c b/third_party/gmp/mpf/inp_str.c
new file mode 100644
index 0000000..c661a79
--- /dev/null
+++ b/third_party/gmp/mpf/inp_str.c
@@ -0,0 +1,92 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_TYPE (alloc_size, char);
+ 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 = __GMP_REALLOCATE_FUNC_TYPE (str, old_alloc_size, alloc_size, char);
+ }
+ 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 = __GMP_REALLOCATE_FUNC_TYPE (str, old_alloc_size, alloc_size, char);
+ }
+ 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/third_party/gmp/mpf/int_p.c b/third_party/gmp/mpf/int_p.c
new file mode 100644
index 0000000..024cfb5
--- /dev/null
+++ b/third_party/gmp/mpf/int_p.c
@@ -0,0 +1,55 @@
+/* mpf_integer_p -- test whether an mpf is an integer */
+
+/*
+Copyright 2001, 2002, 2014-2015 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-impl.h"
+
+
+int
+mpf_integer_p (mpf_srcptr f) __GMP_NOTHROW
+{
+ mp_srcptr fp;
+ mp_exp_t exp;
+ mp_size_t size;
+
+ size = SIZ (f);
+ exp = EXP (f);
+ if (exp <= 0)
+ return (size == 0); /* zero is an integer,
+ others have only fraction limbs */
+ size = ABS (size);
+
+ /* Ignore zeroes at the low end of F. */
+ for (fp = PTR (f); *fp == 0; ++fp)
+ --size;
+
+ /* no fraction limbs */
+ return size <= exp;
+}
diff --git a/third_party/gmp/mpf/iset.c b/third_party/gmp/mpf/iset.c
new file mode 100644
index 0000000..07f9006
--- /dev/null
+++ b/third_party/gmp/mpf/iset.c
@@ -0,0 +1,61 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+ 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/third_party/gmp/mpf/iset_d.c b/third_party/gmp/mpf/iset_d.c
new file mode 100644
index 0000000..2f36240
--- /dev/null
+++ b/third_party/gmp/mpf/iset_d.c
@@ -0,0 +1,41 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+
+ mpf_set_d (r, val);
+}
diff --git a/third_party/gmp/mpf/iset_si.c b/third_party/gmp/mpf/iset_si.c
new file mode 100644
index 0000000..65abb9a
--- /dev/null
+++ b/third_party/gmp/mpf/iset_si.c
@@ -0,0 +1,57 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+
+ 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/third_party/gmp/mpf/iset_str.c b/third_party/gmp/mpf/iset_str.c
new file mode 100644
index 0000000..10acda9
--- /dev/null
+++ b/third_party/gmp/mpf/iset_str.c
@@ -0,0 +1,43 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+
+ return mpf_set_str (r, s, base);
+}
diff --git a/third_party/gmp/mpf/iset_ui.c b/third_party/gmp/mpf/iset_ui.c
new file mode 100644
index 0000000..2c426bf
--- /dev/null
+++ b/third_party/gmp/mpf/iset_ui.c
@@ -0,0 +1,52 @@
+/* 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-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 = __GMP_ALLOCATE_FUNC_LIMBS (prec + 1);
+ 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/third_party/gmp/mpf/mul.c b/third_party/gmp/mpf/mul.c
new file mode 100644
index 0000000..3099461
--- /dev/null
+++ b/third_party/gmp/mpf/mul.c
@@ -0,0 +1,134 @@
+/* mpf_mul -- Multiply two floats.
+
+Copyright 1993, 1994, 1996, 2001, 2005, 2019 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-impl.h"
+
+void
+mpf_mul (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
+{
+ mp_size_t sign_product;
+ mp_size_t prec = r->_mp_prec;
+ mp_size_t rsize;
+ mp_limb_t cy_limb;
+ mp_ptr rp, tp;
+ mp_size_t adj;
+ TMP_DECL;
+
+ if (u == v)
+ {
+ mp_srcptr up;
+ mp_size_t usize;
+
+ usize = u->_mp_size;
+ sign_product = 0;
+
+ usize = ABS (usize);
+
+ up = u->_mp_d;
+ if (usize > prec)
+ {
+ up += usize - prec;
+ usize = prec;
+ }
+
+ if (usize == 0)
+ {
+ r->_mp_size = 0;
+ r->_mp_exp = 0; /* ??? */
+ return;
+ }
+ else
+ {
+ TMP_MARK;
+ rsize = 2 * usize;
+ tp = TMP_ALLOC_LIMBS (rsize);
+
+ mpn_sqr (tp, up, usize);
+ cy_limb = tp[rsize - 1];
+ }
+ }
+ else
+ {
+ mp_srcptr up, vp;
+ mp_size_t usize, vsize;
+
+ 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;
+ return;
+ }
+ else
+ {
+ TMP_MARK;
+ 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/third_party/gmp/mpf/mul_2exp.c b/third_party/gmp/mpf/mul_2exp.c
new file mode 100644
index 0000000..5de7363
--- /dev/null
+++ b/third_party/gmp/mpf/mul_2exp.c
@@ -0,0 +1,132 @@
+/* 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-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/third_party/gmp/mpf/mul_ui.c b/third_party/gmp/mpf/mul_ui.c
new file mode 100644
index 0000000..30da6ae
--- /dev/null
+++ b/third_party/gmp/mpf/mul_ui.c
@@ -0,0 +1,181 @@
+/* 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-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/third_party/gmp/mpf/neg.c b/third_party/gmp/mpf/neg.c
new file mode 100644
index 0000000..d294815
--- /dev/null
+++ b/third_party/gmp/mpf/neg.c
@@ -0,0 +1,61 @@
+/* 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-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/third_party/gmp/mpf/out_str.c b/third_party/gmp/mpf/out_str.c
new file mode 100644
index 0000000..1802d0f
--- /dev/null
+++ b/third_party/gmp/mpf/out_str.c
@@ -0,0 +1,116 @@
+/* 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-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/third_party/gmp/mpf/pow_ui.c b/third_party/gmp/mpf/pow_ui.c
new file mode 100644
index 0000000..8d54dc0
--- /dev/null
+++ b/third_party/gmp/mpf/pow_ui.c
@@ -0,0 +1,83 @@
+/* mpf_pow_ui -- Compute b^e.
+
+Copyright 1998, 1999, 2001, 2012, 2015 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-impl.h"
+#include "longlong.h"
+
+/* This uses a plain right-to-left square-and-multiply algorithm.
+
+ FIXME: When popcount(e) is not too small, it would probably speed things up
+ to use a k-ary sliding window algorithm. */
+
+void
+mpf_pow_ui (mpf_ptr r, mpf_srcptr b, unsigned long int e)
+{
+ mpf_t t;
+ int cnt;
+
+ if (e <= 1)
+ {
+ if (e == 0)
+ mpf_set_ui (r, 1);
+ else
+ mpf_set (r, b);
+ return;
+ }
+
+ count_leading_zeros (cnt, (mp_limb_t) e);
+ cnt = GMP_LIMB_BITS - 1 - cnt;
+
+ /* Increase computation precision as a function of the exponent. Adding
+ log2(popcount(e) + log2(e)) bits should be sufficient, but we add log2(e),
+ i.e. much more. With mpf's rounding of precision to whole limbs, this
+ will be excessive only when limbs are artificially small. */
+ mpf_init2 (t, mpf_get_prec (r) + cnt);
+
+ mpf_set (t, b); /* consume most significant bit */
+ while (--cnt > 0)
+ {
+ mpf_mul (t, t, t);
+ if ((e >> cnt) & 1)
+ mpf_mul (t, t, b);
+ }
+
+ /* Do the last iteration specially in order to save a copy operation. */
+ if (e & 1)
+ {
+ mpf_mul (t, t, t);
+ mpf_mul (r, t, b);
+ }
+ else
+ {
+ mpf_mul (r, t, t);
+ }
+
+ mpf_clear (t);
+}
diff --git a/third_party/gmp/mpf/random2.c b/third_party/gmp/mpf/random2.c
new file mode 100644
index 0000000..2e0163c
--- /dev/null
+++ b/third_party/gmp/mpf/random2.c
@@ -0,0 +1,66 @@
+/* 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-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/third_party/gmp/mpf/reldiff.c b/third_party/gmp/mpf/reldiff.c
new file mode 100644
index 0000000..81f05dd
--- /dev/null
+++ b/third_party/gmp/mpf/reldiff.c
@@ -0,0 +1,64 @@
+/* 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-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/third_party/gmp/mpf/set.c b/third_party/gmp/mpf/set.c
new file mode 100644
index 0000000..382fe86
--- /dev/null
+++ b/third_party/gmp/mpf/set.c
@@ -0,0 +1,55 @@
+/* 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-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/third_party/gmp/mpf/set_d.c b/third_party/gmp/mpf/set_d.c
new file mode 100644
index 0000000..0442f2f
--- /dev/null
+++ b/third_party/gmp/mpf/set_d.c
@@ -0,0 +1,59 @@
+/* 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-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/third_party/gmp/mpf/set_dfl_prec.c b/third_party/gmp/mpf/set_dfl_prec.c
new file mode 100644
index 0000000..9be71c0
--- /dev/null
+++ b/third_party/gmp/mpf/set_dfl_prec.c
@@ -0,0 +1,39 @@
+/* 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-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/third_party/gmp/mpf/set_prc.c b/third_party/gmp/mpf/set_prc.c
new file mode 100644
index 0000000..40c3f0e
--- /dev/null
+++ b/third_party/gmp/mpf/set_prc.c
@@ -0,0 +1,68 @@
+/* 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-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/third_party/gmp/mpf/set_prc_raw.c b/third_party/gmp/mpf/set_prc_raw.c
new file mode 100644
index 0000000..e5c52cc
--- /dev/null
+++ b/third_party/gmp/mpf/set_prc_raw.c
@@ -0,0 +1,39 @@
+/* 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-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/third_party/gmp/mpf/set_q.c b/third_party/gmp/mpf/set_q.c
new file mode 100644
index 0000000..b0b7fe3
--- /dev/null
+++ b/third_party/gmp/mpf/set_q.c
@@ -0,0 +1,118 @@
+/* mpf_set_q (mpf_t rop, mpq_t op) -- Convert the rational op to the float rop.
+
+Copyright 1996, 1999, 2001, 2002, 2004, 2005, 2016 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-impl.h"
+
+
+/* As usual the aim is to produce PREC(r) limbs, with the high non-zero. The
+ basic mpn_div_q 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 mpn_div_q.
+
+ 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. */
+
+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;
+ 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; /* size of intermediate numerator */
+ tp = TMP_ALLOC_LIMBS (tsize + 1); /* +1 for mpn_div_q's scratch */
+
+ if (zeros > 0)
+ {
+ /* pad n with zeros into temporary space */
+ MPN_ZERO (tp, zeros);
+ MPN_COPY (tp+zeros, np, nsize);
+ np = tp; /* mpn_div_q allows this overlap */
+ }
+ else
+ {
+ /* shorten n to get desired qsize */
+ np -= zeros;
+ }
+
+ ASSERT (tsize-dsize+1 == qsize);
+ mpn_div_q (qp, np, tsize, dp, dsize, tp);
+
+ /* 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/third_party/gmp/mpf/set_si.c b/third_party/gmp/mpf/set_si.c
new file mode 100644
index 0000000..23f713d
--- /dev/null
+++ b/third_party/gmp/mpf/set_si.c
@@ -0,0 +1,52 @@
+/* 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-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/third_party/gmp/mpf/set_str.c b/third_party/gmp/mpf/set_str.c
new file mode 100644
index 0000000..c7bfe0b
--- /dev/null
+++ b/third_party/gmp/mpf/set_str.c
@@ -0,0 +1,412 @@
+/* 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, 2019 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-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;
+ const char *expptr;
+ int exp_base;
+ const char *point = GMP_DECIMAL_POINT;
+ size_t pointlen = strlen (point);
+ const unsigned char *digit_value;
+ int incr;
+ size_t n_zeros_skipped;
+
+ 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);
+
+ incr = 0;
+ n_zeros_skipped = 0;
+ dotpos = NULL;
+
+ /* 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;
+ incr |= dig != 0;
+ s += incr; /* Increment after first non-0 digit seen. */
+ if (dotpos != NULL)
+ /* Count skipped zeros between radix point and first non-0
+ digit. */
+ n_zeros_skipped += 1 - incr;
+ }
+ }
+ 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 needs careful testing. Leave disabled for now. */
+ /* 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
+
+ if (str_size == 0)
+ {
+ SIZ(x) = 0;
+ EXP(x) = 0;
+ TMP_FREE;
+ return 0;
+ }
+
+ 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);
+
+ 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 + n_zeros_skipped;
+ 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);
+ TMP_ALLOC_LIMBS_2 (rp, ra, tp, 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/third_party/gmp/mpf/set_ui.c b/third_party/gmp/mpf/set_ui.c
new file mode 100644
index 0000000..bd4ba26
--- /dev/null
+++ b/third_party/gmp/mpf/set_ui.c
@@ -0,0 +1,48 @@
+/* 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-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/third_party/gmp/mpf/set_z.c b/third_party/gmp/mpf/set_z.c
new file mode 100644
index 0000000..f762633
--- /dev/null
+++ b/third_party/gmp/mpf/set_z.c
@@ -0,0 +1,56 @@
+/* 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-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/third_party/gmp/mpf/size.c b/third_party/gmp/mpf/size.c
new file mode 100644
index 0000000..f7a9dbd
--- /dev/null
+++ b/third_party/gmp/mpf/size.c
@@ -0,0 +1,38 @@
+/* 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-impl.h"
+
+size_t
+mpf_size (mpf_srcptr f) __GMP_NOTHROW
+{
+ return __GMP_ABS (f->_mp_size);
+}
diff --git a/third_party/gmp/mpf/sqrt.c b/third_party/gmp/mpf/sqrt.c
new file mode 100644
index 0000000..ffb7c10
--- /dev/null
+++ b/third_party/gmp/mpf/sqrt.c
@@ -0,0 +1,112 @@
+/* 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-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/third_party/gmp/mpf/sqrt_ui.c b/third_party/gmp/mpf/sqrt_ui.c
new file mode 100644
index 0000000..9f91f99
--- /dev/null
+++ b/third_party/gmp/mpf/sqrt_ui.c
@@ -0,0 +1,108 @@
+/* mpf_sqrt_ui -- Compute the square root of an unsigned integer.
+
+Copyright 1993, 1994, 1996, 2000, 2001, 2004, 2005, 2015 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-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 <= 1))
+ {
+ SIZ (r) = EXP (r) = u;
+ *PTR (r) = u;
+ return;
+ }
+
+ TMP_MARK;
+
+ prec = PREC (r);
+ 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 (PTR (r), NULL, tp, rsize);
+
+ SIZ (r) = prec;
+ EXP (r) = 1;
+ TMP_FREE;
+}
diff --git a/third_party/gmp/mpf/sub.c b/third_party/gmp/mpf/sub.c
new file mode 100644
index 0000000..56f26f6
--- /dev/null
+++ b/third_party/gmp/mpf/sub.c
@@ -0,0 +1,395 @@
+/* mpf_sub -- Subtract two floats.
+
+Copyright 1993-1996, 1999-2002, 2004, 2005, 2011, 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 "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 = SIZ (u);
+ vsize = SIZ (v);
+
+ /* 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 = EXP (v);
+ v_negated._mp_d = PTR (v);
+ 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 (EXP (u) < EXP (v))
+ {
+ mpf_srcptr t;
+ t = u; u = v; v = t;
+ negate ^= 1;
+ usize = SIZ (u);
+ vsize = SIZ (v);
+ }
+
+ usize = ABS (usize);
+ vsize = ABS (vsize);
+ up = PTR (u);
+ vp = PTR (v);
+ rp = PTR (r);
+ prec = PREC (r) + 1;
+ exp = EXP (u);
+ ediff = exp - EXP (v);
+
+ /* 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. */
+ /* This loop normally exits immediately. Optimize for that. */
+ while (up[usize - 1] == vp[vsize - 1])
+ {
+ 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;
+ }
+ }
+
+ 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--;
+ }
+ }
+ else 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)
+ {
+ MPN_COPY (tp, up, usize);
+ tp[usize] = 1;
+ rsize = usize + 1;
+ exp++;
+ goto normalized;
+ }
+ if (usize == 0)
+ {
+ cy_limb = mpn_neg (tp, vp, vsize);
+ rsize = vsize;
+ }
+ else 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;
+ size = vsize - usize;
+ cy_limb = mpn_neg (tp, vp, size);
+ cy_limb = mpn_sub_nc (tp + size, up, vp + size, usize, cy_limb);
+ rsize = vsize;
+ }
+ if (cy_limb == 0)
+ {
+ tp[rsize] = 1;
+ rsize++;
+ exp++;
+ goto normalized;
+ }
+ 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;
+ size = vsize - usize;
+ ASSERT_CARRY (mpn_neg (tp, vp, size));
+ mpn_sub_nc (tp + size, up, vp + size, usize, CNST_LIMB (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;
+ rsize = vsize + ediff;
+ size = rsize - usize;
+ ASSERT_CARRY (mpn_neg (tp, vp, size));
+ mpn_sub (tp + size, up, usize, vp + size, usize - ediff);
+ /* Should we use sub_nc then sub_1? */
+ MPN_DECR_U (tp + size, usize, CNST_LIMB (1));
+ }
+ }
+ }
+ else
+ {
+ /* uuuu */
+ /* vv */
+ mp_size_t size, i;
+ size = vsize + ediff - usize;
+ ASSERT_CARRY (mpn_neg (tp, vp, vsize));
+ 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--;
+ }
+ normalized:
+ MPN_COPY (rp, tp, rsize);
+ }
+
+ done:
+ TMP_FREE;
+ if (rsize == 0) {
+ SIZ (r) = 0;
+ EXP (r) = 0;
+ } else {
+ SIZ (r) = negate ? -rsize : rsize;
+ EXP (r) = exp;
+ }
+}
diff --git a/third_party/gmp/mpf/sub_ui.c b/third_party/gmp/mpf/sub_ui.c
new file mode 100644
index 0000000..a23d2a8
--- /dev/null
+++ b/third_party/gmp/mpf/sub_ui.c
@@ -0,0 +1,50 @@
+/* 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-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/third_party/gmp/mpf/swap.c b/third_party/gmp/mpf/swap.c
new file mode 100644
index 0000000..80b2e9b
--- /dev/null
+++ b/third_party/gmp/mpf/swap.c
@@ -0,0 +1,56 @@
+/* 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-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/third_party/gmp/mpf/trunc.c b/third_party/gmp/mpf/trunc.c
new file mode 100644
index 0000000..e9af4a7
--- /dev/null
+++ b/third_party/gmp/mpf/trunc.c
@@ -0,0 +1,74 @@
+/* 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-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/third_party/gmp/mpf/ui_div.c b/third_party/gmp/mpf/ui_div.c
new file mode 100644
index 0000000..d228bd4
--- /dev/null
+++ b/third_party/gmp/mpf/ui_div.c
@@ -0,0 +1,127 @@
+/* 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-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/third_party/gmp/mpf/ui_sub.c b/third_party/gmp/mpf/ui_sub.c
new file mode 100644
index 0000000..58da56b
--- /dev/null
+++ b/third_party/gmp/mpf/ui_sub.c
@@ -0,0 +1,281 @@
+/* mpf_ui_sub -- Subtract a float from an unsigned long int.
+
+Copyright 1993-1996, 2001, 2002, 2005, 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 "gmp-impl.h"
+
+void
+mpf_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v)
+{
+#if 1
+ __mpf_struct uu;
+ mp_limb_t ul;
+
+ if (u == 0)
+ {
+ mpf_neg (r, v);
+ return;
+ }
+
+ ul = u;
+ uu._mp_size = 1;
+ uu._mp_d = &ul;
+ uu._mp_exp = 1;
+ mpf_sub (r, &uu, v);
+
+#else
+ 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;
+ }
+
+ /* Signs are now known to be the same. */
+ ASSERT (vsize > 0);
+ ulimb = u;
+ /* Make U be the operand with the largest exponent. */
+ negate = 1 < v->_mp_exp;
+ prec = r->_mp_prec + negate;
+ rp = r->_mp_d;
+ if (negate)
+ {
+ usize = vsize;
+ vsize = 1;
+ up = v->_mp_d;
+ vp = &ulimb;
+ uexp = v->_mp_exp;
+ ediff = uexp - 1;
+
+ /* If U extends beyond PREC, ignore the part that does. */
+ if (usize > prec)
+ {
+ up += usize - prec;
+ usize = prec;
+ }
+ ASSERT (ediff > 0);
+ }
+ else
+ {
+ vp = v->_mp_d;
+ 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 && ulimb == vp[vsize - 1])
+ {
+ usize = 0;
+ vsize--;
+ uexp = 0;
+ /* Note that V might now have leading zero limbs.
+ In that case we have to adjust uexp. */
+ for (;;)
+ {
+ if (vsize == 0) {
+ rsize = 0;
+ uexp = 0;
+ goto done;
+ }
+ if ( vp[vsize - 1] != 0)
+ break;
+ vsize--, uexp--;
+ }
+ }
+ else
+ {
+ usize = 1;
+ uexp = 1;
+ up = &ulimb;
+ }
+ ASSERT (usize <= prec);
+ }
+
+ if (ediff >= prec)
+ {
+ /* V completely cancelled. */
+ if (rp != up)
+ MPN_COPY (rp, up, usize);
+ rsize = usize;
+ }
+ else
+ {
+ /* If V extends beyond PREC, ignore the part that does.
+ Note that this can make vsize neither zero nor negative. */
+ if (vsize + ediff > prec)
+ {
+ vp += vsize + ediff - prec;
+ vsize = prec - ediff;
+ }
+
+ /* Locate the least significant non-zero limb in (the needed
+ parts of) U and V, to simplify the code below. */
+ ASSERT (vsize > 0);
+ for (;;)
+ {
+ if (vp[0] != 0)
+ break;
+ vp++, vsize--;
+ if (vsize == 0)
+ {
+ MPN_COPY (rp, up, usize);
+ rsize = usize;
+ goto done;
+ }
+ }
+ for (;;)
+ {
+ if (usize == 0)
+ {
+ MPN_COPY (rp, vp, vsize);
+ rsize = vsize;
+ negate ^= 1;
+ goto done;
+ }
+ if (up[0] != 0)
+ break;
+ up++, usize--;
+ }
+
+ ASSERT (usize > 0 && vsize > 0);
+ TMP_MARK;
+
+ tp = TMP_ALLOC_LIMBS (prec);
+
+ /* uuuu | uuuu | uuuu | uuuu | uuuu */
+ /* vvvvvvv | vv | vvvvv | v | vv */
+
+ if (usize > ediff)
+ {
+ /* U and V partially overlaps. */
+ if (ediff == 0)
+ {
+ ASSERT (usize == 1 && vsize >= 1 && ulimb == *up); /* usize is 1>ediff, vsize >= 1 */
+ if (1 < vsize)
+ {
+ /* u */
+ /* vvvvvvv */
+ rsize = vsize;
+ vsize -= 1;
+ /* mpn_cmp (up, vp + vsize - usize, usize) > 0 */
+ if (ulimb > vp[vsize])
+ {
+ tp[vsize] = ulimb - vp[vsize] - 1;
+ ASSERT_CARRY (mpn_neg (tp, vp, vsize));
+ }
+ else
+ {
+ /* vvvvvvv */ /* Swap U and V. */
+ /* u */
+ MPN_COPY (tp, vp, vsize);
+ tp[vsize] = vp[vsize] - ulimb;
+ negate = 1;
+ }
+ }
+ else /* vsize == usize == 1 */
+ {
+ /* u */
+ /* v */
+ rsize = 1;
+ negate = ulimb < vp[0];
+ tp[0] = negate ? vp[0] - ulimb: ulimb - vp[0];
+ }
+ }
+ else
+ {
+ ASSERT (vsize + ediff <= usize);
+ ASSERT (vsize == 1 && usize >= 2 && ulimb == *vp);
+ {
+ /* uuuu */
+ /* v */
+ mp_size_t size;
+ size = usize - ediff - 1;
+ MPN_COPY (tp, up, size);
+ ASSERT_NOCARRY (mpn_sub_1 (tp + size, up + size, usize - size, ulimb));
+ rsize = usize;
+ }
+ /* Other cases are not possible */
+ /* uuuu */
+ /* vvvvv */
+ }
+ }
+ else
+ {
+ /* uuuu */
+ /* vv */
+ mp_size_t size, i;
+ ASSERT_CARRY (mpn_neg (tp, vp, vsize));
+ rsize = vsize + ediff;
+ size = rsize - usize;
+ for (i = vsize; i < size; i++)
+ tp[i] = GMP_NUMB_MAX;
+ ASSERT_NOCARRY (mpn_sub_1 (tp + size, up, usize, CNST_LIMB (1)));
+ }
+
+ /* Full normalize. Optimize later. */
+ while (rsize != 0 && tp[rsize - 1] == 0)
+ {
+ rsize--;
+ uexp--;
+ }
+ MPN_COPY (rp, tp, rsize);
+ TMP_FREE;
+ }
+
+ done:
+ r->_mp_size = negate ? -rsize : rsize;
+ r->_mp_exp = uexp;
+#endif
+}
diff --git a/third_party/gmp/mpf/urandomb.c b/third_party/gmp/mpf/urandomb.c
new file mode 100644
index 0000000..ebe85ab
--- /dev/null
+++ b/third_party/gmp/mpf/urandomb.c
@@ -0,0 +1,68 @@
+/* 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-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;
+}