[Gmp-commit] /var/hg/gmp: 2 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Thu Feb 9 21:44:19 CET 2012
details: /var/hg/gmp/rev/9ae3e277d796
changeset: 14611:9ae3e277d796
user: Torbjorn Granlund <tege at gmplib.org>
date: Thu Feb 09 19:40:08 2012 +0100
description:
Fix off-by-one condition in invert_limb code.
details: /var/hg/gmp/rev/fb26dce883f5
changeset: 14612:fb26dce883f5
user: Torbjorn Granlund <tege at gmplib.org>
date: Thu Feb 09 21:44:14 2012 +0100
description:
Trivial merge.
diffstat:
ChangeLog | 22 +++++++++++++++++++++-
doc/gmp.texi | 23 +++++++++++++++++------
gmp-impl.h | 2 +-
mpn/powerpc32/divrem_2.asm | 10 +++++-----
mpz/fac_ui.c | 2 +-
mpz/oddfac_1.c | 38 +++++++++++++++++++++++++++-----------
tests/refmpn.c | 8 ++++----
7 files changed, 76 insertions(+), 29 deletions(-)
diffs (252 lines):
diff -r 44083fea3ee4 -r fb26dce883f5 ChangeLog
--- a/ChangeLog Fri Feb 03 15:51:26 2012 +0100
+++ b/ChangeLog Thu Feb 09 21:44:14 2012 +0100
@@ -1,3 +1,23 @@
+2012-02-09 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+ * mpz/oddfac_1.c (mpz_oddfac_1): Get ready for n!!
+ * gmp-impl.h (mpz_oddfac_1): Update signature.
+ * mpz/fac_ui.c (mpz_fac_ui): Update call to mpz_oddfac_1.
+
+2012-02-08 Torbjorn Granlund <tege at gmplib.org>
+
+ * mpn/powerpc32/divrem_2.asm: Fix off-by-one condition in invert_limb
+ code.
+
+2012-02-08 Niels Möller <nisse at lysator.liu.se>
+
+ * doc/gmp.texi (mpz_gcdext): Clarified corner cases in cofactor
+ canonicalization.
+
+2012-02-04 Marco Bodrato <bodrato at mail.dm.unipi.it>
+
+ * tests/refmpn.c (refmpn_mul): More conservative allocations.
+
2012-02-03 Torbjorn Granlund <tege at gmplib.org>
* longlong.h (udiv_qrnnd from sdiv_qrnnd): Declare udiv_w_sdiv.
@@ -836,7 +856,7 @@
* configure.in: Support s390x.
- * longlong.h: Add spport for 64-bit s390x.
+ * longlong.h: Add support for 64-bit s390x.
* mpn/s390_64: New directory.
* mpn/s390_64/add_n.asm: New file.
diff -r 44083fea3ee4 -r fb26dce883f5 doc/gmp.texi
--- a/doc/gmp.texi Fri Feb 03 15:51:26 2012 +0100
+++ b/doc/gmp.texi Thu Feb 09 21:44:14 2012 +0100
@@ -15,7 +15,7 @@
arithmetic library, version @value{VERSION}.
Copyright 1991, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
-2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software
+2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Free Software
Foundation, Inc.
Permission is granted to copy, distribute and/or modify this document under
@@ -3551,11 +3551,22 @@
addition set @var{s} and @var{t} to coefficients satisfying
@math{@var{a}@GMPmultiply{}@var{s} + @var{b}@GMPmultiply{}@var{t} = @var{g}}.
The value in @var{g} is always positive, even if one or both of @var{a} and
- at var{b} are negative (or zero if both inputs are zero). The values in @var{s} and @var{t} are chosen such that
-normally, @math{@GMPabs{@var{s}} < @GMPabs{@var{b}} / (2 @var{g})} and
- at math{@GMPabs{@var{t}} < @GMPabs{@var{a}} / (2 @var{g})}. The exceptional
-cases are that @math{@var{s} = sgn(@var{a})} if @math{@var{b} = 0} and
- at math{@var{t} = sgn(@var{b})} if @var{b = 0}.
+ at var{b} are negative (or zero if both inputs are zero). The values in @var{s}
+and @var{t} are chosen such that normally, @math{@GMPabs{@var{s}} <
+ at GMPabs{@var{b}} / (2 @var{g})} and @math{@GMPabs{@var{t}} < @GMPabs{@var{a}}
+/ (2 @var{g})}, and these relations define @var{s} and @var{t} uniquely. There
+are a few exceptional cases:
+
+If @math{@GMPabs{@var{a}} = @GMPabs{@var{b}}}, then @math{@var{s} = 0},
+ at math{@var{t} = sgn(@var{b})}.
+
+Otherwise, @math{@var{s} = sgn(@var{a})} if @math{@var{b} = 0} or
+ at math{@GMPabs{@var{b}} = 2 @var{g}}, and @math{@var{t} = sgn(@var{b})} if
+ at math{@var{a} = 0} or @math{@GMPabs{@var{a}} = 2 @var{g}}.
+
+In all cases, @math{@var{s} = 0} if and only if @math{@var{g} =
+ at GMPabs{@var{b}}}, i.e., if @var{b} divides @var{a} or @math{@var{a} = @var{b}
+= 0}.
If @var{t} is @code{NULL} then that value is not computed.
@end deftypefun
diff -r 44083fea3ee4 -r fb26dce883f5 gmp-impl.h
--- a/gmp-impl.h Fri Feb 03 15:51:26 2012 +0100
+++ b/gmp-impl.h Thu Feb 09 21:44:14 2012 +0100
@@ -1573,7 +1573,7 @@
__GMP_DECLSPEC mp_size_t mpz_prodlimbs __GMP_PROTO ((mpz_ptr, mp_ptr, mp_size_t));
#define mpz_oddfac_1 __gmpz_oddfac_1
-__GMP_DECLSPEC void mpz_oddfac_1 __GMP_PROTO ((mpz_ptr, mp_limb_t));
+__GMP_DECLSPEC void mpz_oddfac_1 __GMP_PROTO ((mpz_ptr, mp_limb_t, unsigned));
#define mpz_inp_str_nowhite __gmpz_inp_str_nowhite
#ifdef _GMP_H_HAVE_FILE
diff -r 44083fea3ee4 -r fb26dce883f5 mpn/powerpc32/divrem_2.asm
--- a/mpn/powerpc32/divrem_2.asm Fri Feb 03 15:51:26 2012 +0100
+++ b/mpn/powerpc32/divrem_2.asm Thu Feb 09 21:44:14 2012 +0100
@@ -1,6 +1,6 @@
dnl PPC-32 mpn_divrem_2 -- Divide an mpn number by a normalized 2-limb number.
-dnl Copyright 2007, 2008 Free Software Foundation, Inc.
+dnl Copyright 2007, 2008, 2012 Free Software Foundation, Inc.
dnl This file is part of the GNU MP Library.
@@ -84,9 +84,9 @@
bge- cr7, L(9)
add r0, r0, r10
cmplw cr7, r0, r10
- cmplw cr6, r0, r6
+ cmplw cr6, r6, r0
addi r31, r31, -1 C q1--
- cror 28, 28, 25
+ crorc 28, 28, 25
bc+ 12, 28, L(9)
addi r31, r31, -1 C q1--
add r0, r0, r10
@@ -101,9 +101,9 @@
bge- cr7, L(13)
add r0, r0, r10
cmplw cr7, r0, r10
- cmplw cr6, r0, r11
+ cmplw cr6, r11, r0
addi r6, r6, -1 C q0--
- cror 28, 28, 25
+ crorc 28, 28, 25
bc+ 12, 28, L(13)
C add r0, r0, r10 C final remainder
addi r6, r6, -1 C q0--
diff -r 44083fea3ee4 -r fb26dce883f5 mpz/fac_ui.c
--- a/mpz/fac_ui.c Fri Feb 03 15:51:26 2012 +0100
+++ b/mpz/fac_ui.c Thu Feb 09 21:44:14 2012 +0100
@@ -92,7 +92,7 @@
else
{
mp_limb_t count;
- mpz_oddfac_1 (x, n);
+ mpz_oddfac_1 (x, n, 0);
popc_limb (count, n);
mpz_mul_2exp (x, x, n - count);
}
diff -r 44083fea3ee4 -r fb26dce883f5 mpz/oddfac_1.c
--- a/mpz/oddfac_1.c Fri Feb 03 15:51:26 2012 +0100
+++ b/mpz/oddfac_1.c Thu Feb 09 21:44:14 2012 +0100
@@ -463,12 +463,18 @@
/* mpz_oddfac_1 computes the odd part of the factorial of the
parameter n. I.e. n! = x 2^a, where x is the returned value: an
odd positive integer.
+
+ If flag != 0 a square is skipped in the DSC part, e.g.
+ if n is odd, n > FAC_DSC_THRESHOLD and flag = 1, x is set to n!!.
+
+ If n is too small, flag is ignored, and an ASSERT can be triggered.
*/
void
-mpz_oddfac_1 (mpz_ptr x, mp_limb_t n)
+mpz_oddfac_1 (mpz_ptr x, mp_limb_t n, unsigned flag)
{
static const mp_limb_t tablef[] = { ONE_LIMB_ODD_FACTORIAL_TABLE };
- static const mp_limb_t tabled[] = { ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE };
+
+ ASSERT (flag == 0 || (n >= numberof (tablef) && ABOVE_THRESHOLD (n, FAC_DSC_THRESHOLD)));
if (n < numberof (tablef))
{
@@ -484,6 +490,7 @@
s = 0;
{
+ static const mp_limb_t tabled[] = { ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE };
mp_limb_t tn;
mp_size_t j;
TMP_SDECL;
@@ -546,13 +553,15 @@
TMP_MARK;
- size = n / GMP_NUMB_BITS;
- ASSERT (primesieve_size (n - 1) <= (size + 4) - (size / 2 + 3));
+ flag--;
+ size = n / GMP_NUMB_BITS + 4;
+ ASSERT (primesieve_size (n - 1) <= size - (size / 2 + 1));
/* 2-multiswing(n) < 2^(n-1)*sqrt(n/pi) < 2^(n+GMP_NUMB_BITS);
one more can be overwritten by mul, another for the sieve */
- MPZ_TMP_INIT (mswing, size + 4);
+ MPZ_TMP_INIT (mswing, size);
/* Put the sieve on the second half, it will be overwritten by the last mswing. */
- sieve = PTR(mswing) + size / 2 + 3;
+ sieve = PTR (mswing) + size / 2 + 1;
+ ASSERT ((SIZ (mswing) = 0) || ALLOC (mswing) == size);
size = (bitwise_primesieve (sieve, n - 1) + 1) / log_n_max (n) + 1;
@@ -564,15 +573,22 @@
TMP_DECL;
s--;
+ ASSERT (ABSIZ (mswing) < ALLOC (mswing) / 2); /* Check: sieve has not been overwritten */
mpz_2multiswing_1 (mswing, n >> s, sieve, factors);
TMP_MARK;
nx = SIZ (x);
- size = nx << 1;
- square = TMP_ALLOC_LIMBS (size);
- mpn_sqr (square, PTR (x), nx);
- ns = SIZ(mswing);
- size -= (square[size - 1] == 0);
+ if (s == flag) {
+ size = nx;
+ square = TMP_ALLOC_LIMBS (size);
+ MPN_COPY (square, PTR (x), nx);
+ } else {
+ size = nx << 1;
+ square = TMP_ALLOC_LIMBS (size);
+ mpn_sqr (square, PTR (x), nx);
+ size -= (square[size - 1] == 0);
+ }
+ ns = SIZ (mswing);
nx = size + ns;
MPZ_REALLOC (x, nx);
ASSERT (ns <= size);
diff -r 44083fea3ee4 -r fb26dce883f5 tests/refmpn.c
--- a/tests/refmpn.c Fri Feb 03 15:51:26 2012 +0100
+++ b/tests/refmpn.c Thu Feb 09 21:44:14 2012 +0100
@@ -2,7 +2,7 @@
of the normal gmp code. Speed isn't a consideration.
Copyright 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006,
-2007, 2008, 2009, 2011 Free Software Foundation, Inc.
+2007, 2008, 2009, 2011, 2012 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -1881,21 +1881,21 @@
if (vn < TOOM4_THRESHOLD)
{
/* In the mpn_toom33_mul range, use mpn_toom22_mul. */
- tn = 2 * vn + mpn_toom22_mul_itch (vn, vn);
+ tn = 3 * vn + mpn_toom22_mul_itch (vn, vn);
tp = refmpn_malloc_limbs (tn);
mpn_toom22_mul (tp, up, vn, vp, vn, tp + 2 * vn);
}
else if (vn < FFT_THRESHOLD)
{
/* In the mpn_toom44_mul range, use mpn_toom33_mul. */
- tn = 2 * vn + mpn_toom33_mul_itch (vn, vn);
+ tn = 3 * vn + mpn_toom33_mul_itch (vn, vn);
tp = refmpn_malloc_limbs (tn);
mpn_toom33_mul (tp, up, vn, vp, vn, tp + 2 * vn);
}
else
{
/* Finally, for the largest operands, use mpn_toom44_mul. */
- tn = 2 * vn + mpn_toom44_mul_itch (vn, vn);
+ tn = 3 * vn + mpn_toom44_mul_itch (vn, vn);
tp = refmpn_malloc_limbs (tn);
mpn_toom44_mul (tp, up, vn, vp, vn, tp + 2 * vn);
}
More information about the gmp-commit
mailing list