diff options
Diffstat (limited to 'mpf')
| -rw-r--r-- | mpf/Makefile.in | 115 | ||||
| -rw-r--r-- | mpf/abs.c | 56 | ||||
| -rw-r--r-- | mpf/add.c | 180 | ||||
| -rw-r--r-- | mpf/add_ui.c | 142 | ||||
| -rw-r--r-- | mpf/clear.c | 35 | ||||
| -rw-r--r-- | mpf/cmp.c | 114 | ||||
| -rw-r--r-- | mpf/cmp_si.c | 98 | ||||
| -rw-r--r-- | mpf/cmp_ui.c | 80 | ||||
| -rw-r--r-- | mpf/configure.in | 12 | ||||
| -rw-r--r-- | mpf/div.c | 143 | ||||
| -rw-r--r-- | mpf/div_2exp.c | 79 | ||||
| -rw-r--r-- | mpf/div_ui.c | 78 | ||||
| -rw-r--r-- | mpf/dump.c | 43 | ||||
| -rw-r--r-- | mpf/eq.c | 121 | ||||
| -rw-r--r-- | mpf/get_prc.c | 34 | ||||
| -rw-r--r-- | mpf/get_str.c | 500 | ||||
| -rw-r--r-- | mpf/init.c | 38 | ||||
| -rw-r--r-- | mpf/init2.c | 41 | ||||
| -rw-r--r-- | mpf/inp_str.c | 89 | ||||
| -rw-r--r-- | mpf/iset.c | 59 | ||||
| -rw-r--r-- | mpf/iset_d.c | 39 | ||||
| -rw-r--r-- | mpf/iset_si.c | 55 | ||||
| -rw-r--r-- | mpf/iset_str.c | 40 | ||||
| -rw-r--r-- | mpf/iset_ui.c | 40 | ||||
| -rw-r--r-- | mpf/mul.c | 94 | ||||
| -rw-r--r-- | mpf/mul_2exp.c | 89 | ||||
| -rw-r--r-- | mpf/mul_ui.c | 72 | ||||
| -rw-r--r-- | mpf/neg.c | 59 | ||||
| -rw-r--r-- | mpf/out_str.c | 89 | ||||
| -rw-r--r-- | mpf/random2.c | 65 | ||||
| -rw-r--r-- | mpf/reldiff.c | 52 | ||||
| -rw-r--r-- | mpf/set.c | 53 | ||||
| -rw-r--r-- | mpf/set_d.c | 97 | ||||
| -rw-r--r-- | mpf/set_dfl_prec.c | 40 | ||||
| -rw-r--r-- | mpf/set_prc.c | 57 | ||||
| -rw-r--r-- | mpf/set_prc_raw.c | 39 | ||||
| -rw-r--r-- | mpf/set_si.c | 51 | ||||
| -rw-r--r-- | mpf/set_str.c | 302 | ||||
| -rw-r--r-- | mpf/set_ui.c | 45 | ||||
| -rw-r--r-- | mpf/size.c | 35 | ||||
| -rw-r--r-- | mpf/sqrt.c | 79 | ||||
| -rw-r--r-- | mpf/sqrt_ui.c | 65 | ||||
| -rw-r--r-- | mpf/sub.c | 402 | ||||
| -rw-r--r-- | mpf/sub_ui.c | 49 | ||||
| -rw-r--r-- | mpf/tests/Makefile.in | 58 | ||||
| -rw-r--r-- | mpf/tests/configure.in | 11 | ||||
| -rw-r--r-- | mpf/tests/ref.c | 184 | ||||
| -rw-r--r-- | mpf/tests/t-add.c | 109 | ||||
| -rw-r--r-- | mpf/tests/t-conv.c | 117 | ||||
| -rw-r--r-- | mpf/tests/t-sqrt.c | 100 | ||||
| -rw-r--r-- | mpf/tests/t-sub.c | 114 | ||||
| -rw-r--r-- | mpf/ui_div.c | 116 | ||||
| -rw-r--r-- | mpf/ui_sub.c | 334 |
53 files changed, 5208 insertions, 0 deletions
diff --git a/mpf/Makefile.in b/mpf/Makefile.in new file mode 100644 index 000000000..0a77fdb63 --- /dev/null +++ b/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 \ + 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_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 \ + 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_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.in $(MPF_OBJS) + rm -f $@ + $(AR) $(AR_FLAGS) $@ $(MPF_OBJS) + +.c.o: + $(CC) -c $(INCLUDES) $(CFLAGS) $(XCFLAGS) $< + +test: libmpf.a + cd tests; $(MAKE) srcdir=$(srcdir) CC=$(CC) SHELL=$(SHELL) + +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_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_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) +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/mpf/abs.c b/mpf/abs.c new file mode 100644 index 000000000..029007a5f --- /dev/null +++ b/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/mpf/add.c b/mpf/add.c new file mode 100644 index 000000000..2db876fc6 --- /dev/null +++ b/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/mpf/add_ui.c b/mpf/add_ui.c new file mode 100644 index 000000000..794f0f1f6 --- /dev/null +++ b/mpf/add_ui.c @@ -0,0 +1,142 @@ +/* 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); + 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 (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/mpf/clear.c b/mpf/clear.c new file mode 100644 index 000000000..beaf4ee3a --- /dev/null +++ b/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/mpf/cmp.c b/mpf/cmp.c new file mode 100644 index 000000000..d440e1152 --- /dev/null +++ b/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/mpf/cmp_si.c b/mpf/cmp_si.c new file mode 100644 index 000000000..01f970812 --- /dev/null +++ b/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/mpf/cmp_ui.c b/mpf/cmp_ui.c new file mode 100644 index 000000000..3a4911ba7 --- /dev/null +++ b/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/mpf/configure.in b/mpf/configure.in new file mode 100644 index 000000000..b6ecf8b3e --- /dev/null +++ b/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/mpf/div.c b/mpf/div.c new file mode 100644 index 000000000..bd2d04dba --- /dev/null +++ b/mpf/div.c @@ -0,0 +1,143 @@ +/* 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 qp, rp, rrp; + mp_size_t usize = u->_mp_size; + mp_size_t vsize = v->_mp_size; + mp_size_t qsize, rsize; + mp_size_t sign_quotient = usize ^ vsize; + unsigned normalization_steps; + mp_limb_t q_limb; + mp_size_t prec; + mp_exp_t rexp; + TMP_DECL (marker); + + prec = r->_mp_prec; + usize = ABS (usize); + vsize = ABS (vsize); + + 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; + + qp = r->_mp_d; + up = u->_mp_d; + vp = v->_mp_d; + + if (vsize > prec) + { + vp += vsize - prec; + vsize = prec; + } + + rsize = vsize + prec; + + rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB); + + if (usize > rsize) + { + up += usize - rsize; + usize = rsize; + rrp = rp; + } + else + { + MPN_ZERO (rp, rsize - usize); + rrp = rp + (rsize - usize); + } + + count_leading_zeros (normalization_steps, vp[vsize - 1]); + + /* Normalize the divisor and the dividend. */ + if (normalization_steps != 0) + { + mp_ptr tp; + mp_limb_t nlimb; + + /* Shift up the divisor setting the most significant bit of + the most significant word. Use temporary storage not to clobber + the original contents of the divisor. */ + tp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + mpn_lshift (tp, vp, vsize, normalization_steps); + vp = tp; + + /* Shift up the dividend, possibly introducing a new most + significant word. Move the shifted dividend in the remainder + meanwhile. */ + nlimb = mpn_lshift (rrp, up, usize, normalization_steps); + if (nlimb != 0) + { + rrp[usize] = nlimb; + rsize++; + rexp++; + } + } + else + { + /* The divisor is already normalized, as required. + Copy it to temporary space if it overlaps with the quotient. */ + if (vp == qp) + { + vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); + MPN_COPY ((mp_ptr) vp, qp, vsize); + } + + /* Move the dividend to the remainder. */ + MPN_COPY (rrp, up, usize); + } + + q_limb = mpn_divmod (qp, rp, rsize, vp, vsize); + + qsize = rsize - vsize; + if (q_limb) + { + qp[qsize] = q_limb; + qsize++; + rexp++; + } + + r->_mp_size = sign_quotient >= 0 ? qsize : -qsize; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/mpf/div_2exp.c b/mpf/div_2exp.c new file mode 100644 index 000000000..d7296254e --- /dev/null +++ b/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/mpf/div_ui.c b/mpf/div_ui.c new file mode 100644 index 000000000..9656395f2 --- /dev/null +++ b/mpf/div_ui.c @@ -0,0 +1,78 @@ +/* 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 qp, rp, rrp; + mp_size_t usize = u->_mp_size; + mp_size_t qsize, rsize; + mp_size_t sign_quotient = usize; + mp_limb_t q_limb; + mp_size_t prec; + mp_exp_t rexp; + TMP_DECL (marker); + + TMP_MARK (marker); + prec = r->_mp_prec; + usize = ABS (usize); + + rsize = prec; + + qp = r->_mp_d; + up = u->_mp_d; + rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB); + if (usize > rsize) + { + up += usize - rsize; + usize = rsize; + rrp = rp; + } + else + { + MPN_ZERO (rp, rsize - usize); + rrp = rp + (rsize - usize); + } + + /* Move the dividend to the remainder. */ + MPN_COPY (rrp, up, usize); + + mpn_divmod_1 (qp, rp, rsize, v); + q_limb = qp[rsize - 1]; + + qsize = rsize - (q_limb == 0); + rexp = u->_mp_exp - (q_limb == 0); + r->_mp_size = sign_quotient >= 0 ? qsize : -qsize; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/mpf/dump.c b/mpf/dump.c new file mode 100644 index 000000000..46d5c05f2 --- /dev/null +++ b/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/mpf/eq.c b/mpf/eq.c new file mode 100644 index 000000000..a474a4dc0 --- /dev/null +++ b/mpf/eq.c @@ -0,0 +1,121 @@ +/* mpf_cmp2 -- 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_cmp2 (mpf_srcptr u, mpf_srcptr v, unsigned long int n_bits) +#else +mpf_cmp2 (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/mpf/get_prc.c b/mpf/get_prc.c new file mode 100644 index 000000000..7f7e41f67 --- /dev/null +++ b/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/mpf/get_str.c b/mpf/get_str.c new file mode 100644 index 000000000..bfee18d07 --- /dev/null +++ b/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/mpf/init.c b/mpf/init.c new file mode 100644 index 000000000..53701218f --- /dev/null +++ b/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/mpf/init2.c b/mpf/init2.c new file mode 100644 index 000000000..a3e5752b2 --- /dev/null +++ b/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/mpf/inp_str.c b/mpf/inp_str.c new file mode 100644 index 000000000..1d3cd4cc2 --- /dev/null +++ b/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/mpf/iset.c b/mpf/iset.c new file mode 100644 index 000000000..c2362e172 --- /dev/null +++ b/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/mpf/iset_d.c b/mpf/iset_d.c new file mode 100644 index 000000000..d09e8b89a --- /dev/null +++ b/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/mpf/iset_si.c b/mpf/iset_si.c new file mode 100644 index 000000000..e67eef4cd --- /dev/null +++ b/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/mpf/iset_str.c b/mpf/iset_str.c new file mode 100644 index 000000000..5c267e7ea --- /dev/null +++ b/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, char *s, int base) +#else +mpf_init_set_str (r, s, base) + mpf_ptr r; + 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/mpf/iset_ui.c b/mpf/iset_ui.c new file mode 100644 index 000000000..c6c24afaf --- /dev/null +++ b/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/mpf/mul.c b/mpf/mul.c new file mode 100644 index 000000000..e5a1f06b3 --- /dev/null +++ b/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/mpf/mul_2exp.c b/mpf/mul_2exp.c new file mode 100644 index 000000000..0ed35a0b1 --- /dev/null +++ b/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/mpf/mul_ui.c b/mpf/mul_ui.c new file mode 100644 index 000000000..8299b4c21 --- /dev/null +++ b/mpf/mul_ui.c @@ -0,0 +1,72 @@ +/* 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; + } + + /* Since we can do it at almost no cost, remove zero limbs at low end of + result. */ + if (up[0] == 0) + up++, size--; + + 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/mpf/neg.c b/mpf/neg.c new file mode 100644 index 000000000..a4139fb90 --- /dev/null +++ b/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/mpf/out_str.c b/mpf/out_str.c new file mode 100644 index 000000000..a1bd328ca --- /dev/null +++ b/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, "@%ld\n", exp); + written += fpret; + } + + TMP_FREE (marker); + return ferror (stream) ? 0 : written; +} diff --git a/mpf/random2.c b/mpf/random2.c new file mode 100644 index 000000000..b09128396 --- /dev/null +++ b/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/mpf/reldiff.c b/mpf/reldiff.c new file mode 100644 index 000000000..97039abc1 --- /dev/null +++ b/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, const mpf_t x, const mpf_t y) +#else +mpf_reldiff (rdiff, x, y) + mpf_t rdiff; + const mpf_t x; + const mpf_t 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/mpf/set.c b/mpf/set.c new file mode 100644 index 000000000..5778b85b0 --- /dev/null +++ b/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/mpf/set_d.c b/mpf/set_d.c new file mode 100644 index 000000000..b82e5f023 --- /dev/null +++ b/mpf/set_d.c @@ -0,0 +1,97 @@ +/* mpf_set_d -- Assign a float from a IEEE 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_set_d (mpf_ptr r, double d) +#else +mpf_set_d (r, d) + mpf_ptr r; + double d; +#endif +{ + mp_ptr rp; + mp_size_t size; + mp_exp_t exp; + mp_limb_t manh, manl; + mp_limb_t man2, man1, man0; + union ieee_double_extract x; + unsigned sc; + + /* ??? This needs more work since it assumes that limbs are 32 bits. + See mpz/set_d for ideas on how to handle 64-bit machines. */ + + if (d == 0) + { + r->_mp_exp = 0; + r->_mp_size = 0; + return; + } + + rp = r->_mp_d; + x.d = d; + + exp = x.s.exp; + sc = (unsigned) (exp + 2) % BITS_PER_MP_LIMB; + + manh = 0x80000000 | (x.s.manh << 11) | (x.s.manl >> 21); + manl = x.s.manl << 11; + if (sc != 0) + { + man2 = manh >> (BITS_PER_MP_LIMB - sc); + man1 = (manl >> (BITS_PER_MP_LIMB - sc)) | (manh << sc); + man0 = manl << sc; + } + else + { + man2 = manh; + man1 = manl; + man0 = 0; + } + + if (man0 == 0) + { + if (man1 == 0) + { + size = 1; + rp[0] = man2; + } + else + { + size = 2; + rp[1] = man2; + rp[0] = man1; + } + } + else + { + size = 3; + rp[2] = man2; + rp[1] = man1; + rp[0] = man0; + } + + r->_mp_exp = (exp + 1) / BITS_PER_MP_LIMB - 1024 / BITS_PER_MP_LIMB + 1; + r->_mp_size = x.s.sig == 0 ? size : -size; +} diff --git a/mpf/set_dfl_prec.c b/mpf/set_dfl_prec.c new file mode 100644 index 000000000..55069e4e4 --- /dev/null +++ b/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/mpf/set_prc.c b/mpf/set_prc.c new file mode 100644 index 000000000..10f2b0671 --- /dev/null +++ b/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/mpf/set_prc_raw.c b/mpf/set_prc_raw.c new file mode 100644 index 000000000..f55188a12 --- /dev/null +++ b/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/mpf/set_si.c b/mpf/set_si.c new file mode 100644 index 000000000..f9b4b55eb --- /dev/null +++ b/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/mpf/set_str.c b/mpf/set_str.c new file mode 100644 index 000000000..2ab9faebe --- /dev/null +++ b/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/mpf/set_ui.c b/mpf/set_ui.c new file mode 100644 index 000000000..ead0498aa --- /dev/null +++ b/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/mpf/size.c b/mpf/size.c new file mode 100644 index 000000000..23a57ec72 --- /dev/null +++ b/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/mpf/sqrt.c b/mpf/sqrt.c new file mode 100644 index 000000000..a0cc4ddb7 --- /dev/null +++ b/mpf/sqrt.c @@ -0,0 +1,79 @@ +/* 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. */ + +/* This code is just correct if "unsigned char" has at least 8 bits. It + doesn't help to use CHAR_BIT from limits.h, as the real problem is + the static arrays. */ + +#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 / usize > 0; /* Divide by zero for negative OP. */ + r->_mp_size = 0; + 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/mpf/sqrt_ui.c b/mpf/sqrt_ui.c new file mode 100644 index 000000000..987ef0776 --- /dev/null +++ b/mpf/sqrt_ui.c @@ -0,0 +1,65 @@ +/* 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. */ + +/* This code is just correct if "unsigned char" has at least 8 bits. It + doesn't help to use CHAR_BIT from limits.h, as the real problem is + the static arrays. */ + +#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/mpf/sub.c b/mpf/sub.c new file mode 100644 index 000000000..b87198263 --- /dev/null +++ b/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/mpf/sub_ui.c b/mpf/sub_ui.c new file mode 100644 index 000000000..780621378 --- /dev/null +++ b/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/mpf/tests/Makefile.in b/mpf/tests/Makefile.in new file mode 100644 index 000000000..62815eba3 --- /dev/null +++ b/mpf/tests/Makefile.in @@ -0,0 +1,58 @@ +# 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 ref.c +TEST_OBJS = t-add.o t-sub.o t-conv.o t-sqrt.o +TESTS = t-add t-sub t-conv t-sqrt + +check: st-add st-sub st-conv st-sqrt + @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 $@ + +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 ref.o $(TEST_LIBS) + $(CC) -o $@ $@.o ref.o $(TEST_LIBS) $(CFLAGS) +t-sqrt: t-sqrt.o ref.o $(TEST_LIBS) + $(CC) -o $@ $@.o ref.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-ref.o: $(srcdir)/t-ref.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/mpf/tests/configure.in b/mpf/tests/configure.in new file mode 100644 index 000000000..319219cd7 --- /dev/null +++ b/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/mpf/tests/ref.c b/mpf/tests/ref.c new file mode 100644 index 000000000..474c699ac --- /dev/null +++ b/mpf/tests/ref.c @@ -0,0 +1,184 @@ +/* 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" + +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); + +void +ref_mpf_add (mpf_t w, const mpf_t u, const mpf_t v) +{ + 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 = 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 = 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 = TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = 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 +ref_mpf_sub (mpf_t w, const mpf_t u, const mpf_t v) +{ + 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 = 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 = 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 = TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + vt = TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB); + wt = 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/mpf/tests/t-add.c b/mpf/tests/t-add.c new file mode 100644 index 000000000..93cf81a25 --- /dev/null +++ b/mpf/tests/t-add.c @@ -0,0 +1,109 @@ +/* Test mpf_neg, 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 + +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); + +main (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 (mpf_t x) +{ + mp_size_t i; + printf (" exp = %d\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08X ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/mpf/tests/t-conv.c b/mpf/tests/t-conv.c new file mode 100644 index 000000000..4467f0653 --- /dev/null +++ b/mpf/tests/t-conv.c @@ -0,0 +1,117 @@ +/* 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 (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 (mpf_t x) +{ + int i; + printf (" exp = %d\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08X ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/mpf/tests/t-sqrt.c b/mpf/tests/t-sqrt.c new file mode 100644 index 000000000..af13836ec --- /dev/null +++ b/mpf/tests/t-sqrt.c @@ -0,0 +1,100 @@ +/* Test mpf_sqrt, mpf_neg, 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 + +main (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 (mpf_t x) +{ + mp_size_t i; + printf (" exp = %d\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08X ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/mpf/tests/t-sub.c b/mpf/tests/t-sub.c new file mode 100644 index 000000000..6ff9bcb86 --- /dev/null +++ b/mpf/tests/t-sub.c @@ -0,0 +1,114 @@ +/* Test mpf_neg, 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 + +void ref_mpf_add (mpf_t, const mpf_t, const mpf_t); +void ref_mpf_sub (mpf_t, const mpf_t, const mpf_t); + +main (int argc, char **argv) +{ + mp_size_t size; + mp_exp_t exp; + int reps = 1000000; + 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 (mpf_t x) +{ + mp_size_t i; + printf (" exp = %d\n", x->_mp_exp); + printf ("size = %d\n", x->_mp_size); + for (i = ABS (x->_mp_size) - 1; i >= 0; i--) + printf ("%08X ", x->_mp_d[i]); + printf ("\n"); + mpf_dump (x); +} diff --git a/mpf/ui_div.c b/mpf/ui_div.c new file mode 100644 index 000000000..5ec4d1c9e --- /dev/null +++ b/mpf/ui_div.c @@ -0,0 +1,116 @@ +/* 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 dividend, mpf_srcptr v) +#else +mpf_ui_div (r, dividend, v) + mpf_ptr r; + unsigned long int dividend; + mpf_srcptr v; +#endif +{ + mp_ptr up, vp; + mp_ptr rp; + mp_size_t usize, vsize, rsize; + mp_size_t abs_vsize; + mp_size_t prec; + unsigned normalization_steps; + mp_limb_t qlimb; + mp_exp_t rexp; + TMP_DECL (marker); + + TMP_MARK (marker); + prec = r->_mp_prec; + rp = r->_mp_d; + + vp = v->_mp_d; + rexp = 1 - v->_mp_exp; + vsize = v->_mp_size; + abs_vsize = ABS (vsize); + + usize = abs_vsize + prec; + up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); + MPN_ZERO (up, usize); + + count_leading_zeros (normalization_steps, vp[abs_vsize - 1]); + + /* Normalize the denominator and the numerator. */ + if (normalization_steps != 0) + { + mp_ptr tp; + 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. */ + tp = (mp_ptr) TMP_ALLOC (abs_vsize * BYTES_PER_MP_LIMB); + mpn_lshift (tp, vp, abs_vsize, normalization_steps); + vp = tp; + + /* Shift up the dividend, possibly introducing a new most + significant word. */ + dividend_high = (mp_limb_t) dividend >> (BITS_PER_MP_LIMB - normalization_steps); + dividend_low = (mp_limb_t) dividend << normalization_steps; + + if (dividend_high != 0) + { + up[usize] = dividend_high; + up[usize - 1] = dividend_low; + rexp++; + } + else + { + usize--; + up[usize] = dividend_low; + } + } + else + { + /* The divisor is already normalized, as required. + Copy it to temporary space if it overlaps with the quotient. */ + if (vp == rp) + { + vp = (mp_ptr) TMP_ALLOC (abs_vsize * BYTES_PER_MP_LIMB); + MPN_COPY (vp, rp, abs_vsize); + } + + usize--; + up[usize] = dividend; + } + + qlimb = mpn_divmod (rp, up, usize + 1, vp, abs_vsize); + rsize = usize + 1 - abs_vsize; + if (qlimb) + { + rp[rsize] = qlimb; + rsize++; + rexp++; + } + r->_mp_size = vsize >= 0 ? rsize : -rsize; + r->_mp_exp = rexp; + TMP_FREE (marker); +} diff --git a/mpf/ui_sub.c b/mpf/ui_sub.c new file mode 100644 index 000000000..acb9210a6 --- /dev/null +++ b/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); +} |