Permalink
Browse files

Merge branch 'trunk' of https://github.com/fredrik-johansson/flint2 i…

…nto trunk
  • Loading branch information...
2 parents e97c9ec + c971aaf commit e59f1b9245a1652ece59a14e5610c9f53335fe50 @TomMD TomMD committed Apr 11, 2012
Showing with 461 additions and 0 deletions.
  1. +5 −0 fmpz.h
  2. +87 −0 fmpz/abs_lbound_ui_2exp.c
  3. +112 −0 fmpz/abs_ubound_ui_2exp.c
  4. +16 −0 fmpz/doc/fmpz.txt
  5. +102 −0 fmpz/test/t-abs_lbound_ui_2exp.c
  6. +139 −0 fmpz/test/t-abs_ubound_ui_2exp.c
View
5 fmpz.h
@@ -561,6 +561,11 @@ static __inline__ void fmpz_set_ui_smod(fmpz_t f, mp_limb_t x, mp_limb_t m)
fmpz_set_si(f, x - m);
}
+mp_limb_t fmpz_abs_ubound_ui_2exp(long * exp, const fmpz_t x, int bits);
+
+mp_limb_t fmpz_abs_lbound_ui_2exp(long * exp, const fmpz_t x, int bits);
+
+
#ifdef __cplusplus
}
#endif
View
@@ -0,0 +1,87 @@
+/*============================================================================
+
+ This file is part of FLINT.
+
+ FLINT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ FLINT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include <mpir.h>
+#include "flint.h"
+#include "ulong_extras.h"
+#include "fmpz.h"
+
+mp_limb_t
+fmpz_abs_lbound_ui_2exp(long * exp, const fmpz_t x, int bits)
+{
+ mp_limb_t m;
+ long shift, e, size;
+ fmpz c = *x;
+
+ if (!COEFF_IS_MPZ(c))
+ {
+ m = FLINT_ABS(c);
+ e = 0;
+ }
+ else
+ {
+ __mpz_struct * z = COEFF_TO_PTR(c);
+ size = z->_mp_size;
+ size = FLINT_ABS(size);
+ e = (size - 1) * FLINT_BITS;
+
+ if (size == 1)
+ {
+ m = z->_mp_d[0];
+ }
+ else /* there are two or more limbs */
+ {
+ /* top limb (which must be nonzero) */
+ m = z->_mp_d[size - 1];
+
+ count_leading_zeros(shift, m);
+ shift = FLINT_BITS - shift - bits;
+ e += shift;
+
+ if (shift >= 0)
+ {
+ m >>= shift;
+ }
+ else
+ {
+ /* read a second limb to get an accurate value */
+ mp_limb_t m2 = z->_mp_d[size - 2];
+ m = (m << (-shift)) | (m2 >> (FLINT_BITS + shift));
+ }
+
+ *exp = e;
+ return m;
+ }
+ }
+
+ count_leading_zeros(shift, m);
+ e += FLINT_BITS - shift - bits;
+ if (e >= 0)
+ m >>= e;
+ else
+ m <<= (-e);
+ *exp = e;
+ return m;
+}
View
@@ -0,0 +1,112 @@
+/*============================================================================
+
+ This file is part of FLINT.
+
+ FLINT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ FLINT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include <mpir.h>
+#include "flint.h"
+#include "ulong_extras.h"
+#include "fmpz.h"
+
+mp_limb_t
+fmpz_abs_ubound_ui_2exp(long * exp, const fmpz_t x, int bits)
+{
+ mp_limb_t m;
+ long shift, e, size;
+ fmpz c = *x;
+
+ if (!COEFF_IS_MPZ(c))
+ {
+ m = FLINT_ABS(c);
+ e = 0;
+ }
+ else
+ {
+ /* mpz */
+ __mpz_struct * z = COEFF_TO_PTR(c);
+ size = z->_mp_size;
+ size = FLINT_ABS(size);
+ e = (size - 1) * FLINT_BITS;
+
+ if (size == 1)
+ {
+ m = z->_mp_d[0];
+ }
+ else /* there are two or more limbs */
+ {
+ /* top limb (which must be nonzero) */
+ m = z->_mp_d[size - 1];
+
+ count_leading_zeros(shift, m);
+ shift = FLINT_BITS - shift - bits;
+ e += shift;
+
+ if (shift >= 0)
+ {
+ /* round up */
+ m = (m >> shift) + 1;
+ }
+ else
+ {
+ /* read a second limb to get an accurate value */
+ mp_limb_t m2 = z->_mp_d[size - 2];
+ m = (m << (-shift)) | (m2 >> (FLINT_BITS + shift));
+ /* round up */
+ m++;
+ }
+
+ /* adding 1 caused overflow to the next power of two */
+ if ((m & (m - 1UL)) == 0UL)
+ {
+ m = 1UL << (bits - 1);
+ e++;
+ }
+
+ *exp = e;
+ return m;
+ }
+ }
+
+ /* single limb, adjust */
+ count_leading_zeros(shift, m);
+ e = FLINT_BITS - shift - bits;
+
+ if (e >= 0)
+ {
+ m = (m >> e) + 1;
+
+ /* overflowed to next power of two */
+ if ((m & (m - 1)) == 0UL)
+ {
+ m = 1UL << (bits - 1);
+ e++;
+ }
+ }
+ else
+ {
+ m <<= (-e);
+ }
+
+ *exp = e;
+ return m;
+}
View
@@ -377,6 +377,22 @@ int fmpz_tstbit(const fmpz_t f, ulong i)
Test bit index~$i$ of $f$ and return $0$ or $1$, accordingly.
+mp_limb_t fmpz_abs_lbound_ui_2exp(long * exp, const fmpz_t x, int bits)
+
+ For nonzero $x$, returns a mantissa $m$ with exactly \code{bits} bits and
+ sets \code{exp} to an exponent $e$, such that $|x| \ge m 2^e$. The number
+ of bits must be between 1 and \code{FLINT_BITS} inclusive.
+ The mantissa is guaranteed to be correctly rounded.
+
+mp_limb_t fmpz_abs_ubound_ui_2exp(long * exp, const fmpz_t x, int bits)
+
+ For nonzero $x$, returns a mantissa $m$ with exactly \code{bits} bits
+ and sets \code{exp} to an exponent $e$, such that $|x| \le m 2^e$.
+ The number of bits must be between 1 and \code{FLINT_BITS} inclusive.
+ The mantissa is either correctly rounded or one unit too large
+ (possibly meaning that the exponent is one too large,
+ if the mantissa is a power of two).
+
*******************************************************************************
Comparison
@@ -0,0 +1,102 @@
+/*=============================================================================
+
+ This file is part of FLINT.
+
+ FLINT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ FLINT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with FLINT; if not, write to the Free Software
+ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+ Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <mpir.h>
+#include "flint.h"
+#include "ulong_extras.h"
+#include "fmpz.h"
+
+static mp_limb_t
+refimpl(long * exp, const fmpz_t x, int bits)
+{
+ fmpz_t t;
+ long xbits;
+ mp_limb_t m;
+
+ xbits = fmpz_bits(x);
+
+ fmpz_init(t);
+ fmpz_abs(t, x);
+
+ if (xbits >= bits)
+ fmpz_tdiv_q_2exp(t, t, xbits - bits);
+ else
+ fmpz_mul_2exp(t, t, bits - xbits);
+
+ m = fmpz_get_ui(t);
+ fmpz_clear(t);
+
+ *exp = xbits - bits;
+
+ return m;
+}
+
+int
+main(void)
+{
+ long iter;
+ flint_rand_t state;
+
+ printf("abs_lbound_ui_2exp....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 100000; iter++)
+ {
+ fmpz_t x;
+ long bits;
+ long exp, yexp;
+ mp_limb_t yman, man;
+
+ fmpz_init(x);
+ fmpz_randtest_not_zero(x, state, 1 + n_randint(state, 400));
+
+ bits = 1 + n_randint(state, FLINT_BITS - 1);
+
+ yman = refimpl(&yexp, x, bits);
+ man = fmpz_abs_lbound_ui_2exp(&exp, x, bits);
+
+ if (FLINT_BIT_COUNT(man) != bits || (man != yman) || (exp != yexp))
+ {
+ printf("FAIL\n");
+ printf("bits = %ld, count = %u\n\n", bits, FLINT_BIT_COUNT(man));
+ printf("x = "); fmpz_print(x); printf("\n\n");
+ printf("bits(x) = %ld\n\n", fmpz_bits(x));
+ printf("man = %lu, exp = %ld\n", man, exp);
+ printf("yman = %lu, yexp = %ld\n", yman, yexp);
+ abort();
+ }
+
+ fmpz_clear(x);
+ }
+
+ flint_randclear(state);
+ _fmpz_cleanup();
+ printf("PASS\n");
+ return 0;
+}
Oops, something went wrong.

0 comments on commit e59f1b9

Please sign in to comment.