diff options
Diffstat (limited to 'gnu/lib/libgmp/mpf')
60 files changed, 6001 insertions, 0 deletions
diff --git a/gnu/lib/libgmp/mpf/Makefile.bsd-wrapper b/gnu/lib/libgmp/mpf/Makefile.bsd-wrapper new file mode 100644 index 00000000000..5e7a3f2db65 --- /dev/null +++ b/gnu/lib/libgmp/mpf/Makefile.bsd-wrapper @@ -0,0 +1,27 @@ +# $OpenBSD: Makefile.bsd-wrapper,v 1.1 1997/06/24 21:28:14 provos Exp $ + +GNUCFLAGS= CFLAGS="${CFLAGS} ${COPTS}" + +all: config.status + ${MAKE} ${GNUCFLAGS} BISON=yacc CC=${CC} LDFLAGS=${LDSTATIC} + +.FORCE: .IGNORE + +config: .FORCE + -rm -f config.cache + /bin/sh ${.CURDIR}/configure --with-gnu-as --with-gnu-ld \ + --prefix=/usr --local-prefix=/usr + +config.status: Makefile.in configure + /bin/sh ${.CURDIR}/configure --with-gnu-as --with-gnu-ld \ + --prefix=/usr --local-prefix=/usr && touch config.status + +clean cleandir: + -@if [ -e Makefile ]; then ${MAKE} distclean; fi + +depend: + # Nothing here so far... + +.include <bsd.obj.mk> +.include <bsd.subdir.mk> +.include <bsd.man.mk> diff --git a/gnu/lib/libgmp/mpf/Makefile.in b/gnu/lib/libgmp/mpf/Makefile.in new file mode 100644 index 00000000000..598fa72fc19 --- /dev/null +++ b/gnu/lib/libgmp/mpf/Makefile.in @@ -0,0 +1,115 @@ +# Makefile for GNU MP/mpf functions +# Copyright (C) 1991, 1993, 1994, 1996 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 the GNU Library General Public License as published by +# the Free Software Foundation; either version 2 of the License, or (at your +# option) any later version. + +# 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 Library General Public +# License for more details. + +# You should have received a copy of the GNU Library General Public License +# along with the GNU MP Library; see the file COPYING.LIB. If not, write to +# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +# MA 02111-1307, USA. + +srcdir = . + +CC = gcc + +CFLAGS = -g -O +AR = ar +AR_FLAGS = rc +SHELL = /bin/sh + +#### host and target specific makefile fragments come in here. +### + +MPF_SRCS = init.c init2.c set.c set_ui.c set_si.c set_str.c set_d.c set_z.c \ + iset.c iset_ui.c iset_si.c iset_str.c iset_d.c clear.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_ui.c cmp_si.c mul_2exp.c div_2exp.c abs.c neg.c set_q.c get_d.c \ + set_dfl_prec.c set_prc.c set_prc_raw.c get_prc.c ui_div.c sqrt_ui.c +MPF_OBJS = init.o init2.o set.o set_ui.o set_si.o set_str.o set_d.o set_z.o \ + iset.o iset_ui.o iset_si.o iset_str.o iset_d.o clear.o get_str.o \ + dump.o size.o eq.o reldiff.o sqrt.o random2.o inp_str.o out_str.o \ + add.o add_ui.o sub.o sub_ui.o ui_sub.o mul.o mul_ui.o div.o div_ui.o \ + cmp.o cmp_ui.o cmp_si.o mul_2exp.o div_2exp.o abs.o neg.o set_q.o get_d.o \ + set_dfl_prec.o set_prc.o set_prc_raw.o get_prc.o ui_div.o sqrt_ui.o + +LATER_OBJS = inp_raw.o out_raw.o random.o pow_ui.o fac_ui.o + +INCLUDES = -I. -I.. -I../mpn -I$(srcdir)/.. + +libmpf.a: Makefile $(MPF_OBJS) + rm -f $@ + $(AR) $(AR_FLAGS) $@ $(MPF_OBJS) + +.c.o: + $(CC) -c $(INCLUDES) $(CFLAGS) $(XCFLAGS) $< + +clean mostlyclean: + rm -f *.o libmpf.a + -cd tests; $(MAKE) $@ +distclean maintainer-clean: clean + rm -f Makefile config.status + -cd tests; $(MAKE) $@ + +Makefile: $(srcdir)/Makefile.in + $(SHELL) ./config.status + +H = $(srcdir)/../gmp.h $(srcdir)/../gmp-impl.h ../mpn/gmp-mparam.h + +abs.o: $(srcdir)/abs.c $(H) +add.o: $(srcdir)/add.c $(H) +add_ui.o: $(srcdir)/add_ui.c $(H) +clear.o: $(srcdir)/clear.c $(H) +cmp.o: $(srcdir)/cmp.c $(H) +cmp_si.o: $(srcdir)/cmp_si.c $(H) +cmp_ui.o: $(srcdir)/cmp_ui.c $(H) +eq.o: $(srcdir)/eq.c $(H) +div.o: $(srcdir)/div.c $(H) $(srcdir)/../longlong.h +div_2exp.o: $(srcdir)/div_2exp.c $(H) +div_ui.o: $(srcdir)/div_ui.c $(H) $(srcdir)/../longlong.h +dump.o: $(srcdir)/dump.c $(H) +get_d.o: $(srcdir)/get_d.c $(H) +get_prc.o: $(srcdir)/get_prc.c $(H) +get_str.o: $(srcdir)/get_str.c $(H) $(srcdir)/../longlong.h +init.o: $(srcdir)/init.c $(H) +init2.o: $(srcdir)/init2.c $(H) +inp_str.o: $(srcdir)/inp_str.c $(H) +iset.o: $(srcdir)/iset.c $(H) +iset_d.o: $(srcdir)/iset_d.c $(H) +iset_si.o: $(srcdir)/iset_si.c $(H) +iset_str.o: $(srcdir)/iset_str.c $(H) +iset_ui.o: $(srcdir)/iset_ui.c $(H) +mul.o: $(srcdir)/mul.c $(H) +mul_2exp.o: $(srcdir)/mul_2exp.c $(H) +mul_ui.o: $(srcdir)/mul_ui.c $(H) +neg.o: $(srcdir)/neg.c $(H) +out_str.o: $(srcdir)/out_str.c $(H) +random2.o: $(srcdir)/random2.c $(H) +reldiff.o: $(srcdir)/reldiff.c $(H) +set.o: $(srcdir)/set.c $(H) +set_d.o: $(srcdir)/set_d.c $(H) +set_dfl_prec.o: $(srcdir)/set_dfl_prec.c $(H) +set_prc.o: $(srcdir)/set_prc.c $(H) +set_prc_raw.o: $(srcdir)/set_prc_raw.c $(H) +set_q.o: $(srcdir)/set_q.c $(H) +set_si.o: $(srcdir)/set_si.c $(H) +set_str.o: $(srcdir)/set_str.c $(H) $(srcdir)/../longlong.h +set_ui.o: $(srcdir)/set_ui.c $(H) +set_z.o: $(srcdir)/set_z.c $(H) +size.o: $(srcdir)/size.c $(H) +sqrt.o: $(srcdir)/sqrt.c $(H) +sqrt_ui.o: $(srcdir)/sqrt_ui.c $(H) +sub.o: $(srcdir)/sub.c $(H) +sub_ui.o: $(srcdir)/sub_ui.c $(H) +ui_div.o: $(srcdir)/ui_div.c $(H) $(srcdir)/../longlong.h +ui_sub.o: $(srcdir)/ui_sub.c $(H) diff --git a/gnu/lib/libgmp/mpf/abs.c b/gnu/lib/libgmp/mpf/abs.c new file mode 100644 index 00000000000..029007a5f87 --- /dev/null +++ b/gnu/lib/libgmp/mpf/abs.c @@ -0,0 +1,56 @@ +/* mpf_abs -- Compute the absolute value of a float. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_abs (mpf_ptr r, mpf_srcptr u) +#else +mpf_abs (r, u) + mpf_ptr r; + mpf_srcptr u; +#endif +{ + 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/gnu/lib/libgmp/mpf/add.c b/gnu/lib/libgmp/mpf/add.c new file mode 100644 index 00000000000..2db876fc657 --- /dev/null +++ b/gnu/lib/libgmp/mpf/add.c @@ -0,0 +1,180 @@ +/* mpf_add -- Add two floats. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_add (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +#else +mpf_add (r, u, v) + mpf_ptr r; + mpf_srcptr u; + mpf_srcptr v; +#endif +{ + 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 (marker); + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (usize == 0) + { + mpf_set (r, v); + return; + } + if (vsize == 0) + { + mpf_set (r, u); + return; + } + + /* 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 (marker); + + /* 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 = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB); + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (tp != up) + MPN_COPY (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 (marker); +} diff --git a/gnu/lib/libgmp/mpf/add_ui.c b/gnu/lib/libgmp/mpf/add_ui.c new file mode 100644 index 00000000000..f181df48e62 --- /dev/null +++ b/gnu/lib/libgmp/mpf/add_ui.c @@ -0,0 +1,151 @@ +/* mpf_add_ui -- Add a float and an unsigned integer. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_add_ui (mpf_ptr sum, mpf_srcptr u, unsigned long int v) +#else +mpf_add_ui (sum, u, v) + mpf_ptr sum; + mpf_srcptr u; + unsigned long int v; +#endif +{ + 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 (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 (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/gnu/lib/libgmp/mpf/clear.c b/gnu/lib/libgmp/mpf/clear.c new file mode 100644 index 00000000000..beaf4ee3a3c --- /dev/null +++ b/gnu/lib/libgmp/mpf/clear.c @@ -0,0 +1,35 @@ +/* mpf_clear -- de-allocate the space occupied by the dynamic digit space of + an integer. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_clear (mpf_ptr m) +#else +mpf_clear (m) + mpf_ptr m; +#endif +{ + (*_mp_free_func) (m->_mp_d, (m->_mp_prec + 1) * BYTES_PER_MP_LIMB); +} diff --git a/gnu/lib/libgmp/mpf/cmp.c b/gnu/lib/libgmp/mpf/cmp.c new file mode 100644 index 00000000000..d440e1152bd --- /dev/null +++ b/gnu/lib/libgmp/mpf/cmp.c @@ -0,0 +1,114 @@ +/* mpf_cmp -- Compare two floats. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpf_cmp (mpf_srcptr u, mpf_srcptr v) +#else +mpf_cmp (u, v) + mpf_srcptr u; + mpf_srcptr v; +#endif +{ + mp_srcptr up, vp; + mp_size_t usize, vsize; + mp_exp_t uexp, vexp; + int cmp; + int usign; + + uexp = u->_mp_exp; + vexp = v->_mp_exp; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* 1. Are the signs different? */ + if ((usize ^ vsize) >= 0) + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + /* vsize >= 0 */ + return -(vsize != 0); + if (vsize == 0) + /* usize >= 0 */ + return usize != 0; + /* Fall out. */ + } + else + { + /* Either U or V is negative, but not both. */ + return usize >= 0 ? 1 : -1; + } + + /* U and V have the same sign and are both non-zero. */ + + usign = usize >= 0 ? 1 : -1; + + /* 2. Are the exponents different? */ + if (uexp > vexp) + return usign; + if (uexp < vexp) + return -usign; + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + +#define STRICT_MPF_NORMALIZATION 0 +#if ! STRICT_MPF_NORMALIZATION + /* Ignore zeroes at the low end of U and V. */ + while (up[0] == 0) + { + up++; + usize--; + } + while (vp[0] == 0) + { + vp++; + vsize--; + } +#endif + + if (usize > vsize) + { + cmp = mpn_cmp (up + usize - vsize, vp, vsize); + if (cmp == 0) + return usign; + } + else if (vsize > usize) + { + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp == 0) + return -usign; + } + else + { + cmp = mpn_cmp (up, vp, usize); + if (cmp == 0) + return 0; + } + return cmp > 0 ? usign : -usign; +} diff --git a/gnu/lib/libgmp/mpf/cmp_si.c b/gnu/lib/libgmp/mpf/cmp_si.c new file mode 100644 index 00000000000..01f970812c9 --- /dev/null +++ b/gnu/lib/libgmp/mpf/cmp_si.c @@ -0,0 +1,98 @@ +/* mpf_cmp_si -- Compare a float with a signed integer. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpf_cmp_si (mpf_srcptr u, long int vslimb) +#else +mpf_cmp_si (u, vslimb) + mpf_srcptr u; + long int vslimb; +#endif +{ + mp_srcptr up; + mp_size_t usize; + mp_exp_t uexp; + int usign; + + uexp = u->_mp_exp; + usize = u->_mp_size; + + /* 1. Are the signs different? */ + if ((usize < 0) == (vslimb < 0)) /* don't use xor, type size may differ */ + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + /* vslimb >= 0 */ + return -(vslimb != 0); + if (vslimb == 0) + /* usize >= 0 */ + return usize != 0; + /* Fall out. */ + } + else + { + /* Either U or V is negative, but not both. */ + return usize >= 0 ? 1 : -1; + } + + /* U and V have the same sign and are both non-zero. */ + + usign = usize >= 0 ? 1 : -1; + + /* 2. Are the exponents different (V's exponent == 1)? */ + if (uexp > 1) + return usign; + if (uexp < 1) + return -usign; + + usize = ABS (usize); + vslimb = ABS (vslimb); + + up = u->_mp_d; + +#define STRICT_MPF_NORMALIZATION 0 +#if ! STRICT_MPF_NORMALIZATION + /* Ignore zeroes at the low end of U and V. */ + while (*up == 0) + { + up++; + usize--; + } +#endif + + /* 3. 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 > 1) + return usign; + + /* 4. Compare the mantissas. */ + if (*up > vslimb) + return usign; + else if (*up < vslimb) + return -usign; + + /* Wow, we got zero even if we tried hard to avoid it. */ + return 0; +} diff --git a/gnu/lib/libgmp/mpf/cmp_ui.c b/gnu/lib/libgmp/mpf/cmp_ui.c new file mode 100644 index 00000000000..3a4911ba716 --- /dev/null +++ b/gnu/lib/libgmp/mpf/cmp_ui.c @@ -0,0 +1,80 @@ +/* mpf_cmp_ui -- Compare a float with an unsigned integer. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpf_cmp_ui (mpf_srcptr u, unsigned long int vlimb) +#else +mpf_cmp_ui (u, vlimb) + mpf_srcptr u; + unsigned long int vlimb; +#endif +{ + mp_srcptr up; + mp_size_t usize; + mp_exp_t uexp; + + uexp = u->_mp_exp; + usize = u->_mp_size; + + /* 1. Is U negative? */ + if (usize < 0) + return -1; + /* We rely on usize being non-negative in the code that follows. */ + + if (vlimb == 0) + return usize != 0; + + /* 2. Are the exponents different (V's exponent == 1)? */ + if (uexp > 1) + return 1; + if (uexp < 1) + return -1; + + up = u->_mp_d; + +#define STRICT_MPF_NORMALIZATION 0 +#if ! STRICT_MPF_NORMALIZATION + /* Ignore zeroes at the low end of U. */ + while (*up == 0) + { + up++; + usize--; + } +#endif + + /* 3. 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 > 1) + return 1; + + /* 4. Compare the mantissas. */ + if (*up > vlimb) + return 1; + else if (*up < vlimb) + return -1; + + /* Wow, we got zero even if we tried hard to avoid it. */ + return 0; +} diff --git a/gnu/lib/libgmp/mpf/configure.in b/gnu/lib/libgmp/mpf/configure.in new file mode 100644 index 00000000000..b6ecf8b3e8e --- /dev/null +++ b/gnu/lib/libgmp/mpf/configure.in @@ -0,0 +1,12 @@ +# This file is a shell script fragment that supplies the information +# necessary for a configure script to process the program in +# this directory. For more information, look at ../configure. + +configdirs=tests +srctrigger=add_ui.c +srcname="GNU Multi-Precision library/mpf" + +# per-host: + +# per-target: + diff --git a/gnu/lib/libgmp/mpf/div.c b/gnu/lib/libgmp/mpf/div.c new file mode 100644 index 00000000000..9ba6922cd63 --- /dev/null +++ b/gnu/lib/libgmp/mpf/div.c @@ -0,0 +1,144 @@ +/* mpf_div -- Divide two floats. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +#else +mpf_div (r, u, v) + mpf_ptr r; + mpf_srcptr u; + mpf_srcptr v; +#endif +{ + mp_srcptr up, vp; + mp_ptr rp, tp, rtp; + mp_size_t usize, vsize; + mp_size_t rsize, tsize; + mp_size_t sign_quotient; + mp_size_t prec; + unsigned normalization_steps; + mp_limb_t q_limb; + mp_exp_t rexp; + TMP_DECL (marker); + + usize = u->_mp_size; + vsize = v->_mp_size; + sign_quotient = usize ^ vsize; + usize = ABS (usize); + vsize = ABS (vsize); + prec = r->_mp_prec; + + if (vsize == 0) + vsize = 1 / vsize; /* divide by zero as directed */ + if (usize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK (marker); + rexp = u->_mp_exp - v->_mp_exp; + + rp = r->_mp_d; + up = u->_mp_d; + vp = v->_mp_d; + + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + + tsize = vsize + prec; + tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB); + + if (usize > tsize) + { + up += usize - tsize; + usize = tsize; + rtp = tp; + } + else + { + MPN_ZERO (tp, tsize - usize); + rtp = tp + (tsize - usize); + } + + count_leading_zeros (normalization_steps, vp[vsize - 1]); + + /* Normalize the divisor and the dividend. */ + if (normalization_steps != 0) + { + mp_ptr tmp; + mp_limb_t nlimb; + + /* Shift up the divisor setting the most significant bit of + the most significant limb. Use temporary storage not to clobber + the original contents of the divisor. */ + tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + mpn_lshift (tmp, vp, vsize, normalization_steps); + vp = tmp; + + /* Shift up the dividend, possibly introducing a new most + significant word. Move the shifted dividend in the remainder + at the same time. */ + nlimb = mpn_lshift (rtp, up, usize, normalization_steps); + if (nlimb != 0) + { + rtp[usize] = nlimb; + tsize++; + rexp++; + } + } + else + { + /* The divisor is already normalized, as required. + Copy it to temporary space if it overlaps with the quotient. */ + if (vp - rp <= tsize - vsize) + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + MPN_COPY (tmp, vp, vsize); + vp = (mp_srcptr) tmp; + } + + /* Move the dividend to the remainder. */ + MPN_COPY (rtp, up, usize); + } + + q_limb = mpn_divmod (rp, tp, tsize, vp, vsize); + rsize = tsize - vsize; + if (q_limb) + { + rp[rsize] = q_limb; + rsize++; + rexp++; + } + + r->_mp_size = sign_quotient >= 0 ? rsize : -rsize; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/div_2exp.c b/gnu/lib/libgmp/mpf/div_2exp.c new file mode 100644 index 00000000000..d7296254e9e --- /dev/null +++ b/gnu/lib/libgmp/mpf/div_2exp.c @@ -0,0 +1,79 @@ +/* mpf_div_2exp -- Divide a float by 2^n. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_div_2exp (mpf_ptr r, mpf_srcptr u, unsigned long int exp) +#else +mpf_div_2exp (r, u, exp) + mpf_ptr r; + mpf_srcptr u; + unsigned long int exp; +#endif +{ + 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 (usize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + abs_usize = ABS (usize); + up = u->_mp_d; + + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + } + + if (exp % BITS_PER_MP_LIMB == 0) + { + if (rp != up) + MPN_COPY (rp, up, abs_usize); + r->_mp_exp = uexp - exp / BITS_PER_MP_LIMB; + } + else + { + /* Use mpn_lshift since mpn_rshift operates upwards, and we therefore + would clobber part of U before using that part, when R == U. */ + mp_limb_t cy_limb; + cy_limb = mpn_lshift (rp, up, abs_usize, -exp % BITS_PER_MP_LIMB); + rp[abs_usize] = cy_limb; + cy_limb = cy_limb != 0; + + abs_usize += cy_limb; + r->_mp_exp = uexp - exp / BITS_PER_MP_LIMB - 1 + cy_limb; + } + r->_mp_size = usize >= 0 ? abs_usize : -abs_usize; +} diff --git a/gnu/lib/libgmp/mpf/div_ui.c b/gnu/lib/libgmp/mpf/div_ui.c new file mode 100644 index 00000000000..050a6f0b70e --- /dev/null +++ b/gnu/lib/libgmp/mpf/div_ui.c @@ -0,0 +1,91 @@ +/* mpf_div_ui -- Divide a float with an unsigned integer. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpf_div_ui (mpf_ptr r, mpf_srcptr u, unsigned long int v) +#else +mpf_div_ui (r, u, v) + mpf_ptr r; + mpf_srcptr u; + unsigned long int v; +#endif +{ + 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 (marker); + + usize = u->_mp_size; + sign_quotient = usize; + usize = ABS (usize); + prec = r->_mp_prec; + + if (v == 0) + v = 1 / v; /* divide by zero as directed */ + if (usize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK (marker); + + rp = r->_mp_d; + up = u->_mp_d; + + tsize = 1 + prec; + tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB); + + 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 (marker); +} diff --git a/gnu/lib/libgmp/mpf/dump.c b/gnu/lib/libgmp/mpf/dump.c new file mode 100644 index 00000000000..46d5c05f20a --- /dev/null +++ b/gnu/lib/libgmp/mpf/dump.c @@ -0,0 +1,43 @@ +/* mpf_dump -- Dump a float to stdout. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_dump (mpf_srcptr u) +#else +mpf_dump (u) + mpf_srcptr u; +#endif +{ + 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); + (*_mp_free_func) (str, 0);/* ??? broken alloc interface, pass what size ??? */ +} diff --git a/gnu/lib/libgmp/mpf/eq.c b/gnu/lib/libgmp/mpf/eq.c new file mode 100644 index 00000000000..e6630ad6e9f --- /dev/null +++ b/gnu/lib/libgmp/mpf/eq.c @@ -0,0 +1,121 @@ +/* mpf_eq -- Compare two floats up to a specified bit #. + +Copyright (C) 1993, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpf_eq (mpf_srcptr u, mpf_srcptr v, unsigned long int n_bits) +#else +mpf_eq (u, v, n_bits) + mpf_srcptr u; + mpf_srcptr v; + unsigned long int n_bits; +#endif +{ + mp_srcptr up, vp; + mp_size_t usize, vsize, size, i; + mp_exp_t uexp, vexp; + int usign; + + uexp = u->_mp_exp; + vexp = v->_mp_exp; + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* 1. Are the signs different? */ + if ((usize ^ vsize) >= 0) + { + /* U and V are both non-negative or both negative. */ + if (usize == 0) + 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. */ + + usign = usize >= 0 ? 1 : -1; + + /* 2. Are the exponents different? */ + if (uexp > vexp) + return 0; /* ??? handle (uexp = vexp + 1) */ + if (vexp > uexp) + return 0; /* ??? handle (vexp = uexp + 1) */ + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + + /* Ignore zeroes at the low end of U and V. */ + while (up[0] == 0) + { + up++; + usize--; + } + while (vp[0] == 0) + { + vp++; + vsize--; + } + + if (usize > vsize) + { + if (vsize * BITS_PER_MP_LIMB < n_bits) + return 0; /* surely too different */ + size = vsize; + } + else if (vsize > usize) + { + if (usize * BITS_PER_MP_LIMB < n_bits) + return 0; /* surely too different */ + size = usize; + } + else + { + size = usize; + } + + if (size > (n_bits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) + size = (n_bits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; + + up += usize - size; + vp += vsize - size; + + for (i = size - 1; i >= 0; i--) + { + if (up[i] != vp[i]) + return 0; + } + + return 1; +} diff --git a/gnu/lib/libgmp/mpf/get_d.c b/gnu/lib/libgmp/mpf/get_d.c new file mode 100644 index 00000000000..a10b366a952 --- /dev/null +++ b/gnu/lib/libgmp/mpf/get_d.c @@ -0,0 +1,54 @@ +/* double mpf_get_d (mpf_t src) -- Return the double approximation to SRC. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +double +#if __STDC__ +mpf_get_d (mpf_srcptr src) +#else +mpf_get_d (src) + mpf_srcptr src; +#endif +{ + double res; + mp_size_t size, i, n_limbs_to_use; + int negative; + mp_ptr qp; + + size = SIZ(src); + if (size == 0) + return 0.0; + + negative = size < 0; + size = ABS (size); + qp = PTR(src); + + res = qp[size - 1]; + n_limbs_to_use = MIN (LIMBS_PER_DOUBLE, size); + for (i = 2; i <= n_limbs_to_use; i++) + res = res * MP_BASE_AS_DOUBLE + qp[size - i]; + + res = __gmp_scale2 (res, (EXP(src) - n_limbs_to_use) * BITS_PER_MP_LIMB); + + return negative ? -res : res; +} diff --git a/gnu/lib/libgmp/mpf/get_prc.c b/gnu/lib/libgmp/mpf/get_prc.c new file mode 100644 index 00000000000..7f7e41f673f --- /dev/null +++ b/gnu/lib/libgmp/mpf/get_prc.c @@ -0,0 +1,34 @@ +/* mpf_get_prec(x) -- Return the precision in bits of x. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +unsigned long int +#if __STDC__ +mpf_get_prec (mpf_srcptr x) +#else +mpf_get_prec (x) + mpf_srcptr x; +#endif +{ + return (unsigned long int) x->_mp_prec * BITS_PER_MP_LIMB - BITS_PER_MP_LIMB; +} diff --git a/gnu/lib/libgmp/mpf/get_str.c b/gnu/lib/libgmp/mpf/get_str.c new file mode 100644 index 00000000000..bfee18d076b --- /dev/null +++ b/gnu/lib/libgmp/mpf/get_str.c @@ -0,0 +1,500 @@ +/* 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 (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +/* + New algorithm for converting fractions (951019): + 0. Call the fraction to convert F. + 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e., + [exp * BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is + the number of limbs between the limb point and the most significant + non-zero limb. Call this result n. + 2. Compute B^n. + 3. F*B^n will now be just below 1, which can be converted easily. (Just + multiply by B repeatedly, and see the digits fall out as integers.) + We should interrupt the conversion process of F*B^n as soon as the number + of digits requested have been generated. + + New algorithm for converting integers (951019): + 0. Call the integer to convert I. + 1. Compute [exp * log(2^BITS_PER_MP_LIMB)/log(B)], i.e., + [exp BITS_PER_MP_LIMB * __mp_bases[B].chars_per_bit_exactly]. Exp is + the number of limbs between the limb point and the least significant + non-zero limb. Call this result n. + 2. Compute B^n. + 3. I/B^n can be converted easily. (Just divide by B repeatedly. In GMP, + this is best done by calling mpn_get_str.) + Note that converting I/B^n could yield more digits than requested. For + efficiency, the variable n above should be set larger in such cases, to + kill all undesired digits in the division in step 3. +*/ + +char * +#if __STDC__ +mpf_get_str (char *digit_ptr, mp_exp_t *exp, int base, size_t n_digits, mpf_srcptr u) +#else +mpf_get_str (digit_ptr, exp, base, n_digits, u) + char *digit_ptr; + mp_exp_t *exp; + int base; + size_t n_digits; + mpf_srcptr u; +#endif +{ + mp_size_t usize; + mp_exp_t uexp; + unsigned char *str; + size_t str_size; + char *num_to_text; + long i; /* should be size_t */ + mp_ptr rp; + mp_limb_t big_base; + size_t digits_computed_so_far; + int dig_per_u; + mp_srcptr up; + unsigned char *tstr; + mp_exp_t exp_in_base; + TMP_DECL (marker); + + TMP_MARK (marker); + usize = u->_mp_size; + uexp = u->_mp_exp; + + if (base >= 0) + { + if (base == 0) + base = 10; + num_to_text = "0123456789abcdefghijklmnopqrstuvwxyz"; + } + else + { + base = -base; + num_to_text = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; + } + + /* Don't compute more digits than U can accurately represent. + Also, if 0 digits were requested, give *exactly* as many digits + as can be accurately represented. */ + { + size_t max_digits = (((u->_mp_prec - 1) * BITS_PER_MP_LIMB) + * __mp_bases[base].chars_per_bit_exactly); + if (n_digits == 0 || n_digits > max_digits) + n_digits = max_digits; + } + + if (digit_ptr == 0) + { + /* We didn't get a string from the user. Allocate one (and return + a pointer to it) with space for `-' and terminating null. */ + digit_ptr = (char *) (*_mp_allocate_func) (n_digits + 2); + } + + if (usize == 0) + { + *exp = 0; + *digit_ptr = 0; + return digit_ptr; + } + + str = (unsigned char *) digit_ptr; + + /* Allocate temporary digit space. We can't put digits directly in the user + area, since we almost always generate more digits than requested. */ + tstr = (unsigned char *) TMP_ALLOC (n_digits + 3 * BITS_PER_MP_LIMB); + + if (usize < 0) + { + *digit_ptr = '-'; + str++; + usize = -usize; + } + + digits_computed_so_far = 0; + + if (uexp > usize) + { + /* The number has just an integral part. */ + mp_size_t rsize; + mp_size_t exp_in_limbs; + mp_size_t msize; + mp_ptr tp, xp, mp; + int cnt; + mp_limb_t cy; + mp_size_t start_str; + mp_size_t n_limbs; + + n_limbs = 2 + ((mp_size_t) (n_digits / __mp_bases[base].chars_per_bit_exactly) + / BITS_PER_MP_LIMB); + + /* Compute n such that [u/B^n] contains (somewhat) more than n_digits + digits. (We compute less than that only if that is an exact number, + i.e., exp is small enough.) */ + + exp_in_limbs = uexp; + + if (n_limbs >= exp_in_limbs) + { + /* The number is so small that we convert the entire number. */ + exp_in_base = 0; + rp = (mp_ptr) TMP_ALLOC (exp_in_limbs * BYTES_PER_MP_LIMB); + MPN_ZERO (rp, exp_in_limbs - usize); + MPN_COPY (rp + (exp_in_limbs - usize), u->_mp_d, usize); + rsize = exp_in_limbs; + } + else + { + exp_in_limbs -= n_limbs; + exp_in_base = (((exp_in_limbs * BITS_PER_MP_LIMB - 1)) + * __mp_bases[base].chars_per_bit_exactly); + + rsize = exp_in_limbs + 1; + rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + rp[0] = base; + rsize = 1; + + count_leading_zeros (cnt, exp_in_base); + for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) + { + mpn_mul_n (tp, rp, rp, rsize); + rsize = 2 * rsize; + rsize -= tp[rsize - 1] == 0; + xp = tp; tp = rp; rp = xp; + + if (((exp_in_base >> i) & 1) != 0) + { + cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base); + rp[rsize] = cy; + rsize += cy != 0; + } + } + + mp = u->_mp_d; + msize = usize; + + { + mp_ptr qp; + mp_limb_t qflag; + mp_size_t xtra; + if (msize < rsize) + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC ((rsize+1)* BYTES_PER_MP_LIMB); + MPN_ZERO (tmp, rsize - msize); + MPN_COPY (tmp + rsize - msize, mp, msize); + mp = tmp; + msize = rsize; + } + else + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC ((msize+1)* BYTES_PER_MP_LIMB); + MPN_COPY (tmp, mp, msize); + mp = tmp; + } + count_leading_zeros (cnt, rp[rsize - 1]); + cy = 0; + if (cnt != 0) + { + mpn_lshift (rp, rp, rsize, cnt); + cy = mpn_lshift (mp, mp, msize, cnt); + if (cy) + mp[msize++] = cy; + } + + { + mp_size_t qsize = n_limbs + (cy != 0); + qp = (mp_ptr) TMP_ALLOC ((qsize + 1) * BYTES_PER_MP_LIMB); + xtra = qsize - (msize - rsize); + qflag = mpn_divrem (qp, xtra, mp, msize, rp, rsize); + qp[qsize] = qflag; + rsize = qsize + qflag; + rp = qp; + } + } + } + + str_size = mpn_get_str (tstr, base, rp, rsize); + + if (str_size > n_digits + 3 * BITS_PER_MP_LIMB) + abort (); + + start_str = 0; + while (tstr[start_str] == 0) + start_str++; + + for (i = start_str; i < str_size; i++) + { + tstr[digits_computed_so_far++] = tstr[i]; + if (digits_computed_so_far > n_digits) + break; + } + exp_in_base = exp_in_base + str_size - start_str; + goto finish_up; + } + + exp_in_base = 0; + + if (uexp > 0) + { + /* The number has an integral part, convert that first. + If there is a fractional part too, it will be handled later. */ + mp_size_t start_str; + + rp = (mp_ptr) TMP_ALLOC (uexp * BYTES_PER_MP_LIMB); + up = u->_mp_d + usize - uexp; + MPN_COPY (rp, up, uexp); + + str_size = mpn_get_str (tstr, base, rp, uexp); + + start_str = 0; + while (tstr[start_str] == 0) + start_str++; + + for (i = start_str; i < str_size; i++) + { + tstr[digits_computed_so_far++] = tstr[i]; + if (digits_computed_so_far > n_digits) + { + exp_in_base = str_size - start_str; + goto finish_up; + } + } + + exp_in_base = str_size - start_str; + /* Modify somewhat and fall out to convert fraction... */ + usize -= uexp; + uexp = 0; + } + + if (usize <= 0) + goto finish_up; + + /* Convert the fraction. */ + { + mp_size_t rsize, msize; + mp_ptr rp, tp, xp, mp; + int cnt; + mp_limb_t cy; + mp_exp_t nexp; + + big_base = __mp_bases[base].big_base; + dig_per_u = __mp_bases[base].chars_per_limb; + + /* Hack for correctly (although not efficiently) converting to bases that + are powers of 2. If we deem it important, we could handle powers of 2 + by shifting and masking (just like mpn_get_str). */ + if (big_base < 10) /* logarithm of base when power of two */ + { + int logbase = big_base; + if (dig_per_u * logbase == BITS_PER_MP_LIMB) + dig_per_u--; + big_base = (mp_limb_t) 1 << (dig_per_u * logbase); + /* fall out to general code... */ + } + +#if 0 + if (0 && uexp == 0) + { + rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + up = u->_mp_d; + MPN_COPY (rp, up, usize); + rsize = usize; + nexp = 0; + } + else + {} +#endif + uexp = -uexp; + if (u->_mp_d[usize - 1] == 0) + cnt = 0; + else + count_leading_zeros (cnt, u->_mp_d[usize - 1]); + + nexp = ((uexp * BITS_PER_MP_LIMB) + cnt) + * __mp_bases[base].chars_per_bit_exactly; + + if (nexp == 0) + { + rp = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB); + up = u->_mp_d; + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + rsize = uexp + 2; + rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + rp[0] = base; + rsize = 1; + + count_leading_zeros (cnt, nexp); + for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) + { + mpn_mul_n (tp, rp, rp, rsize); + rsize = 2 * rsize; + rsize -= tp[rsize - 1] == 0; + xp = tp; tp = rp; rp = xp; + + if (((nexp >> i) & 1) != 0) + { + cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base); + rp[rsize] = cy; + rsize += cy != 0; + } + } + + /* Did our multiplier (base^nexp) cancel with uexp? */ +#if 0 + if (uexp != rsize) + { + do + { + cy = mpn_mul_1 (rp, rp, rsize, big_base); + nexp += dig_per_u; + } + while (cy == 0); + rp[rsize++] = cy; + } +#endif + mp = u->_mp_d; + msize = usize; + + tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB); + if (rsize > msize) + cy = mpn_mul (tp, rp, rsize, mp, msize); + else + cy = mpn_mul (tp, mp, msize, rp, rsize); + rsize += msize; + rsize -= cy == 0; + rp = tp; + + /* If we already output digits (for an integral part) pad + leading zeros. */ + if (digits_computed_so_far != 0) + for (i = 0; i < nexp; i++) + tstr[digits_computed_so_far++] = 0; + } + + while (digits_computed_so_far <= n_digits) + { + /* For speed: skip trailing zeroes. */ + if (rp[0] == 0) + { + rp++; + rsize--; + if (rsize == 0) + { + n_digits = digits_computed_so_far; + break; + } + } + + cy = mpn_mul_1 (rp, rp, rsize, big_base); + if (digits_computed_so_far == 0 && cy == 0) + { + abort (); + nexp += dig_per_u; + continue; + } + /* Convert N1 from BIG_BASE to a string of digits in BASE + using single precision operations. */ + { + unsigned char *s = tstr + digits_computed_so_far + dig_per_u; + for (i = dig_per_u - 1; i >= 0; i--) + { + *--s = cy % base; + cy /= base; + } + } + digits_computed_so_far += dig_per_u; + } + if (exp_in_base == 0) + exp_in_base = -nexp; + } + + finish_up: + + /* We can have at most one leading 0. Remove it. */ + if (tstr[0] == 0) + { + tstr++; + digits_computed_so_far--; + exp_in_base--; + } + + /* We should normally have computed too many digits. Round the result + at the point indicated by n_digits. */ + if (digits_computed_so_far > n_digits) + { + /* Round the result. */ + if (tstr[n_digits] * 2 >= base) + { + digits_computed_so_far = n_digits; + for (i = n_digits - 1; i >= 0; i--) + { + unsigned int x; + x = ++(tstr[i]); + if (x < base) + goto rounded_ok; + digits_computed_so_far--; + } + tstr[0] = 1; + digits_computed_so_far = 1; + exp_in_base++; + rounded_ok:; + } + } + + /* 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 (i.e., 0.125 in base 10). */ + if (n_digits > digits_computed_so_far) + n_digits = digits_computed_so_far; + + /* Remove trailing 0. There can be many zeros. */ + while (n_digits != 0 && tstr[n_digits - 1] == 0) + n_digits--; + + /* Translate to ascii and null-terminate. */ + for (i = 0; i < n_digits; i++) + *str++ = num_to_text[tstr[i]]; + *str = 0; + *exp = exp_in_base; + TMP_FREE (marker); + return digit_ptr; +} + +#if COPY_THIS_TO_OTHER_PLACES + /* Use this expression in lots of places in the library instead of the + count_leading_zeros+expression that is used currently. This expression + is much more accurate and will save odles of memory. */ + rsize = ((mp_size_t) (exp_in_base / __mp_bases[base].chars_per_bit_exactly) + + BITS_PER_MP_LIMB) / BITS_PER_MP_LIMB; +#endif diff --git a/gnu/lib/libgmp/mpf/init.c b/gnu/lib/libgmp/mpf/init.c new file mode 100644 index 00000000000..53701218f1c --- /dev/null +++ b/gnu/lib/libgmp/mpf/init.c @@ -0,0 +1,38 @@ +/* mpf_init() -- Make a new multiple precision number with value 0. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init (mpf_ptr r) +#else +mpf_init (r) + mpf_ptr r; +#endif +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + r->_mp_size = 0; + r->_mp_exp = 0; +} diff --git a/gnu/lib/libgmp/mpf/init2.c b/gnu/lib/libgmp/mpf/init2.c new file mode 100644 index 00000000000..a3e5752b24f --- /dev/null +++ b/gnu/lib/libgmp/mpf/init2.c @@ -0,0 +1,41 @@ +/* mpf_init2() -- Make a new multiple precision number with value 0. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init2 (mpf_ptr r, unsigned long int prec_in_bits) +#else +mpf_init2 (r, prec_in_bits) + mpf_ptr r; + unsigned long int prec_in_bits; +#endif +{ + mp_size_t prec; + + prec = (MAX (53, prec_in_bits) + 2 * BITS_PER_MP_LIMB - 1)/BITS_PER_MP_LIMB; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + r->_mp_size = 0; + r->_mp_exp = 0; +} diff --git a/gnu/lib/libgmp/mpf/inp_str.c b/gnu/lib/libgmp/mpf/inp_str.c new file mode 100644 index 00000000000..1d3cd4cc2b5 --- /dev/null +++ b/gnu/lib/libgmp/mpf/inp_str.c @@ -0,0 +1,89 @@ +/* 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 (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include <ctype.h> +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpf_inp_str (mpf_ptr rop, FILE *stream, int base) +#else +mpf_inp_str (rop, stream, base) + mpf_ptr rop; + FILE *stream; + int base; +#endif +{ + char *str; + size_t alloc_size, str_size; + int c; + size_t retval; + size_t nread; + + if (stream == 0) + stream = stdin; + + alloc_size = 100; + str = (char *) (*_mp_allocate_func) (alloc_size); + str_size = 0; + nread = 0; + + /* Skip whitespace. */ + do + { + c = getc (stream); + nread++; + } + while (isspace (c)); + + for (;;) + { + if (str_size >= alloc_size) + { + size_t old_alloc_size = alloc_size; + alloc_size = alloc_size * 3 / 2; + str = (char *) (*_mp_reallocate_func) (str, old_alloc_size, alloc_size); + } + if (c == EOF || isspace (c)) + break; + str[str_size++] = c; + c = getc (stream); + } + ungetc (c, stream); + + if (str_size >= alloc_size) + { + size_t old_alloc_size = alloc_size; + alloc_size = alloc_size * 3 / 2; + str = (char *) (*_mp_reallocate_func) (str, old_alloc_size, alloc_size); + } + str[str_size] = 0; + + retval = mpf_set_str (rop, str, base); + if (retval == -1) + return 0; /* error */ + + (*_mp_free_func) (str, alloc_size); + return str_size + nread; +} diff --git a/gnu/lib/libgmp/mpf/iset.c b/gnu/lib/libgmp/mpf/iset.c new file mode 100644 index 00000000000..c2362e1729a --- /dev/null +++ b/gnu/lib/libgmp/mpf/iset.c @@ -0,0 +1,59 @@ +/* mpf_init_set -- Initialize a float and assign it from another float. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init_set (mpf_ptr r, mpf_srcptr s) +#else +mpf_init_set (r, s) + mpf_ptr r; + mpf_srcptr s; +#endif +{ + mp_ptr rp, sp; + mp_size_t ssize, size; + mp_size_t prec; + + prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + 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; + } + + MPN_COPY (rp, sp, size); + + r->_mp_exp = s->_mp_exp; + r->_mp_size = ssize >= 0 ? size : -size; +} diff --git a/gnu/lib/libgmp/mpf/iset_d.c b/gnu/lib/libgmp/mpf/iset_d.c new file mode 100644 index 00000000000..d09e8b89a5f --- /dev/null +++ b/gnu/lib/libgmp/mpf/iset_d.c @@ -0,0 +1,39 @@ +/* mpf_init_set_d -- Initialize a float and assign it from a double. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init_set_d (mpf_ptr r, double val) +#else +mpf_init_set_d (r, val) + mpf_ptr r; + double val; +#endif +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + + mpf_set_d (r, val); +} diff --git a/gnu/lib/libgmp/mpf/iset_si.c b/gnu/lib/libgmp/mpf/iset_si.c new file mode 100644 index 00000000000..e67eef4cdb5 --- /dev/null +++ b/gnu/lib/libgmp/mpf/iset_si.c @@ -0,0 +1,55 @@ +/* mpf_init_set_si() -- Initialize a float and assign it from a signed int. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init_set_si (mpf_ptr r, long int val) +#else +mpf_init_set_si (r, val) + mpf_ptr r; + long int val; +#endif +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + + if (val > 0) + { + r->_mp_d[0] = val; + r->_mp_size = 1; + r->_mp_exp = 1; + } + else if (val < 0) + { + r->_mp_d[0] = -val; + r->_mp_size = -1; + r->_mp_exp = 1; + } + else + { + r->_mp_size = 0; + r->_mp_exp = 0; + } +} diff --git a/gnu/lib/libgmp/mpf/iset_str.c b/gnu/lib/libgmp/mpf/iset_str.c new file mode 100644 index 00000000000..40bfe7cf27e --- /dev/null +++ b/gnu/lib/libgmp/mpf/iset_str.c @@ -0,0 +1,40 @@ +/* mpf_init_set_str -- Initialize a float and assign it from a string. + +Copyright (C) 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +int +#if __STDC__ +mpf_init_set_str (mpf_ptr r, const char *s, int base) +#else +mpf_init_set_str (r, s, base) + mpf_ptr r; + const char *s; + int base; +#endif +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + + return mpf_set_str (r, s, base); +} diff --git a/gnu/lib/libgmp/mpf/iset_ui.c b/gnu/lib/libgmp/mpf/iset_ui.c new file mode 100644 index 00000000000..c6c24afaf1a --- /dev/null +++ b/gnu/lib/libgmp/mpf/iset_ui.c @@ -0,0 +1,40 @@ +/* mpf_init_set_ui() -- Initialize a float and assign it from an unsigned int. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_init_set_ui (mpf_ptr r, unsigned long int val) +#else +mpf_init_set_ui (r, val) + mpf_ptr r; + unsigned long int val; +#endif +{ + mp_size_t prec = __gmp_default_fp_limb_precision; + r->_mp_d = (mp_ptr) (*_mp_allocate_func) ((prec + 1) * BYTES_PER_MP_LIMB); + r->_mp_prec = prec; + r->_mp_d[0] = val; + r->_mp_size = val != 0; + r->_mp_exp = val != 0; +} diff --git a/gnu/lib/libgmp/mpf/mul.c b/gnu/lib/libgmp/mpf/mul.c new file mode 100644 index 00000000000..8f1a2a9aca1 --- /dev/null +++ b/gnu/lib/libgmp/mpf/mul.c @@ -0,0 +1,94 @@ +/* mpf_mul -- Multiply two floats. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_mul (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +#else +mpf_mul (r, u, v) + mpf_ptr r; + mpf_srcptr u; + mpf_srcptr v; +#endif +{ + mp_srcptr up, vp; + mp_size_t usize, vsize; + mp_size_t sign_product; + mp_size_t prec = r->_mp_prec; + TMP_DECL (marker); + + TMP_MARK (marker); + usize = u->_mp_size; + vsize = v->_mp_size; + sign_product = usize ^ vsize; + + usize = ABS (usize); + vsize = ABS (vsize); + + up = u->_mp_d; + vp = v->_mp_d; + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + + if (usize == 0 || vsize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; /* ??? */ + } + else + { + mp_size_t rsize; + mp_limb_t cy_limb; + mp_ptr rp, tp; + mp_size_t adj; + + rsize = usize + vsize; + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + 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 (marker); +} diff --git a/gnu/lib/libgmp/mpf/mul_2exp.c b/gnu/lib/libgmp/mpf/mul_2exp.c new file mode 100644 index 00000000000..0ed35a0b1c1 --- /dev/null +++ b/gnu/lib/libgmp/mpf/mul_2exp.c @@ -0,0 +1,89 @@ +/* mpf_mul_2exp -- Multiply a float by 2^n. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_mul_2exp (mpf_ptr r, mpf_srcptr u, unsigned long int exp) +#else +mpf_mul_2exp (r, u, exp) + mpf_ptr r; + mpf_srcptr u; + unsigned long int exp; +#endif +{ + 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 (usize == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + abs_usize = ABS (usize); + up = u->_mp_d; + + if (abs_usize > prec) + { + up += abs_usize - prec; + abs_usize = prec; + } + + if (exp % BITS_PER_MP_LIMB == 0) + { + if (rp != up) + MPN_COPY (rp, up, abs_usize); + r->_mp_size = usize >= 0 ? abs_usize : -abs_usize; + r->_mp_exp = uexp + exp / BITS_PER_MP_LIMB; + } + else + { + mp_limb_t cy_limb; + if (r != u) + { + cy_limb = mpn_lshift (rp, up, abs_usize, exp % BITS_PER_MP_LIMB); + rp[abs_usize] = cy_limb; + cy_limb = cy_limb != 0; + } + else + { + /* Use mpn_rshift since mpn_lshift operates downwards, and we + therefore would clobber part of U before using that part. */ + cy_limb = mpn_rshift (rp + 1, up, abs_usize, -exp % BITS_PER_MP_LIMB); + rp[0] = cy_limb; + cy_limb = rp[abs_usize] != 0; + } + + abs_usize += cy_limb; + r->_mp_size = usize >= 0 ? abs_usize : -abs_usize; + r->_mp_exp = uexp + exp / BITS_PER_MP_LIMB + cy_limb; + } +} diff --git a/gnu/lib/libgmp/mpf/mul_ui.c b/gnu/lib/libgmp/mpf/mul_ui.c new file mode 100644 index 00000000000..fc87cb50750 --- /dev/null +++ b/gnu/lib/libgmp/mpf/mul_ui.c @@ -0,0 +1,74 @@ +/* mpf_mul_ui -- Multiply a float and an unsigned integer. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_mul_ui (mpf_ptr r, mpf_srcptr u, unsigned long int v) +#else +mpf_mul_ui (r, u, v) + mpf_ptr r; + mpf_srcptr u; + unsigned long int v; +#endif +{ + mp_srcptr up; + mp_size_t usize; + mp_size_t size; + mp_size_t prec = r->_mp_prec; + mp_limb_t cy_limb; + mp_ptr rp; + + usize = u->_mp_size; + size = ABS (usize); + + rp = r->_mp_d; + up = u->_mp_d; + if (size > prec) + { + up += size - prec; + size = prec; + } + +#if 0 + /* Since we can do it at almost no cost, remove zero limbs at low end of + result. */ + if (up[0] == 0) + up++, size--; +#endif + + if (size == 0 || v == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; /* ??? */ + } + else + { + cy_limb = mpn_mul_1 (rp, up, size, (mp_limb_t) v); + 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/gnu/lib/libgmp/mpf/neg.c b/gnu/lib/libgmp/mpf/neg.c new file mode 100644 index 00000000000..a4139fb9054 --- /dev/null +++ b/gnu/lib/libgmp/mpf/neg.c @@ -0,0 +1,59 @@ +/* mpf_neg -- Negate a float. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_neg (mpf_ptr r, mpf_srcptr u) +#else +mpf_neg (r, u) + mpf_ptr r; + mpf_srcptr u; +#endif +{ + 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/gnu/lib/libgmp/mpf/out_str.c b/gnu/lib/libgmp/mpf/out_str.c new file mode 100644 index 00000000000..2a7a954efde --- /dev/null +++ b/gnu/lib/libgmp/mpf/out_str.c @@ -0,0 +1,89 @@ +/* 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 (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpf_out_str (FILE *stream, int base, size_t n_digits, mpf_srcptr op) +#else +mpf_out_str (stream, base, n_digits, op) + FILE *stream; + int base; + size_t n_digits; + mpf_srcptr op; +#endif +{ + char *str; + mp_exp_t exp; + size_t written; + TMP_DECL (marker); + + TMP_MARK (marker); + + if (base == 0) + base = 10; + if (n_digits == 0) + n_digits = (((op->_mp_prec - 1) * BITS_PER_MP_LIMB) + * __mp_bases[base].chars_per_bit_exactly); + + if (stream == 0) + stream = stdout; + + 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; + } + + fwrite ("0.", 1, 2, stream); + written += 2; + + /* 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 (marker); + return ferror (stream) ? 0 : written; +} diff --git a/gnu/lib/libgmp/mpf/random2.c b/gnu/lib/libgmp/mpf/random2.c new file mode 100644 index 00000000000..8e67b11facb --- /dev/null +++ b/gnu/lib/libgmp/mpf/random2.c @@ -0,0 +1,65 @@ +/* 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 (C) 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +#if defined (__hpux) || defined (alpha__) || defined (__svr4__) || defined (__SVR4) +/* HPUX lacks random(). DEC OSF/1 1.2 random() returns a double. */ +long mrand48 (); +static inline long +random () +{ + return mrand48 (); +} +#else +long random (); +#endif + +void +#if __STDC__ +mpf_random2 (mpf_ptr x, mp_size_t size, mp_exp_t exp) +#else +mpf_random2 (x, size, exp) + mpf_ptr x; + mp_size_t size; + mp_exp_t exp; +#endif +{ + mp_size_t asize; + mp_size_t prec = x->_mp_prec; + + asize = ABS (size); + if (asize != 0) + { + if (asize > prec + 1) + asize = prec + 1; + + mpn_random2 (x->_mp_d, asize); + } + + if (exp != 0) + exp = random () % (2 * exp) - exp; + x->_mp_exp = asize == 0 ? 0 : exp; + x->_mp_size = size < 0 ? -asize : asize; +} diff --git a/gnu/lib/libgmp/mpf/reldiff.c b/gnu/lib/libgmp/mpf/reldiff.c new file mode 100644 index 00000000000..c8464660753 --- /dev/null +++ b/gnu/lib/libgmp/mpf/reldiff.c @@ -0,0 +1,52 @@ +/* mpf_reldiff -- Generate the relative difference of two floats. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_reldiff (mpf_t rdiff, mpf_srcptr x, mpf_srcptr y) +#else +mpf_reldiff (rdiff, x, y) + mpf_t rdiff; + mpf_srcptr x; + mpf_srcptr y; +#endif +{ + if (mpf_cmp_ui (x, 0) == 0) + { + mpf_set_ui (rdiff, (unsigned long int) (mpf_sgn (y) != 0)); + } + else + { + mpf_t d; + mp_limb_t tmp_limb[2]; + + d->_mp_prec = 1; + d->_mp_d = tmp_limb; + + mpf_sub (d, x, y); + mpf_abs (d, d); + mpf_div (rdiff, d, x); + } +} + diff --git a/gnu/lib/libgmp/mpf/set.c b/gnu/lib/libgmp/mpf/set.c new file mode 100644 index 00000000000..5778b85b093 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set.c @@ -0,0 +1,53 @@ +/* mpf_set -- Assign a float from another float. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set (mpf_ptr r, mpf_srcptr u) +#else +mpf_set (r, u) + mpf_ptr r; + mpf_srcptr u; +#endif +{ + 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; + } + + MPN_COPY (rp, up, asize); + r->_mp_exp = u->_mp_exp; + r->_mp_size = size >= 0 ? asize : -asize; +} diff --git a/gnu/lib/libgmp/mpf/set_d.c b/gnu/lib/libgmp/mpf/set_d.c new file mode 100644 index 00000000000..43eb1577457 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_d.c @@ -0,0 +1,47 @@ +/* mpf_set_d -- Assign a float from a IEEE double. + +Copyright (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_d (mpf_ptr r, double d) +#else +mpf_set_d (r, d) + mpf_ptr r; + double d; +#endif +{ + int negative; + + if (d == 0) + { + SIZ(r) = 0; + EXP(r) = 0; + return; + } + negative = d < 0; + d = ABS (d); + + EXP(r) = __gmp_extract_double (PTR(r), d); + SIZ(r) = negative ? -LIMBS_PER_DOUBLE : LIMBS_PER_DOUBLE; +} diff --git a/gnu/lib/libgmp/mpf/set_dfl_prec.c b/gnu/lib/libgmp/mpf/set_dfl_prec.c new file mode 100644 index 00000000000..55069e4e4f1 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_dfl_prec.c @@ -0,0 +1,40 @@ +/* mpf_set_default_prec -- + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +mp_size_t __gmp_default_fp_limb_precision + = (53 + 2 * BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB; + +void +#if __STDC__ +mpf_set_default_prec (unsigned long int prec_in_bits) +#else +mpf_set_default_prec (prec_in_bits) + unsigned long int prec_in_bits; +#endif +{ + mp_size_t prec; + + prec = (MAX (53, prec_in_bits) + 2 * BITS_PER_MP_LIMB - 1)/BITS_PER_MP_LIMB; + __gmp_default_fp_limb_precision = prec; +} diff --git a/gnu/lib/libgmp/mpf/set_prc.c b/gnu/lib/libgmp/mpf/set_prc.c new file mode 100644 index 00000000000..10f2b0671f5 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_prc.c @@ -0,0 +1,57 @@ +/* mpf_set_prec(x) -- Change the precision of x. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_prec (mpf_ptr x, unsigned long int prec_in_bits) +#else +mpf_set_prec (x, prec_in_bits) + mpf_ptr x; + unsigned long int prec_in_bits; +#endif +{ + mp_size_t prec; + mp_size_t size = ABS (x->_mp_size); + + prec = (MAX (53, prec_in_bits) + 2 * BITS_PER_MP_LIMB - 1)/BITS_PER_MP_LIMB; + + /* We want the most significant limbs, so move the limbs down if we are + about to truncate the value. */ + if (size > prec + 1) + { + mp_size_t offset = size - (prec + 1); + mp_ptr xp = x->_mp_d; + + MPN_COPY (xp, xp + offset, prec + 1); + } + + x->_mp_d = (mp_ptr) (*_mp_reallocate_func) + (x->_mp_d, + (x->_mp_prec + 1) * BYTES_PER_MP_LIMB, (prec + 1) * BYTES_PER_MP_LIMB); + x->_mp_prec = prec; + + /* If the precision decreased, truncate the number. */ + if (size > prec + 1) + x->_mp_size = x->_mp_size >= 0 ? (prec + 1) : -(prec + 1); +} diff --git a/gnu/lib/libgmp/mpf/set_prc_raw.c b/gnu/lib/libgmp/mpf/set_prc_raw.c new file mode 100644 index 00000000000..f55188a1248 --- /dev/null +++ b/gnu/lib/libgmp/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 (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_prec_raw (mpf_ptr x, unsigned long int prec_in_bits) +#else +mpf_set_prec_raw (x, prec_in_bits) + mpf_ptr x; + unsigned long int prec_in_bits; +#endif +{ + mp_size_t prec; + prec = (MAX (53, prec_in_bits) + 2 * BITS_PER_MP_LIMB - 1)/BITS_PER_MP_LIMB; + x->_mp_prec = prec; +} diff --git a/gnu/lib/libgmp/mpf/set_q.c b/gnu/lib/libgmp/mpf/set_q.c new file mode 100644 index 00000000000..69d81e33e43 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_q.c @@ -0,0 +1,170 @@ +/* mpf_set_q (mpf_t rop, mpq_t op) -- Convert the rational op to the float rop. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +/* Algorithm: + 1. Develop >= n bits of src.num / src.den, where n is the number of bits + in a double. This (partial) division will use all bits from the + denominator. + 2. Use the remainder to determine how to round the result. + 3. Assign the integral result to a temporary double. + 4. Scale the temporary double, and return the result. + + An alternative algorithm, that would be faster: + 0. Let n be somewhat larger than the number of significant bits in a double. + 1. Extract the most significant n bits of the denominator, and an equal + number of bits from the numerator. + 2. Interpret the extracted numbers as integers, call them a and b + respectively, and develop n bits of the fractions ((a + 1) / b) and + (a / (b + 1)) using mpn_divrem. + 3. If the computed values are identical UP TO THE POSITION WE CARE ABOUT, + we are done. If they are different, repeat the algorithm from step 1, + but first let n = n * 2. + 4. If we end up using all bits from the numerator and denominator, fall + back to the first algorithm above. + 5. Just to make life harder, The computation of a + 1 and b + 1 above + might give carry-out... Needs special handling. It might work to + subtract 1 in both cases instead. +*/ + +void +#if __STDC__ +mpf_set_q (mpf_t r, mpq_srcptr q) +#else +mpf_set_q (r, q) + mpf_t r; + mpq_srcptr q; +#endif +{ + mp_ptr np, dp; + mp_ptr rp; + mp_size_t nsize, dsize; + mp_size_t qsize, rsize; + mp_size_t sign_quotient; + unsigned normalization_steps; + mp_limb_t qlimb; + mp_ptr qp; + mp_size_t prec; + mp_exp_t exp; + TMP_DECL (marker); + + nsize = SIZ (&q->_mp_num); + dsize = SIZ (&q->_mp_den); + + if (nsize == 0) + { + SIZ (r) = 0; + EXP (r) = 0; + return; + } + + prec = PREC (r) + 1; + + TMP_MARK (marker); + + qp = PTR (r); + + sign_quotient = nsize ^ dsize; + nsize = ABS (nsize); + dsize = ABS (dsize); + np = PTR (&q->_mp_num); + dp = PTR (&q->_mp_den); + + exp = nsize - dsize; + + if (nsize > prec) + { + np += nsize - prec; + nsize = prec; + } + if (dsize > prec) + { + dp += dsize - prec; + dsize = prec; + } + + rsize = MAX (nsize, dsize); + rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB); + + count_leading_zeros (normalization_steps, dp[dsize - 1]); + + /* Normalize the denominator, i.e. make its most significant bit set by + shifting it NORMALIZATION_STEPS bits to the left. Also shift the + numerator the same number of steps (to keep the quotient the same!). */ + if (normalization_steps != 0) + { + mp_ptr tp; + mp_limb_t nlimb; + + /* Shift up the denominator setting the most significant bit of + the most significant limb. Use temporary storage not to clobber + the original contents of the denominator. */ + tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB); + mpn_lshift (tp, dp, dsize, normalization_steps); + dp = tp; + + if (rsize != nsize) + { + MPN_ZERO (rp, rsize - nsize); + nlimb = mpn_lshift (rp + (rsize - nsize), + np, nsize, normalization_steps); + } + else + { + nlimb = mpn_lshift (rp, np, nsize, normalization_steps); + } + if (nlimb != 0) + { + rp[rsize] = nlimb; + rsize++; + exp++; + } + } + else + { + if (rsize != nsize) + { + MPN_ZERO (rp, rsize - nsize); + MPN_COPY (rp + (rsize - nsize), np, nsize); + } + else + { + MPN_COPY (rp, np, rsize); + } + } + + qlimb = mpn_divrem (qp, prec - 1 - (rsize - dsize), rp, rsize, dp, dsize); + qsize = prec - 1; + if (qlimb) + { + qp[qsize] = qlimb; + qsize++; + exp++; + } + + EXP (r) = exp; + SIZ (r) = qsize; + + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/set_si.c b/gnu/lib/libgmp/mpf/set_si.c new file mode 100644 index 00000000000..f9b4b55eb25 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_si.c @@ -0,0 +1,51 @@ +/* mpf_set_si() -- Assign a float from a signed int. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_si (mpf_ptr x, long int val) +#else +mpf_set_si (x, val) + mpf_ptr x; + long int val; +#endif +{ + if (val > 0) + { + x->_mp_d[0] = val; + x->_mp_size = 1; + x->_mp_exp = 1; + } + else if (val < 0) + { + x->_mp_d[0] = -val; + x->_mp_size = -1; + x->_mp_exp = 1; + } + else + { + x->_mp_size = 0; + x->_mp_exp = 0; + } +} diff --git a/gnu/lib/libgmp/mpf/set_str.c b/gnu/lib/libgmp/mpf/set_str.c new file mode 100644 index 00000000000..2ab9faebe48 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_str.c @@ -0,0 +1,302 @@ +/* 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 (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <string.h> +#include <ctype.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +long int strtol _PROTO ((const char *, char **ptr, int)); + +static int +digit_value_in_base (c, base) + int c; + int base; +{ + int digit; + + if (isdigit (c)) + digit = c - '0'; + else if (islower (c)) + digit = c - 'a' + 10; + else if (isupper (c)) + digit = c - 'A' + 10; + else + return -1; + + if (digit < base) + return digit; + return -1; +} + +int +#if __STDC__ +mpf_set_str (mpf_ptr x, const char *str, int base) +#else +mpf_set_str (x, str, base) + mpf_ptr x; + char *str; + int base; +#endif +{ + size_t str_size; + char *s, *begs; + size_t i; + mp_size_t xsize; + int c; + int negative; + char *dotpos = 0; + int expflag; + int decimal_exponent_flag; + TMP_DECL (marker); + + TMP_MARK (marker); + + c = *str; + + /* Skip whitespace. */ + while (isspace (c)) + c = *++str; + + negative = 0; + if (c == '-') + { + negative = 1; + c = *++str; + } + + decimal_exponent_flag = base < 0; + base = ABS (base); + + if (digit_value_in_base (c, base == 0 ? 10 : base) < 0) + return -1; /* error if no digits */ + + /* If BASE is 0, try to find out the base by looking at the initial + characters. */ + if (base == 0) + { + base = 10; +#if 0 + if (c == '0') + { + base = 8; + c = *++str; + if (c == 'x' || c == 'X') + base = 16; + } +#endif + } + + expflag = 0; + str_size = strlen (str); + for (i = 0; i < str_size; i++) + { + c = str[i]; + if (c == '@' || (base <= 10 && (c == 'e' || c == 'E'))) + { + expflag = 1; + str_size = i; + break; + } + + } + + s = begs = (char *) TMP_ALLOC (str_size + 1); + + for (i = 0; i < str_size; i++) + { + c = *str; + if (!isspace (c)) + { + int dig; + + if (c == '.') + { + if (dotpos != 0) + { + TMP_FREE (marker); + return -1; + } + dotpos = s; + } + else + { + dig = digit_value_in_base (c, base); + if (dig < 0) + { + TMP_FREE (marker); + return -1; + } + *s++ = dig; + } + } + c = *++str; + } + + str_size = s - begs; + + xsize = str_size / __mp_bases[base].chars_per_limb + 2; + { + long exp_in_base; + mp_size_t rsize, msize; + int cnt, i; + mp_ptr mp, xp, tp, rp; + mp_limb_t cy; + mp_exp_t exp_in_limbs; + mp_size_t prec = x->_mp_prec; + int divflag; + mp_size_t xxx = 0; + + mp = (mp_ptr) TMP_ALLOC (xsize * BYTES_PER_MP_LIMB); + msize = mpn_set_str (mp, (unsigned char *) begs, str_size, base); + + if (msize == 0) + { + x->_mp_size = 0; + x->_mp_exp = 0; + TMP_FREE (marker); + return 0; + } + + if (expflag != 0) + exp_in_base = strtol (str + 1, (char **) 0, + decimal_exponent_flag ? 10 : base); + else + exp_in_base = 0; + if (dotpos != 0) + exp_in_base -= s - dotpos; + divflag = exp_in_base < 0; + exp_in_base = ABS (exp_in_base); + + if (exp_in_base == 0) + { + MPN_COPY (x->_mp_d, mp, msize); + x->_mp_size = negative ? -msize : msize; + x->_mp_exp = msize; + TMP_FREE (marker); + return 0; + } + +#if 1 + rsize = (((mp_size_t) (exp_in_base / __mp_bases[base].chars_per_bit_exactly)) + / BITS_PER_MP_LIMB + 3); +#else + count_leading_zeros (cnt, (mp_limb_t) base); + rsize = exp_in_base - cnt * exp_in_base / BITS_PER_MP_LIMB + 1; +#endif + rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + rp[0] = base; + rsize = 1; + + count_leading_zeros (cnt, exp_in_base); + + for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--) + { + mpn_mul_n (tp, rp, rp, rsize); + rsize = 2 * rsize; + rsize -= tp[rsize - 1] == 0; + xp = tp; tp = rp; rp = xp; + + if (((exp_in_base >> i) & 1) != 0) + { + cy = mpn_mul_1 (rp, rp, rsize, (mp_limb_t) base); + rp[rsize] = cy; + rsize += cy != 0; + } + } + + if (rsize > prec) + { + xxx += rsize - prec; + rp += rsize - prec; + rsize = prec; + } +#if 0 + if (msize > prec) + { + xxx -= msize - prec; + mp += msize - prec; + msize = prec; + } +#endif + if (divflag) + { + mp_ptr qp; + mp_limb_t qflag; + mp_size_t xtra; + if (msize <= rsize) + { + /* Allocate extra limb for current divrem sematics. */ + mp_ptr tmp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (tmp, rsize - msize); + MPN_COPY (tmp + rsize - msize, mp, msize); + mp = tmp; + xxx += rsize - msize; + msize = rsize; + } + count_leading_zeros (cnt, rp[rsize - 1]); + if (cnt != 0) + { + mpn_lshift (rp, rp, rsize, cnt); + cy = mpn_lshift (mp, mp, msize, cnt); + if (cy) + mp[msize++] = cy; + } + qp = (mp_ptr) TMP_ALLOC ((prec + 1) * BYTES_PER_MP_LIMB); + xtra = prec - (msize - rsize); + qflag = mpn_divrem (qp, xtra, mp, msize, rp, rsize); + qp[prec] = qflag; + tp = qp; + rsize = prec + qflag; + exp_in_limbs = rsize - xtra - xxx; + } + else + { + tp = (mp_ptr) TMP_ALLOC ((rsize + msize) * BYTES_PER_MP_LIMB); + if (rsize > msize) + mpn_mul (tp, rp, rsize, mp, msize); + else + mpn_mul (tp, mp, msize, rp, rsize); + rsize += msize; + rsize -= tp[rsize - 1] == 0; + exp_in_limbs = rsize + xxx; + + if (rsize > prec) + { + xxx = rsize - prec; + tp += rsize - prec; + rsize = prec; + exp_in_limbs += 0; + } + } + + MPN_COPY (x->_mp_d, tp, rsize); + x->_mp_size = negative ? -rsize : rsize; + x->_mp_exp = exp_in_limbs; + TMP_FREE (marker); + return 0; + } +} diff --git a/gnu/lib/libgmp/mpf/set_ui.c b/gnu/lib/libgmp/mpf/set_ui.c new file mode 100644 index 00000000000..ead0498aa73 --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_ui.c @@ -0,0 +1,45 @@ +/* mpf_set_ui() -- Assign a float from an unsigned int. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_ui (mpf_ptr x, unsigned long int val) +#else +mpf_set_ui (x, val) + mpf_ptr x; + unsigned long int val; +#endif +{ + if (val != 0) + { + x->_mp_d[0] = val; + x->_mp_size = 1; + x->_mp_exp = 1; + } + else + { + x->_mp_size = 0; + x->_mp_exp = 0; + } +} diff --git a/gnu/lib/libgmp/mpf/set_z.c b/gnu/lib/libgmp/mpf/set_z.c new file mode 100644 index 00000000000..4fff86edc3b --- /dev/null +++ b/gnu/lib/libgmp/mpf/set_z.c @@ -0,0 +1,54 @@ +/* mpf_set_z -- Assign a float from an integer. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_set_z (mpf_ptr r, mpz_srcptr u) +#else +mpf_set_z (r, u) + mpf_ptr r; + mpz_srcptr u; +#endif +{ + 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; + } + + MPN_COPY (rp, up, asize); + SIZ (r) = size >= 0 ? asize : -asize; +} diff --git a/gnu/lib/libgmp/mpf/size.c b/gnu/lib/libgmp/mpf/size.c new file mode 100644 index 00000000000..23a57ec72c9 --- /dev/null +++ b/gnu/lib/libgmp/mpf/size.c @@ -0,0 +1,35 @@ +/* mpf_size(x) -- return the number of limbs currently used by the + value of the float X. + +Copyright (C) 1993, 1994, 1995 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +size_t +#if __STDC__ +mpf_size (mpf_srcptr x) +#else +mpf_size (x) + mpf_srcptr x; +#endif +{ + return ABS (x->_mp_size); +} diff --git a/gnu/lib/libgmp/mpf/sqrt.c b/gnu/lib/libgmp/mpf/sqrt.c new file mode 100644 index 00000000000..6d115e41112 --- /dev/null +++ b/gnu/lib/libgmp/mpf/sqrt.c @@ -0,0 +1,75 @@ +/* mpf_sqrt -- Compute the square root of a float. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_sqrt (mpf_ptr r, mpf_srcptr u) +#else +mpf_sqrt (r, u) + mpf_ptr r; + mpf_srcptr u; +#endif +{ + mp_size_t usize; + mp_ptr up, tp; + mp_size_t prec; + mp_exp_t tsize, rexp; + TMP_DECL (marker); + + usize = u->_mp_size; + if (usize <= 0) + { + usize = 1 - 1 / (usize == 0); /* Divide by zero for negative OP. */ + r->_mp_size = usize; /* cheat flow by using usize here */ + r->_mp_exp = 0; + return; + } + + TMP_MARK (marker); + + prec = r->_mp_prec; + rexp = (u->_mp_exp + 1) >> 1; /* round towards -inf */ + tsize = 2 * prec + (u->_mp_exp & 1); + + up = u->_mp_d; + tp = (mp_ptr) TMP_ALLOC (tsize * BYTES_PER_MP_LIMB); + + 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); + + r->_mp_size = (tsize + 1) / 2; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/sqrt_ui.c b/gnu/lib/libgmp/mpf/sqrt_ui.c new file mode 100644 index 00000000000..421aa021993 --- /dev/null +++ b/gnu/lib/libgmp/mpf/sqrt_ui.c @@ -0,0 +1,61 @@ +/* mpf_sqrt_ui -- Compute the square root of an unsigned integer. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_sqrt_ui (mpf_ptr r, unsigned long int u) +#else +mpf_sqrt_ui (r, u) + mpf_ptr r; + unsigned long int u; +#endif +{ + mp_size_t rsize; + mp_ptr tp; + mp_size_t prec; + TMP_DECL (marker); + + if (u == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK (marker); + + prec = r->_mp_prec; + rsize = 2 * prec + 1; + + tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB); + + MPN_ZERO (tp, rsize - 1); + tp[rsize - 1] = u; + + mpn_sqrtrem (r->_mp_d, NULL, tp, rsize); + + r->_mp_size = prec + 1; + r->_mp_exp = 1; + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/sub.c b/gnu/lib/libgmp/mpf/sub.c new file mode 100644 index 00000000000..b8719826300 --- /dev/null +++ b/gnu/lib/libgmp/mpf/sub.c @@ -0,0 +1,402 @@ +/* mpf_sub -- Subtract two floats. + +Copyright (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_sub (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) +#else +mpf_sub (r, u, v) + mpf_ptr r; + mpf_srcptr u; + mpf_srcptr v; +#endif +{ + 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 (marker); + + usize = u->_mp_size; + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (usize == 0) + { + mpf_neg (r, v); + return; + } + if (vsize == 0) + { + mpf_set (r, u); + return; + } + + /* If signs of U and V are different, perform addition. */ + if ((usize ^ vsize) < 0) + { + __mpf_struct v_negated; + v_negated._mp_size = -vsize; + v_negated._mp_exp = v->_mp_exp; + v_negated._mp_d = v->_mp_d; + mpf_add (r, u, &v_negated); + return; + } + + TMP_MARK (marker); + + /* Signs are now known to be the same. */ + negate = usize < 0; + + /* Make U be the operand with the largest exponent. */ + if (u->_mp_exp < v->_mp_exp) + { + mpf_srcptr t; + t = u; u = v; v = t; + negate ^= 1; + usize = u->_mp_size; + vsize = v->_mp_size; + } + + usize = ABS (usize); + vsize = ABS (vsize); + up = u->_mp_d; + vp = v->_mp_d; + rp = r->_mp_d; + prec = r->_mp_prec + 1; + exp = u->_mp_exp; + ediff = u->_mp_exp - v->_mp_exp; + + /* If ediff is 0 or 1, we might have a situation where the operands are + extremely close. We need to scan the operands from the most significant + end ignore the initial parts that are equal. */ + if (ediff <= 1) + { + if (ediff == 0) + { + /* Skip leading limbs in U and V that are equal. */ + if (up[usize - 1] == vp[vsize - 1]) + { + /* This loop normally exits immediately. Optimize for that. */ + do + { + usize--; + vsize--; + exp--; + + if (usize == 0) + { + rsize = vsize; + tp = (mp_ptr) vp; + negate ^= 1; + goto normalize; + } + if (vsize == 0) + { + rsize = usize; + tp = (mp_ptr) up; + goto normalize; + } + } + while (up[usize - 1] == vp[vsize - 1]); + } + + if (up[usize - 1] < vp[vsize - 1]) + { + /* For simplicity, swap U and V. Note that since the loop above + wouldn't have exited unless up[usize - 1] and vp[vsize - 1] + were non-equal, this if-statement catches all cases where U + is smaller than V. */ + { mp_srcptr tp = up; up = vp; vp = tp; } + { mp_size_t tsize = usize; usize = vsize; vsize = tsize; } + 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] != ~(mp_limb_t) 0 + || (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] == ~(mp_limb_t) 0) + { + usize--; + vsize--; + exp--; + } + + if (usize == 0) + { + while (vsize != 0 && vp[vsize - 1] == ~(mp_limb_t) 0) + { + vsize--; + exp--; + } + } + + if (usize > prec - 1) + { + up += usize - (prec - 1); + usize = prec - 1; + } + if (vsize > prec - 1) + { + vp += vsize - (prec - 1); + vsize = prec - 1; + } + + tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB); + { + mp_limb_t cy_limb; + if (vsize == 0) + { + mp_size_t size, i; + size = usize; + for (i = 0; i < size; i++) + tp[i] = up[i]; + tp[size] = 1; + rsize = size + 1; + exp++; + goto normalize; + } + if (usize == 0) + { + mp_size_t size, i; + size = vsize; + for (i = 0; i < size; i++) + tp[i] = ~vp[i]; + cy_limb = 1 - mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1); + rsize = vsize; + if (cy_limb == 0) + { + tp[rsize] = 1; + rsize++; + exp++; + } + goto normalize; + } + if (usize >= vsize) + { + /* uuuu */ + /* vv */ + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + cy_limb = mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else /* (usize < vsize) */ + { + /* uuuu */ + /* vvvvvvv */ + mp_size_t size, i; + size = vsize - usize; + for (i = 0; i < size; i++) + tp[i] = ~vp[i]; + cy_limb = mpn_sub_n (tp + size, up, vp + size, usize); + cy_limb+= mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + cy_limb-= mpn_add_1 (tp, tp, vsize, (mp_limb_t) 1); + rsize = vsize; + } + if (cy_limb == 0) + { + tp[rsize] = 1; + rsize++; + exp++; + } + goto normalize; + } + } + +general_case: + /* If U extends beyond PREC, ignore the part that does. */ + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + + /* If V extends beyond PREC, ignore the part that does. + Note that this may make vsize negative. */ + if (vsize + ediff > prec) + { + vp += vsize + ediff - prec; + vsize = prec - ediff; + } + + /* Allocate temp space for the result. Allocate + just vsize + ediff later??? */ + tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB); + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (tp != up) + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + /* Locate the least significant non-zero limb in (the needed + parts of) U and V, to simplify the code below. */ + for (;;) + { + if (vsize == 0) + { + MPN_COPY (rp, up, usize); + rsize = usize; + goto done; + } + if (vp[0] != 0) + break; + vp++, vsize--; + } + for (;;) + { + if (usize == 0) + { + MPN_COPY (rp, vp, vsize); + rsize = vsize; + negate ^= 1; + goto done; + } + if (up[0] != 0) + break; + up++, usize--; + } + + /* uuuu | uuuu | uuuu | uuuu | uuuu */ + /* vvvvvvv | vv | vvvvv | v | vv */ + + if (usize > ediff) + { + /* U and V partially overlaps. */ + if (ediff == 0) + { + /* Have to compare the leading limbs of u and v + to determine whether to compute u - v or v - u. */ + if (usize >= vsize) + { + /* uuuu */ + /* vv */ + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else /* (usize < vsize) */ + { + /* uuuu */ + /* vvvvvvv */ + mp_size_t size, i; + size = vsize - usize; + tp[0] = -vp[0]; + for (i = 1; i < size; i++) + tp[i] = ~vp[i]; + mpn_sub_n (tp + size, up, vp + size, usize); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize; + } + } + else + { + if (vsize + ediff <= usize) + { + /* uuuu */ + /* v */ + mp_size_t size; + size = usize - ediff - vsize; + MPN_COPY (tp, up, size); + mpn_sub (tp + size, up + size, usize - size, vp, vsize); + rsize = usize; + } + else + { + /* uuuu */ + /* vvvvv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0]; + for (i = 1; i < size; i++) + tp[i] = ~vp[i]; + mpn_sub (tp + size, up, usize, vp + size, usize - ediff); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize + ediff; + } + } + } + else + { + /* uuuu */ + /* vv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0]; + for (i = 1; i < vsize; i++) + tp[i] = ~vp[i]; + for (i = vsize; i < size; i++) + tp[i] = ~(mp_limb_t) 0; + mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1); + rsize = size + usize; + } + + normalize: + /* Full normalize. Optimize later. */ + while (rsize != 0 && tp[rsize - 1] == 0) + { + rsize--; + exp--; + } + MPN_COPY (rp, tp, rsize); + } + + done: + r->_mp_size = negate ? -rsize : rsize; + r->_mp_exp = exp; + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/sub_ui.c b/gnu/lib/libgmp/mpf/sub_ui.c new file mode 100644 index 00000000000..7806213788d --- /dev/null +++ b/gnu/lib/libgmp/mpf/sub_ui.c @@ -0,0 +1,49 @@ +/* mpf_sub_ui -- Subtract an unsigned integer from a float. + +Copyright (C) 1993, 1994, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_sub_ui (mpf_ptr sum, mpf_srcptr u, unsigned long int v) +#else +mpf_sub_ui (sum, u, v) + mpf_ptr sum; + mpf_srcptr u; + unsigned long int v; +#endif +{ + __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/gnu/lib/libgmp/mpf/tests/Makefile.in b/gnu/lib/libgmp/mpf/tests/Makefile.in new file mode 100644 index 00000000000..e108b9e6800 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/Makefile.in @@ -0,0 +1,75 @@ +# Makefile for mpf/tests for GNU MP + +srcdir = . + +CC = gcc + +TEST_LIBS = ../../libgmp.a +INCLUDES = -I../../mpn -I$(srcdir)/../.. +CFLAGS = -g -O + +.c.o: + $(CC) -c $(INCLUDES) $(CFLAGS) $(XCFLAGS) $< + +TEST_SRCS = t-add.c t-sub.c t-conv.c t-sqrt.c t-muldiv.c t-dm2exp.c t-reuse.c ref.c +TEST_OBJS = t-add.o t-sub.o t-conv.o t-sqrt.o t-muldiv.o t-dm2exp.c reuse.o +TESTS = t-add t-sub t-conv t-sqrt t-muldiv t-dm2exp reuse + +check: Makefile st-add st-sub st-conv st-sqrt st-muldiv st-dm2exp st-reuse + @echo "The tests passed." + +st-add: t-add + ./t-add + touch $@ +st-sub: t-sub + ./t-sub + touch $@ +st-conv: t-conv + ./t-conv + touch $@ +st-sqrt: t-sqrt + ./t-sqrt + touch $@ +st-muldiv: t-muldiv + ./t-muldiv + touch $@ +st-dm2exp: t-dm2exp + ./t-dm2exp + touch $@ +st-reuse: reuse + ./reuse + touch $@ + +H = $(srcdir)/../../gmp.h $(srcdir)/../../gmp-impl.h \ + $(srcdir)/../../urandom.h ../../mpn/gmp-mparam.h + +t-add: t-add.o ref.o $(TEST_LIBS) + $(CC) -o $@ $@.o ref.o $(TEST_LIBS) $(CFLAGS) +t-sub: t-sub.o ref.o $(TEST_LIBS) + $(CC) -o $@ $@.o ref.o $(TEST_LIBS) $(CFLAGS) +t-conv: t-conv.o $(TEST_LIBS) + $(CC) -o $@ $@.o $(TEST_LIBS) $(CFLAGS) +t-sqrt: t-sqrt.o $(TEST_LIBS) + $(CC) -o $@ $@.o $(TEST_LIBS) $(CFLAGS) +t-muldiv: t-muldiv.o $(TEST_LIBS) + $(CC) -o $@ $@.o $(TEST_LIBS) $(CFLAGS) +t-dm2exp: t-dm2exp.o $(TEST_LIBS) + $(CC) -o $@ $@.o $(TEST_LIBS) $(CFLAGS) +reuse: reuse.o $(TEST_LIBS) + $(CC) -o $@ $@.o $(TEST_LIBS) $(CFLAGS) + +t-add.o: $(srcdir)/t-add.c +t-sub.o: $(srcdir)/t-sub.c +t-conv.o: $(srcdir)/t-conv.c +t-sqrt.o: $(srcdir)/t-sqrt.c +t-muldiv.o: $(srcdir)/t-muldiv.c +t-dm2exp.o: $(srcdir)/t-dm2exp.c +reuse.o: $(srcdir)/reuse.c + +clean mostlyclean: + rm -f *.o st-* $(TESTS) +distclean maintainer-clean: clean + rm -f Makefile config.status + +Makefile: $(srcdir)/Makefile.in + $(SHELL) ./config.status diff --git a/gnu/lib/libgmp/mpf/tests/configure.in b/gnu/lib/libgmp/mpf/tests/configure.in new file mode 100644 index 00000000000..319219cd75a --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/configure.in @@ -0,0 +1,11 @@ +# This file is a shell script that supplies the information necessary +# to tailor a template configure script into the configure script +# appropriate for this directory. For more information, check any +# existing configure script. + +srctrigger=t-add.c +srcname="gmp/mpf/tests" + +# per-host: + +# per-target: diff --git a/gnu/lib/libgmp/mpf/tests/ref.c b/gnu/lib/libgmp/mpf/tests/ref.c new file mode 100644 index 00000000000..b2c3e6c2011 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/ref.c @@ -0,0 +1,203 @@ +/* Reference floating point routines. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +#if __STDC__ +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); +#else +void ref_mpf_add (); +void ref_mpf_sub (); +#endif + +void +#if __STDC__ +ref_mpf_add (mpf_t w, const mpf_t u, const mpf_t v) +#else +ref_mpf_add (w, u, v) + mpf_t w; + const mpf_t u; + const mpf_t v; +#endif +{ + mp_size_t hi, lo, size; + mp_ptr ut, vt, wt; + int neg; + mp_exp_t exp; + mp_limb_t cy; + TMP_DECL (mark); + + TMP_MARK (mark); + + if (SIZ (u) == 0) + { + size = ABSIZ (v); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (v), size); + exp = EXP (v); + neg = SIZ (v) < 0; + goto done; + } + if (SIZ (v) == 0) + { + size = ABSIZ (u); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (u), size); + exp = EXP (u); + neg = SIZ (u) < 0; + goto done; + } + if ((SIZ (u) ^ SIZ (v)) < 0) + { + mpf_t tmp; + SIZ (tmp) = -SIZ (v); + EXP (tmp) = EXP (v); + PTR (tmp) = PTR (v); + ref_mpf_sub (w, u, tmp); + return; + } + neg = SIZ (u) < 0; + + /* Compute the significance of the hi and lo end of the result. */ + hi = MAX (EXP (u), EXP (v)); + lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); + size = hi - lo; + ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (ut, size); + MPN_ZERO (vt, size); + {int off; + off = size + (EXP (u) - hi) - ABSIZ (u); + MPN_COPY (ut + off, PTR (u), ABSIZ (u)); + off = size + (EXP (v) - hi) - ABSIZ (v); + MPN_COPY (vt + off, PTR (v), ABSIZ (v)); + } + + cy = mpn_add_n (wt, ut, vt, size); + wt[size] = cy; + size += cy; + exp = hi + cy; + +done: + if (size > PREC (w)) + { + wt += size - PREC (w); + size = PREC (w); + } + MPN_COPY (PTR (w), wt, size); + SIZ (w) = neg == 0 ? size : -size; + EXP (w) = exp; + TMP_FREE (mark); +} + +void +#if __STDC__ +ref_mpf_sub (mpf_t w, const mpf_t u, const mpf_t v) +#else +ref_mpf_sub (w, u, v) + mpf_t w; + const mpf_t u; + const mpf_t v; +#endif +{ + mp_size_t hi, lo, size; + mp_ptr ut, vt, wt; + int neg; + mp_exp_t exp; + TMP_DECL (mark); + + TMP_MARK (mark); + + if (SIZ (u) == 0) + { + size = ABSIZ (v); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (v), size); + exp = EXP (v); + neg = SIZ (v) > 0; + goto done; + } + if (SIZ (v) == 0) + { + size = ABSIZ (u); + wt = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB); + MPN_COPY (wt, PTR (u), size); + exp = EXP (u); + neg = SIZ (u) < 0; + goto done; + } + if ((SIZ (u) ^ SIZ (v)) < 0) + { + mpf_t tmp; + SIZ (tmp) = -SIZ (v); + EXP (tmp) = EXP (v); + PTR (tmp) = PTR (v); + ref_mpf_add (w, u, tmp); + if (SIZ (u) < 0) + mpf_neg (w, w); + return; + } + neg = SIZ (u) < 0; + + /* Compute the significance of the hi and lo end of the result. */ + hi = MAX (EXP (u), EXP (v)); + lo = MIN (EXP (u) - ABSIZ (u), EXP (v) - ABSIZ (v)); + size = hi - lo; + ut = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (ut, size); + MPN_ZERO (vt, size); + {int off; + off = size + (EXP (u) - hi) - ABSIZ (u); + MPN_COPY (ut + off, PTR (u), ABSIZ (u)); + off = size + (EXP (v) - hi) - ABSIZ (v); + MPN_COPY (vt + off, PTR (v), ABSIZ (v)); + } + + if (mpn_cmp (ut, vt, size) >= 0) + mpn_sub_n (wt, ut, vt, size); + else + { + mpn_sub_n (wt, vt, ut, size); + neg ^= 1; + } + exp = hi; + while (size != 0 && wt[size - 1] == 0) + { + size--; + exp--; + } + +done: + if (size > PREC (w)) + { + wt += size - PREC (w); + size = PREC (w); + } + MPN_COPY (PTR (w), wt, size); + SIZ (w) = neg == 0 ? size : -size; + EXP (w) = exp; + TMP_FREE (mark); +} diff --git a/gnu/lib/libgmp/mpf/tests/reuse.c b/gnu/lib/libgmp/mpf/tests/reuse.c new file mode 100644 index 00000000000..6829cc46de1 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/reuse.c @@ -0,0 +1,186 @@ +/* Test that routines allow reusing a source variable as destination. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +#ifndef EXPO +#define EXPO 32 +#endif + +#if __STDC__ +typedef void (*dss_func) (mpf_ptr, mpf_srcptr, mpf_srcptr); +#else +typedef void (*dss_func) (); +#endif + +dss_func dss_funcs[] = +{ + mpf_div, mpf_add, mpf_mul, mpf_sub, +}; + +char *dss_func_names[] = +{ + "mpf_div", "mpf_add", "mpf_mul", "mpf_sub", +}; + +#if __STDC__ +typedef void (*dsi_func) (mpf_ptr, mpf_srcptr, unsigned long int); +#else +typedef void (*dsi_func) (); +#endif + +dsi_func dsi_funcs[] = +{ + mpf_div_ui, mpf_add_ui, mpf_mul_ui, mpf_sub_ui, +}; + +char *dsi_func_names[] = +{ + "mpf_div_ui", "mpf_add_ui", "mpf_mul_ui", "mpf_sub_ui", +}; + +#if __STDC__ +typedef void (*dis_func) (mpf_ptr, unsigned long int, mpf_srcptr); +#else +typedef void (*dis_func) (); +#endif + +dis_func dis_funcs[] = +{ + mpf_ui_div, mpf_ui_sub, +}; + +char *dis_func_names[] = +{ + "mpf_ui_div", "mpf_ui_sub", +}; + +main (argc, argv) + int argc; + char **argv; +{ + int i; + int pass, reps = 100000; + mpf_t in1, in2, out1; + unsigned long int in1i, in2i; + mpf_t res1, res2, res3; + mp_size_t bprec = 100; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init (in1); + mpf_init (in2); + mpf_init (out1); + mpf_init (res1); + mpf_init (res2); + mpf_init (res3); + + for (pass = 1; pass <= reps; pass++) + { + mpf_random2 (in1, urandom () % SIZE - SIZE/2, urandom () % EXPO); + mpf_random2 (in2, urandom () % SIZE - SIZE/2, urandom () % EXPO); + + for (i = 0; i < sizeof (dss_funcs) / sizeof (dss_func); i++) + { + /* Don't divide by 0. */ + if (i == 0 && mpf_cmp_ui (in2, 0) == 0) + continue; + + (dss_funcs[i]) (res1, in1, in2); + + mpf_set (out1, in1); + (dss_funcs[i]) (out1, out1, in2); + mpf_set (res2, out1); + + mpf_set (out1, in2); + (dss_funcs[i]) (out1, in1, out1); + mpf_set (res3, out1); + + if (mpf_cmp (res1, res2) != 0) + dump_abort (dss_func_names[i], res1, res2); + if (mpf_cmp (res1, res3) != 0) + dump_abort (dss_func_names[i], res1, res3); + } + + in2i = urandom (); + for (i = 0; i < sizeof (dsi_funcs) / sizeof (dsi_func); i++) + { + /* Don't divide by 0. */ + if (i == 0 && in2i == 0) + continue; + + (dsi_funcs[i]) (res1, in1, in2i); + + mpf_set (out1, in1); + (dsi_funcs[i]) (out1, out1, in2i); + mpf_set (res2, out1); + + if (mpf_cmp (res1, res2) != 0) + dump_abort (dsi_func_names[i], res1, res2); + } + + in1i = urandom (); + for (i = 0; i < sizeof (dis_funcs) / sizeof (dis_func); i++) + { + /* Don't divide by 0. */ + if (i == 0 && mpf_cmp_ui (in2, 0) == 0) + continue; + + (dis_funcs[i]) (res1, in1i, in2); + + mpf_set (out1, in2); + (dis_funcs[i]) (out1, in1i, in2); + mpf_set (res2, out1); + + if (mpf_cmp (res1, res2) != 0) + dump_abort (dis_func_names[i], res1, res2); + } + + } + + exit (0); +} + +dump_abort (name, res1, res2) + char *name; + mpf_t res1, res2; +{ + printf ("failure in %s:", name); + oo (res1); + oo (res2); + abort (); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} + +#if 0 +void mpf_div_2exp _PROTO ((mpf_ptr, mpf_srcptr, unsigned long int)); +void mpf_mul_2exp _PROTO ((mpf_ptr, mpf_srcptr, unsigned long int)); + +void mpf_abs _PROTO ((mpf_ptr, mpf_srcptr)); +void mpf_sqrt _PROTO ((mpf_ptr, mpf_srcptr)); +void mpf_neg _PROTO ((mpf_ptr, mpf_srcptr)); +#endif diff --git a/gnu/lib/libgmp/mpf/tests/t-add.c b/gnu/lib/libgmp/mpf/tests/t-add.c new file mode 100644 index 00000000000..6f985bb1853 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-add.c @@ -0,0 +1,117 @@ +/* Test mpf_add. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +#if __STDC__ +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); +#else +void ref_mpf_add (); +void ref_mpf_sub (); +#endif + +main (argc, argv) + int argc; + char **argv; +{ + mp_size_t size; + mp_exp_t exp; + int reps = 100000; + int i; + mpf_t u, v, w, wref; + mp_size_t bprec = 100; + mpf_t rerr, max_rerr, limit_rerr; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, bprec); +#if VERBOSE + mpf_dump (limit_rerr); +#endif + mpf_init (rerr); + mpf_init_set_ui (max_rerr, 0); + + mpf_init (u); + mpf_init (v); + mpf_init (w); + mpf_init (wref); + for (i = 0; i < reps; i++) + { + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (u, size, exp); + + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (v, size, exp); + + mpf_add (w, u, v); + ref_mpf_add (wref, u, v); + + mpf_reldiff (rerr, w, wref); + if (mpf_cmp (rerr, max_rerr) > 0) + { + mpf_set (max_rerr, rerr); +#if VERBOSE + mpf_dump (max_rerr); +#endif + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf ("wref = "); mpf_dump (wref); + printf (" w = "); mpf_dump (w); + abort (); + } + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/tests/t-conv.c b/gnu/lib/libgmp/mpf/tests/t-conv.c new file mode 100644 index 00000000000..00d461d064c --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-conv.c @@ -0,0 +1,120 @@ +/* Test mpf_get_str and mpf_set_str. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 10 +#endif + +#ifndef EXPO +#define EXPO 20 +#endif + +main (argc, argv) + int argc; + char **argv; +{ + mpf_t x, y; + int reps = 20000; + int i; + mp_size_t bprec = 100; + mpf_t d, rerr, max_rerr, limit_rerr; + char *str; + long bexp; + long size, exp; + int base; + char buf[SIZE * BITS_PER_MP_LIMB + 5]; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, bprec); +#if VERBOSE + mpf_dump (limit_rerr); +#endif + mpf_init (rerr); + mpf_init_set_ui (max_rerr, 0); + + mpf_init (x); + mpf_init (y); + mpf_init (d); + + for (i = 0; i < reps; i++) + { + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % EXPO; + mpf_random2 (x, size, exp); + + base = urandom () % 35 + 2; + + str = mpf_get_str (0, &bexp, base, 0, x); + + if (str[0] == '-') + sprintf (buf, "-0.%s@%ld", str + 1, bexp); + else + sprintf (buf, "0.%s@%ld", str, bexp); + + mpf_set_str (y, buf, -base); + free (str); + + mpf_reldiff (rerr, x, y); + if (mpf_cmp (rerr, max_rerr) > 0) + { + mpf_set (max_rerr, rerr); +#if VERBOSE + mpf_dump (max_rerr); +#endif + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR after %d tests\n", i); + printf ("base = %d\n", base); + printf (" x = "); mpf_dump (x); + printf (" y = "); mpf_dump (y); + abort (); + } + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + int i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/tests/t-dm2exp.c b/gnu/lib/libgmp/mpf/tests/t-dm2exp.c new file mode 100644 index 00000000000..a8c4e8f519e --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-dm2exp.c @@ -0,0 +1,101 @@ +/* Test mpf_div, mpf_div_2exp, mpf_mul_2exp. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +main (argc, argv) + int argc; + char **argv; +{ + int reps = 500000; + int i; + mpf_t u, v, w1, w2; + mp_size_t bprec = 100; + mpf_t rerr, limit_rerr; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init (rerr); + mpf_init (limit_rerr); + + mpf_init (u); + mpf_init (v); + mpf_init (w1); + mpf_init (w2); + + for (i = 0; i < reps; i++) + { + mp_size_t res_prec; + unsigned long int pow2; + + res_prec = urandom () % (bprec + 100); + mpf_set_prec (w1, res_prec); + mpf_set_prec (w2, res_prec); + + mpf_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, res_prec); + + pow2 = urandom () % 0x10000; + mpf_set_ui (v, 1); + mpf_mul_2exp (v, v, pow2); + + mpf_div_2exp (w1, u, pow2); + mpf_div (w2, u, v); + mpf_reldiff (rerr, w1, w2); + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR in mpf_div or mpf_div_2exp after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf (" w1 = "); mpf_dump (w1); + printf (" w2 = "); mpf_dump (w2); + abort (); + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/tests/t-muldiv.c b/gnu/lib/libgmp/mpf/tests/t-muldiv.c new file mode 100644 index 00000000000..d1721e8db78 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-muldiv.c @@ -0,0 +1,161 @@ +/* Test mpf_mul, mpf_div, mpf_ui_div, and mpf_div_ui. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +main (argc, argv) + int argc; + char **argv; +{ + mp_size_t size; + mp_exp_t exp; + int reps = 100000; + int i; + mpf_t u, v, w, x; + mp_size_t bprec = 100; + mpf_t rerr, limit_rerr; + mp_limb_t ulimb, vlimb; + int single_flag; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init (rerr); + mpf_init (limit_rerr); + + mpf_init (u); + mpf_init (v); + mpf_init (w); + mpf_init (x); + + for (i = 0; i < reps; i++) + { + mp_size_t res_prec; + + res_prec = urandom () % (bprec + 100); + mpf_set_prec (w, res_prec); + mpf_set_prec (x, res_prec); + + mpf_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, res_prec); + + single_flag = 0; + + if ((urandom () & 1) != 0) + { + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (u, size, exp); + } + else + { + ulimb = urandom (); + mpf_set_ui (u, ulimb); + single_flag = 1; + } + + if ((urandom () & 1) != 0) + { + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (v, size, exp); + } + else + { + vlimb = urandom (); + mpf_set_ui (v, vlimb); + single_flag = 2; + } + + if (mpf_sgn (v) == 0) + continue; + + mpf_div (w, u, v); + mpf_mul (x, w, v); + mpf_reldiff (rerr, u, x); + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR in mpf_mul or mpf_div after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf (" x = "); mpf_dump (x); + printf (" w = "); mpf_dump (w); + abort (); + } + + if (single_flag == 2) + { + mpf_div_ui (x, u, vlimb); + mpf_reldiff (rerr, w, x); + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR in mpf_div or mpf_div_ui after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf (" x = "); mpf_dump (x); + printf (" w = "); mpf_dump (w); + abort (); + } + } + + if (single_flag == 1) + { + mpf_ui_div (x, ulimb, v); + mpf_reldiff (rerr, w, x); + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR in mpf_div or mpf_ui_div after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf (" x = "); mpf_dump (x); + printf (" w = "); mpf_dump (w); + abort (); + } + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/tests/t-sqrt.c b/gnu/lib/libgmp/mpf/tests/t-sqrt.c new file mode 100644 index 00000000000..5c72ab5fefe --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-sqrt.c @@ -0,0 +1,103 @@ +/* Test mpf_sqrt, mpf_mul. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +main (argc, argv) + int argc; + char **argv; +{ + mp_size_t size; + mp_exp_t exp; + int reps = 100000; + int i; + mpf_t x, y, y2; + mp_size_t bprec = 100; + mpf_t rerr, max_rerr, limit_rerr; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, bprec); +#if VERBOSE + mpf_dump (limit_rerr); +#endif + mpf_init (rerr); + mpf_init_set_ui (max_rerr, 0); + + mpf_init (x); + mpf_init (y); + mpf_init (y2); + for (i = 0; i < reps; i++) + { + size = urandom () % SIZE; + exp = urandom () % SIZE; + mpf_random2 (x, size, exp); + + mpf_sqrt (y, x); + mpf_mul (y2, y, y); + + mpf_reldiff (rerr, x, y2); + if (mpf_cmp (rerr, max_rerr) > 0) + { + mpf_set (max_rerr, rerr); +#if VERBOSE + mpf_dump (max_rerr); +#endif + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR after %d tests\n", i); + printf (" x = "); mpf_dump (x); + printf (" y = "); mpf_dump (y); + printf (" y2 = "); mpf_dump (y2); + abort (); + } + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/tests/t-sub.c b/gnu/lib/libgmp/mpf/tests/t-sub.c new file mode 100644 index 00000000000..69188665577 --- /dev/null +++ b/gnu/lib/libgmp/mpf/tests/t-sub.c @@ -0,0 +1,122 @@ +/* Test mpf_sub. + +Copyright (C) 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include <stdio.h> +#include "gmp.h" +#include "gmp-impl.h" +#include "urandom.h" + +#ifndef SIZE +#define SIZE 16 +#endif + +#if __STDC__ +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); +#else +void ref_mpf_add (); +void ref_mpf_sub (); +#endif + +main (argc, argv) + int argc; + char **argv; +{ + mp_size_t size; + mp_exp_t exp; + int reps = 500000; + int i; + mpf_t u, v, w, wref; + mp_size_t bprec = 100; + mpf_t rerr, max_rerr, limit_rerr; + + if (argc > 1) + { + reps = strtol (argv[1], 0, 0); + if (argc > 2) + bprec = strtol (argv[2], 0, 0); + } + + mpf_set_default_prec (bprec); + + mpf_init_set_ui (limit_rerr, 1); + mpf_div_2exp (limit_rerr, limit_rerr, bprec); +#if VERBOSE + mpf_dump (limit_rerr); +#endif + mpf_init (rerr); + mpf_init_set_ui (max_rerr, 0); + + mpf_init (u); + mpf_init (v); + mpf_init (w); + mpf_init (wref); + for (i = 0; i < reps; i++) + { + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (u, size, exp); + + size = urandom () % (2 * SIZE) - SIZE; + exp = urandom () % SIZE; + mpf_random2 (v, size, exp); + + if ((urandom () & 1) != 0) + mpf_add_ui (u, v, 1); + else if ((urandom () & 1) != 0) + mpf_sub_ui (u, v, 1); + + mpf_sub (w, u, v); + ref_mpf_sub (wref, u, v); + + mpf_reldiff (rerr, w, wref); + if (mpf_cmp (rerr, max_rerr) > 0) + { + mpf_set (max_rerr, rerr); +#if VERBOSE + mpf_dump (max_rerr); +#endif + if (mpf_cmp (rerr, limit_rerr) > 0) + { + printf ("ERROR after %d tests\n", i); + printf (" u = "); mpf_dump (u); + printf (" v = "); mpf_dump (v); + printf ("wref = "); mpf_dump (wref); + printf (" w = "); mpf_dump (w); + abort (); + } + } + } + + exit (0); +} + +oo (x) + mpf_t x; +{ + mp_size_t i; + printf (" exp = %ld\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08lX ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/gnu/lib/libgmp/mpf/ui_div.c b/gnu/lib/libgmp/mpf/ui_div.c new file mode 100644 index 00000000000..a49a952048b --- /dev/null +++ b/gnu/lib/libgmp/mpf/ui_div.c @@ -0,0 +1,131 @@ +/* mpf_ui_div -- Divide an unsigned integer with a float. + +Copyright (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +void +#if __STDC__ +mpf_ui_div (mpf_ptr r, unsigned long int u, mpf_srcptr v) +#else +mpf_ui_div (r, u, v) + mpf_ptr r; + unsigned long int u; + mpf_srcptr v; +#endif +{ + mp_srcptr vp; + mp_ptr rp, tp; + mp_size_t vsize; + mp_size_t rsize, tsize; + mp_size_t sign_quotient; + mp_size_t prec; + unsigned normalization_steps; + mp_limb_t q_limb; + mp_exp_t rexp; + TMP_DECL (marker); + + vsize = v->_mp_size; + sign_quotient = vsize; + vsize = ABS (vsize); + prec = r->_mp_prec; + + if (vsize == 0) + vsize = 1 / vsize; /* divide by zero as directed */ + if (u == 0) + { + r->_mp_size = 0; + r->_mp_exp = 0; + return; + } + + TMP_MARK (marker); + rexp = 1 - v->_mp_exp; + + rp = r->_mp_d; + vp = v->_mp_d; + + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + + tsize = vsize + prec; + tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (tp, tsize); + + count_leading_zeros (normalization_steps, vp[vsize - 1]); + + /* Normalize the divisor and the dividend. */ + if (normalization_steps != 0) + { + mp_ptr tmp; + mp_limb_t dividend_high, dividend_low; + + /* Shift up the divisor setting the most significant bit of + the most significant limb. Use temporary storage not to clobber + the original contents of the divisor. */ + tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + mpn_lshift (tmp, vp, vsize, normalization_steps); + vp = tmp; + + /* Shift up the dividend, possibly introducing a new most + significant word. */ + dividend_high = (mp_limb_t) u >> (BITS_PER_MP_LIMB - normalization_steps); + dividend_low = (mp_limb_t) u << normalization_steps; + + tp[tsize - 1] = dividend_low; + if (dividend_high != 0) + { + tp[tsize] = dividend_high; + tsize++; + rexp++; + } + } + else + { + /* The divisor is already normalized, as required. + Copy it to temporary space if it overlaps with the quotient. */ + if (vp - rp <= tsize - vsize) + { + mp_ptr tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + MPN_COPY (tmp, vp, vsize); + vp = (mp_srcptr) tmp; + } + + tp[tsize - 1] = u; + } + + q_limb = mpn_divmod (rp, tp, tsize, vp, vsize); + rsize = tsize - vsize; + if (q_limb) + { + rp[rsize] = q_limb; + rsize++; + rexp++; + } + + r->_mp_size = sign_quotient >= 0 ? rsize : -rsize; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/gnu/lib/libgmp/mpf/ui_sub.c b/gnu/lib/libgmp/mpf/ui_sub.c new file mode 100644 index 00000000000..acb9210a641 --- /dev/null +++ b/gnu/lib/libgmp/mpf/ui_sub.c @@ -0,0 +1,334 @@ +/* mpf_ui_sub -- Subtract a float from an unsigned long int. + +Copyright (C) 1993, 1994, 1995, 1996 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 the GNU Library General Public License as published by +the Free Software Foundation; either version 2 of the License, or (at your +option) any later version. + +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 Library General Public +License for more details. + +You should have received a copy of the GNU Library General Public License +along with the GNU MP Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" + +void +#if __STDC__ +mpf_ui_sub (mpf_ptr r, unsigned long int u, mpf_srcptr v) +#else +mpf_ui_sub (r, u, v) + mpf_ptr r; + unsigned long int u; + mpf_srcptr v; +#endif +{ + 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 (marker); + + vsize = v->_mp_size; + + /* Handle special cases that don't work in generic code below. */ + if (u == 0) + { + mpf_neg (r, v); + return; + } + if (vsize == 0) + { + mpf_set_ui (r, u); + return; + } + + /* If signs of U and V are different, perform addition. */ + if (vsize < 0) + { + __mpf_struct v_negated; + v_negated._mp_size = -vsize; + v_negated._mp_exp = v->_mp_exp; + v_negated._mp_d = v->_mp_d; + mpf_add_ui (r, &v_negated, u); + return; + } + + TMP_MARK (marker); + + /* Signs are now known to be the same. */ + + ulimb = u; + /* Make U be the operand with the largest exponent. */ + if (1 < v->_mp_exp) + { + negate = 1; + usize = ABS (vsize); + vsize = 1; + up = v->_mp_d; + vp = &ulimb; + rp = r->_mp_d; + prec = r->_mp_prec + 1; + uexp = v->_mp_exp; + ediff = uexp - 1; + } + else + { + negate = 0; + usize = 1; + vsize = ABS (vsize); + up = &ulimb; + vp = v->_mp_d; + rp = r->_mp_d; + prec = r->_mp_prec; + uexp = 1; + ediff = 1 - v->_mp_exp; + } + + /* Ignore leading limbs in U and V that are equal. Doing + this helps increase the precision of the result. */ + if (ediff == 0) + { + /* This loop normally exits immediately. Optimize for that. */ + for (;;) + { + usize--; + vsize--; + if (up[usize] != vp[vsize]) + break; + uexp--; + if (usize == 0) + goto Lu0; + if (vsize == 0) + goto Lv0; + } + usize++; + vsize++; + /* Note that either operand (but not both operands) might now have + leading zero limbs. It matters only that U is unnormalized if + vsize is now zero, and vice versa. And it is only in that case + that we have to adjust uexp. */ + if (vsize == 0) + Lv0: + while (usize != 0 && up[usize - 1] == 0) + usize--, uexp--; + if (usize == 0) + Lu0: + while (vsize != 0 && vp[vsize - 1] == 0) + vsize--, uexp--; + } + + /* If U extends beyond PREC, ignore the part that does. */ + if (usize > prec) + { + up += usize - prec; + usize = prec; + } + + /* If V extends beyond PREC, ignore the part that does. + Note that this may make vsize negative. */ + if (vsize + ediff > prec) + { + vp += vsize + ediff - prec; + vsize = prec - ediff; + } + + /* Allocate temp space for the result. Allocate + just vsize + ediff later??? */ + tp = (mp_ptr) TMP_ALLOC (prec * BYTES_PER_MP_LIMB); + + if (ediff >= prec) + { + /* V completely cancelled. */ + if (tp != up) + MPN_COPY (rp, up, usize); + rsize = usize; + } + else + { + /* Locate the least significant non-zero limb in (the needed + parts of) U and V, to simplify the code below. */ + for (;;) + { + if (vsize == 0) + { + MPN_COPY (rp, up, usize); + rsize = usize; + goto done; + } + if (vp[0] != 0) + break; + vp++, vsize--; + } + for (;;) + { + if (usize == 0) + { + MPN_COPY (rp, vp, vsize); + rsize = vsize; + negate ^= 1; + goto done; + } + if (up[0] != 0) + break; + up++, usize--; + } + + /* uuuu | uuuu | uuuu | uuuu | uuuu */ + /* vvvvvvv | vv | vvvvv | v | vv */ + + if (usize > ediff) + { + /* U and V partially overlaps. */ + if (ediff == 0) + { + /* Have to compare the leading limbs of u and v + to determine whether to compute u - v or v - u. */ + if (usize > vsize) + { + /* uuuu */ + /* vv */ + int cmp; + cmp = mpn_cmp (up + usize - vsize, vp, vsize); + if (cmp >= 0) + { + mp_size_t size; + size = usize - vsize; + MPN_COPY (tp, up, size); + mpn_sub_n (tp + size, up + size, vp, vsize); + rsize = usize; + } + else + { + /* vv */ /* Swap U and V. */ + /* uuuu */ + mp_size_t size, i; + size = usize - vsize; + tp[0] = -up[0]; + for (i = 1; i < size; i++) + tp[i] = ~up[i]; + mpn_sub_n (tp + size, vp, up + size, vsize); + mpn_sub_1 (tp + size, tp + size, vsize, (mp_limb_t) 1); + negate ^= 1; + rsize = usize; + } + } + else if (usize < vsize) + { + /* uuuu */ + /* vvvvvvv */ + int cmp; + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp > 0) + { + mp_size_t size, i; + size = vsize - usize; + tp[0] = -vp[0]; + for (i = 1; i < size; i++) + tp[i] = ~vp[i]; + mpn_sub_n (tp + size, up, vp + size, usize); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize; + } + else + { + /* vvvvvvv */ /* Swap U and V. */ + /* uuuu */ + /* This is the only place we can get 0.0. */ + mp_size_t size; + size = vsize - usize; + MPN_COPY (tp, vp, size); + mpn_sub_n (tp + size, vp + size, up, usize); + negate ^= 1; + rsize = vsize; + } + } + else + { + /* uuuu */ + /* vvvv */ + int cmp; + cmp = mpn_cmp (up, vp + vsize - usize, usize); + if (cmp > 0) + { + mpn_sub_n (tp, up, vp, usize); + rsize = usize; + } + else + { + mpn_sub_n (tp, vp, up, usize); + negate ^= 1; + rsize = usize; + /* can give zero */ + } + } + } + else + { + if (vsize + ediff <= usize) + { + /* uuuu */ + /* v */ + mp_size_t size; + size = usize - ediff - vsize; + MPN_COPY (tp, up, size); + mpn_sub (tp + size, up + size, usize - size, vp, vsize); + rsize = usize; + } + else + { + /* uuuu */ + /* vvvvv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0]; + for (i = 1; i < size; i++) + tp[i] = ~vp[i]; + mpn_sub (tp + size, up, usize, vp + size, usize - ediff); + mpn_sub_1 (tp + size, tp + size, usize, (mp_limb_t) 1); + rsize = vsize + ediff; + } + } + } + else + { + /* uuuu */ + /* vv */ + mp_size_t size, i; + size = vsize + ediff - usize; + tp[0] = -vp[0]; + for (i = 1; i < vsize; i++) + tp[i] = ~vp[i]; + for (i = vsize; i < size; i++) + tp[i] = ~(mp_limb_t) 0; + mpn_sub_1 (tp + size, up, usize, (mp_limb_t) 1); + rsize = size + usize; + } + + /* Full normalize. Optimize later. */ + while (rsize != 0 && tp[rsize - 1] == 0) + { + rsize--; + uexp--; + } + MPN_COPY (rp, tp, rsize); + } + + done: + r->_mp_size = negate ? -rsize : rsize; + r->_mp_exp = uexp; + TMP_FREE (marker); +} |