[Gmp-commit] /var/hg/gmp: Use integer fields for mp_bases logarithm tables.
mercurial at gmplib.org
mercurial at gmplib.org
Sun Aug 7 21:28:04 CEST 2011
details: /var/hg/gmp/rev/88ab17a91350
changeset: 14222:88ab17a91350
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Aug 07 21:28:02 2011 +0200
description:
Use integer fields for mp_bases logarithm tables.
diffstat:
ChangeLog | 14 ++++++
gen-bases.c | 98 +++++++++++++++++++++++++++++++++++++++--------
gmp-impl.h | 47 ++++++++++++++++++++--
mpf/get_str.c | 11 ++--
mpf/out_str.c | 3 +-
mpf/set_str.c | 9 +--
mpn/generic/get_str.c | 7 ++-
mpn/generic/sizeinbase.c | 8 ++-
mpq/get_str.c | 10 ++--
mpz/inp_str.c | 8 +-
mpz/out_str.c | 8 ++-
mpz/set_str.c | 8 +-
printf/doprntf.c | 3 +-
tune/speed.h | 7 +-
tune/tuneup.c | 3 +-
15 files changed, 182 insertions(+), 62 deletions(-)
diffs (truncated from 591 to 300 lines):
diff -r 24aef6efb68d -r 88ab17a91350 ChangeLog
--- a/ChangeLog Thu Aug 04 19:17:06 2011 +0200
+++ b/ChangeLog Sun Aug 07 21:28:02 2011 +0200
@@ -1,3 +1,17 @@
+2011-08-07 Torbjorn Granlund <tege at gmplib.org>
+
+ * gmp-impl.h (struct bases): Add log2b and logb2 field, remove
+ chars_per_limb_exactly field.
+ (DIGITS_IN_BASE_FROM_BITS): New.
+ (DIGITS_IN_BASE_PER_LIMB): New.
+ (LIMBS_PER_DIGIT_IN_BASE): New.
+ * gen-bases.c: Generate log2b and logb2 fields; do not generate
+ chars_per_limb_exactly field.
+ * mpf/get_str.c mpf/out_str.c mpf/set_str.c mpn/generic/get_str.c
+ mpn/generic/sizeinbase.c mpq/get_str.c mpz/inp_str.c mpz/out_str.c
+ mpz/set_str.c printf/doprntf.c tune/speed.h tune/tuneup.c:
+ Use new macros.
+
2011-08-03 Torbjorn Granlund <tege at gmplib.org>
* dumbmp.c (mpz_tdiv_qr): Correctly handle dividend value being equal
diff -r 24aef6efb68d -r 88ab17a91350 gen-bases.c
--- a/gen-bases.c Thu Aug 04 19:17:06 2011 +0200
+++ b/gen-bases.c Sun Aug 07 21:28:02 2011 +0200
@@ -1,7 +1,7 @@
/* Generate mp_bases data.
-Copyright 1991, 1993, 1994, 1996, 2000, 2002, 2004 Free Software Foundation,
-Inc.
+Copyright 1991, 1993, 1994, 1996, 2000, 2002, 2004, 2011 Free Software
+Foundation, Inc.
This file is part of the GNU MP Library.
@@ -24,7 +24,6 @@
int chars_per_limb;
-double chars_per_bit_exactly;
mpz_t big_base;
int normalization_steps;
mpz_t big_base_inverted;
@@ -59,8 +58,6 @@
chars_per_limb++;
}
- chars_per_bit_exactly = 0.69314718055994530942 / log ((double) base);
-
mpz_ui_pow_ui (big_base, (long) base, (long) chars_per_limb);
normalization_steps = limb_bits - mpz_sizeinbase (big_base, 2);
@@ -97,11 +94,61 @@
printf ("#define MP_BASES_NORMALIZATION_STEPS_10 %d\n", normalization_steps);
}
+
+#define EXTRA 16
+
+/* Compute log(2)/log(b) as a fixnum. */
+void
+mp_2logb (mpz_t r, int bi, int prec)
+{
+ mpz_t t, t2, two, b;
+ int i;
+
+ mpz_init_set_ui (t, 1);
+ mpz_mul_2exp (t, t, prec+EXTRA);
+
+ mpz_init (t2);
+
+ mpz_init_set_ui (two, 2);
+ mpz_mul_2exp (two, two, prec+EXTRA);
+
+ mpz_set_ui (r, 0);
+
+ mpz_init_set_ui (b, bi);
+ mpz_mul_2exp (b, b, prec+EXTRA);
+
+ for (i = prec-1; i >= 0; i--)
+ {
+ mpz_mul_2exp (b, b, prec+EXTRA);
+ mpz_sqrt (b, b);
+
+ mpz_mul (t2, t, b);
+ mpz_tdiv_q_2exp (t2, t2, prec+EXTRA);
+
+ if (mpz_cmp (t2, two) < 0) /* not too large? */
+ {
+ mpz_setbit (r, i); /* set next less significant bit */
+ mpz_set (t, t2); /* new value acceptable */
+ }
+ }
+
+ mpz_clear (t);
+ mpz_clear (t2);
+ mpz_clear (two);
+ mpz_clear (b);
+}
+
void
table (int limb_bits, int nail_bits)
{
int numb_bits = limb_bits - nail_bits;
int base;
+ mpz_t r, t, logb2, log2b;
+
+ mpz_init (r);
+ mpz_init (t);
+ mpz_init (logb2);
+ mpz_init (log2b);
printf ("/* This file generated by gen-bases.c - DO NOT EDIT. */\n");
printf ("\n");
@@ -113,30 +160,45 @@
printf ("#endif\n");
printf ("\n");
puts ("const struct bases mp_bases[257] =\n{");
- puts (" /* 0 */ { 0, 0.0, 0 },");
- puts (" /* 1 */ { 0, 1e37, 0 },");
+ puts (" /* 0 */ { 0, 0.0, 0, 0, 0, 0 },");
+ puts (" /* 1 */ { 0, 1e37, 0, 0, 0, 0 },");
for (base = 2; base <= 256; base++)
{
generate (limb_bits, nail_bits, base);
+ mp_2logb (r, base, limb_bits + 8);
+ mpz_tdiv_q_2exp (logb2, r, 8);
+ mpz_set_ui (t, 1);
+ mpz_mul_2exp (t, t, 2*limb_bits + 5);
+ mpz_sub_ui (t, t, 1);
+ mpz_add_ui (r, r, 1);
+ mpz_tdiv_q (log2b, t, r);
printf (" /* %3u */ { ", base);
if (POW2_P (base))
{
- printf ("%u, %.16f, 0x%x },\n",
- chars_per_limb, chars_per_bit_exactly, ulog2 (base) - 1);
+ mpz_set_ui (big_base, ulog2 (base) - 1);
+ mpz_set_ui (big_base_inverted, 0);
}
- else
- {
- printf ("%u, %.16f, CNST_LIMB(0x",
- chars_per_limb, chars_per_bit_exactly);
- mpz_out_str (stdout, 16, big_base);
- printf ("), CNST_LIMB(0x");
- mpz_out_str (stdout, 16, big_base_inverted);
- printf (") },\n");
- }
+
+ printf ("%u,", chars_per_limb);
+ printf (" CNST_LIMB(0x");
+ mpz_out_str (stdout, 16, logb2);
+ printf ("), CNST_LIMB(0x");
+ mpz_out_str (stdout, 16, log2b);
+ printf ("), CNST_LIMB(0x");
+ mpz_out_str (stdout, 16, big_base);
+ printf ("), CNST_LIMB(0x");
+ mpz_out_str (stdout, 16, big_base_inverted);
+ printf (") },\n");
}
puts ("};");
+
+ mpz_clear (r);
+ mpz_clear (t);
+ mpz_clear (logb2);
+ mpz_clear (log2b);
+
}
int
diff -r 24aef6efb68d -r 88ab17a91350 gmp-impl.h
--- a/gmp-impl.h Thu Aug 04 19:17:06 2011 +0200
+++ b/gmp-impl.h Sun Aug 07 21:28:02 2011 +0200
@@ -2594,7 +2594,10 @@
int chars_per_limb;
/* log(2)/log(conversion_base) */
- double chars_per_bit_exactly;
+ mp_limb_t logb2;
+
+ /* log(conversion_base)/log(2) */
+ mp_limb_t log2b;
/* base**chars_per_limb, i.e. the biggest number that fits a word, built by
factors of base. Exception: For 2, 4, 8, etc, big_base is log2(base),
@@ -2611,6 +2614,20 @@
__GMP_DECLSPEC extern const struct bases mp_bases[257];
+/* Compute the number of digits in base for nbits bits, making sure the result
+ is never too small. The +1 in umul_ppmm makes the rounded-down log value to
+ be a rounded-up value. This make the full ph,,dummy product be rounded up,
+ but we chop the low half, meaning that we round down, but just once. The
+ return value is then incremented, to make sure it is never too small.
+ Caveat: The table value for base 2 will be all-bit-set, so things will break
+ in that case. Consider alternatives that allow all bases. */
+#define DIGITS_IN_BASE_FROM_BITS(res, nbits, b) \
+ do { \
+ mp_limb_t ph, dummy; \
+ umul_ppmm (ph, dummy, mp_bases[b].logb2 + 1, nbits); \
+ res = ph + 1; \
+ } while (0)
+
/* For power of 2 bases this is exact. For other bases the result is either
exact or one too big.
@@ -2644,8 +2661,9 @@
(result) = (__totbits + __lb_base - 1) / __lb_base; \
} \
else \
- (result) = (size_t) \
- (__totbits * mp_bases[base].chars_per_bit_exactly) + 1; \
+ { \
+ DIGITS_IN_BASE_FROM_BITS (result, __totbits, base); \
+ } \
} \
} while (0)
@@ -4031,6 +4049,24 @@
__GMP_DECLSPEC extern mp_size_t __gmp_default_fp_limb_precision;
+/* Compute the number of base-b digits corresponding to nlimbs limbs, rounding
+ down. */
+#define DIGITS_IN_BASE_PER_LIMB(res, nlimbs, b) \
+ do { \
+ mp_limb_t ph, dummy; \
+ umul_ppmm (ph, dummy, mp_bases[b].logb2, GMP_NUMB_BITS * (nlimbs)); \
+ res = ph; \
+ } while (0)
+
+/* Compute the number of limbs corresponding to ndigits base-b digits, rounding
+ up. */
+#define LIMBS_PER_DIGIT_IN_BASE(res, ndigits, b) \
+ do { \
+ mp_limb_t ph, dummy; \
+ umul_ppmm (ph, dummy, mp_bases[base].log2b, ndigits); \
+ res = 8 * ph / GMP_NUMB_BITS + 2; \
+ } while (0)
+
/* Set n to the number of significant digits an mpf of the given _mp_prec
field, in the given base. This is a rounded up value, designed to ensure
@@ -4046,9 +4082,10 @@
#define MPF_SIGNIFICANT_DIGITS(n, base, prec) \
do { \
+ size_t rawn; \
ASSERT (base >= 2 && base < numberof (mp_bases)); \
- (n) = 2 + (size_t) ((((size_t) (prec) - 1) * GMP_NUMB_BITS) \
- * mp_bases[(base)].chars_per_bit_exactly); \
+ DIGITS_IN_BASE_PER_LIMB (rawn, (prec) - 1, base); \
+ n = rawn + 2; \
} while (0)
diff -r 24aef6efb68d -r 88ab17a91350 mpf/get_str.c
--- a/mpf/get_str.c Thu Aug 04 19:17:06 2011 +0200
+++ b/mpf/get_str.c Sun Aug 07 21:28:02 2011 +0200
@@ -4,8 +4,8 @@
example, the number 3.1416 would be returned as "31416" in DIGIT_PTR and
1 in EXP.
-Copyright 1993, 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2005, 2006 Free
-Software Foundation, Inc.
+Copyright 1993, 1994, 1995, 1996, 1997, 2000, 2001, 2002, 2003, 2005, 2006, 2011
+Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -178,8 +178,7 @@
conversion.) */
tstr = (unsigned char *) TMP_ALLOC (n_digits + 2 * GMP_LIMB_BITS + 3);
- n_limbs_needed = 2 + (mp_size_t)
- (n_digits / (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly));
+ LIMBS_PER_DIGIT_IN_BASE (n_limbs_needed, n_digits, base);
if (ue <= n_limbs_needed)
{
@@ -188,7 +187,7 @@
unsigned long e;
n_more_limbs_needed = n_limbs_needed - ue;
- e = (unsigned long) n_more_limbs_needed * (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly);
+ DIGITS_IN_BASE_PER_LIMB (e, n_more_limbs_needed, base);
if (un > n_limbs_needed)
{
@@ -225,7 +224,7 @@
mp_ptr dummyp, xp;
n_less_limbs_needed = ue - n_limbs_needed;
- e = (unsigned long) n_less_limbs_needed * (GMP_NUMB_BITS * mp_bases[base].chars_per_bit_exactly);
+ DIGITS_IN_BASE_PER_LIMB (e, n_less_limbs_needed, base);
if (un > n_limbs_needed)
{
More information about the gmp-commit
mailing list