diff options
| -rw-r--r-- | mpn/generic/get_str.c | 225 | ||||
| -rw-r--r-- | mpn/generic/set_str.c | 506 |
2 files changed, 391 insertions, 340 deletions
diff --git a/mpn/generic/get_str.c b/mpn/generic/get_str.c index 8ea256ac3..c52e2fffc 100644 --- a/mpn/generic/get_str.c +++ b/mpn/generic/get_str.c @@ -1,32 +1,32 @@ /* mpn_get_str -- Convert a MSIZE long limb vector pointed to by MPTR to a printable string in STR in base BASE. -Copyright 1991, 1992, 1993, 1994, 1996, 2000, 2001, 2002 Free Software -Foundation, Inc. +Copyright 1991, 1992, 1993, 1994, 1996, 2000, 2001, 2002, 2004, 2006 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 Lesser General Public License as published by -the Free Software Foundation; either version 2.1 of the License, or (at your +The GNU MP Library is free software; you can redistribute it and/or modify it +under the terms of the GNU Lesser General Public License as published by the +Free Software Foundation; either version 2.1 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 Lesser General Public -License for more details. +WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License +for more details. -You should have received a copy of the GNU Lesser General Public License -along with the GNU MP Library; see the file COPYING.LIB. If not, write to -the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, -MA 02110-1301, USA. */ +You should have received a copy of the GNU Lesser General Public License along +with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free +Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, +USA. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" /* Conversion of U {up,un} to a string in base b. Internally, we convert to - base B = b^m, the largest power of b that fits a limb. Basic algorithms: + base B = b^m, the largest power of b that fits a limb. Basic algorithms: A) Divide U repeatedly by B, generating a quotient and remainder, until the quotient becomes zero. The remainders hold the converted digits. Digits @@ -37,11 +37,10 @@ MA 02110-1301, USA. */ come out from left to right. (Currently not used herein, except for in code for converting single limbs to individual digits.) - C) Compute B^1, B^2, B^4, ..., B^(2^s), for s such that B^(2^s) > sqrt(U). - Then divide U by B^(2^k), generating an integer quotient and remainder. - Recursively convert the quotient, then the remainder, using the - precomputed powers. Digits come out from left to right. (Used in - mpn_dc_get_str.) + C) Compute B^1, B^2, B^4, ..., B^s, for s such that B^s just above sqrt(U). + Then divide U by B^s, generating quotient and remainder. Recursively + convert the quotient, then the remainder, using the precomputed powers. + Digits come out from left to right. (Used in mpn_dc_get_str.) When using algorithm C, algorithm B might be suitable for basecase code, since the required b^g power will be readily accessible. @@ -135,37 +134,26 @@ MA 02110-1301, USA. */ #define GET_STR_PRECOMPUTE_THRESHOLD 30 #endif -struct powers -{ - size_t digits_in_base; - mp_ptr p; - mp_size_t n; /* mpz_struct uses int for sizes, but not mpn! */ - int base; -}; -typedef struct powers powers_t; - -/* Convert {UP,UN} to a string with a base as represented in POWTAB, and put - the string in STR. Generate LEN characters, possibly padding with zeros to - the left. If LEN is zero, generate as many characters as required. - Return a pointer immediately after the last digit of the result string. - Complexity is O(UN^2); intended for small conversions. */ +/* Convert {up,un} to a string in base base, and put the result in str. + Generate len characters, possibly padding with zeros to the left. If len is + zero, generate as many characters as required. Return a pointer immediately + after the last digit of the result string. Complexity is O(un^2); intended + for small conversions. */ static unsigned char * mpn_sb_get_str (unsigned char *str, size_t len, - mp_ptr up, mp_size_t un, - const powers_t *powtab) + mp_ptr up, mp_size_t un, int base) { mp_limb_t rl, ul; unsigned char *s; - int base; size_t l; /* Allocate memory for largest possible string, given that we only get here for operands with un < GET_STR_PRECOMPUTE_THRESHOLD and that the smallest base is 3. 7/11 is an approximation to 1/log2(3). */ #if TUNE_PROGRAM_BUILD -#define BUF_ALLOC (GET_STR_THRESHOLD_LIMIT * BITS_PER_MP_LIMB * 7 / 11) +#define BUF_ALLOC (GET_STR_THRESHOLD_LIMIT * GMP_LIMB_BITS * 7 / 11) #else -#define BUF_ALLOC (GET_STR_PRECOMPUTE_THRESHOLD * BITS_PER_MP_LIMB * 7 / 11) +#define BUF_ALLOC (GET_STR_PRECOMPUTE_THRESHOLD * GMP_LIMB_BITS * 7 / 11) #endif unsigned char buf[BUF_ALLOC]; #if TUNE_PROGRAM_BUILD @@ -174,7 +162,6 @@ mpn_sb_get_str (unsigned char *str, size_t len, mp_limb_t rp[GET_STR_PRECOMPUTE_THRESHOLD]; #endif - base = powtab->base; if (base == 10) { /* Special case code for base==10 so that the compiler has a chance to @@ -194,10 +181,10 @@ mpn_sb_get_str (unsigned char *str, size_t len, un -= rp[un] == 0; frac = (rp[0] + 1) << GMP_NAIL_BITS; s -= MP_BASES_CHARS_PER_LIMB_10; - i = MP_BASES_CHARS_PER_LIMB_10; #if HAVE_HOST_CPU_FAMILY_x86 /* The code below turns out to be a bit slower for x86 using gcc. Use plain code. */ + i = MP_BASES_CHARS_PER_LIMB_10; do { umul_ppmm (digit, frac, frac, 10); @@ -212,31 +199,28 @@ mpn_sb_get_str (unsigned char *str, size_t len, { umul_ppmm (digit, frac, frac, 10); *s++ = digit; - i--; } if (MP_BASES_NORMALIZATION_STEPS_10 <= 1) { umul_ppmm (digit, frac, frac, 10); *s++ = digit; - i--; } if (MP_BASES_NORMALIZATION_STEPS_10 <= 2) { umul_ppmm (digit, frac, frac, 10); *s++ = digit; - i--; } if (MP_BASES_NORMALIZATION_STEPS_10 <= 3) { umul_ppmm (digit, frac, frac, 10); *s++ = digit; - i--; } + i = MP_BASES_CHARS_PER_LIMB_10 - (4-MP_BASES_NORMALIZATION_STEPS_10); frac = (frac + 0xf) >> 4; do { frac *= 10; - digit = frac >> (BITS_PER_MP_LIMB - 4); + digit = frac >> (GMP_LIMB_BITS - 4); *s++ = digit; frac &= (~(mp_limb_t) 0) >> 4; } @@ -323,7 +307,7 @@ mpn_dc_get_str (unsigned char *str, size_t len, if (un < GET_STR_DC_THRESHOLD) { if (un != 0) - str = mpn_sb_get_str (str, len, up, un, powtab); + str = mpn_sb_get_str (str, len, up, un, powtab->base); else { while (len != 0) @@ -352,8 +336,12 @@ mpn_dc_get_str (unsigned char *str, size_t len, mpn_tdiv_qr (qp, rp, 0L, up, un, pwp, pwn); qn = un - pwn; qn += qp[qn] != 0; /* quotient size */ + + ASSERT (qn < pwn || (qn == pwn && mpn_cmp (qp, pwp, pwn) < 0)); + if (len != 0) len = len - powtab->digits_in_base; + str = mpn_dc_get_str (str, len, qp, qn, powtab - 1, tmp + un - pwn + 1); str = mpn_dc_get_str (str, powtab->digits_in_base, rp, pwn, powtab - 1, tmp); } @@ -362,6 +350,27 @@ mpn_dc_get_str (unsigned char *str, size_t len, } +/* + Example sequences: + + 1 1 2 3 5 9 17 (worst case for un = 33,34) + 1 1 1 2s + 2a 4s + 3a 8s + 5a 16s + 9a + 17a + + 1 1 2 4 8 16 32 (best case for un = 63,64) + 1 1 2s 4s 8s 16s 32s + + 1 1 2 4 7 14 28 (medium) + 1 1 1 3sa + 2a 7sa + 4a 17s + 28s +*/ + /* There are no leading zeros on the digits generated at str, but that's not currently a documented feature. */ @@ -371,12 +380,13 @@ mpn_get_str (unsigned char *str, int base, mp_ptr up, mp_size_t un) mp_ptr powtab_mem, powtab_mem_ptr; mp_limb_t big_base; size_t digits_in_base; - powers_t powtab[30]; + powers_t powtab[GMP_LIMB_BITS]; int pi; mp_size_t n; mp_ptr p, t; size_t out_len; mp_ptr tmp; + TMP_DECL; /* Special case zero, as the code below doesn't handle it. */ if (un == 0) @@ -437,18 +447,16 @@ mpn_get_str (unsigned char *str, int base, mp_ptr up, mp_size_t un) /* General case. The base is not a power of 2. */ if (un < GET_STR_PRECOMPUTE_THRESHOLD) - { - struct powers ptab[1]; - ptab[0].base = base; - return mpn_sb_get_str (str, (size_t) 0, up, un, ptab) - str; - } + return mpn_sb_get_str (str, (size_t) 0, up, un, base) - str; + + TMP_MARK; /* Allocate one large block for the powers of big_base. With the current scheme, we need to allocate twice as much as would be possible if a minimal set of powers were generated. */ -#define POWTAB_ALLOC_SIZE (2 * un + 30) -#define TMP_ALLOC_SIZE (un + 30) - powtab_mem = __GMP_ALLOCATE_FUNC_LIMBS (POWTAB_ALLOC_SIZE); +#define POWTAB_ALLOC_SIZE (un + 2 * GMP_LIMB_BITS) +#define TMP_ALLOC_SIZE (un) + powtab_mem = TMP_BALLOC_LIMBS (POWTAB_ALLOC_SIZE); powtab_mem_ptr = powtab_mem; /* Compute a table of powers: big_base^1, big_base^2, big_base^4, ..., @@ -458,43 +466,84 @@ mpn_get_str (unsigned char *str, int base, mp_ptr up, mp_size_t un) big_base = __mp_bases[base].big_base; digits_in_base = __mp_bases[base].chars_per_limb; - powtab[0].base = base; /* FIXME: hack for getting base to mpn_sb_get_str */ - powtab[1].p = &big_base; - powtab[1].n = 1; - powtab[1].digits_in_base = digits_in_base; - powtab[1].base = base; - powtab[2].p = &big_base; - powtab[2].n = 1; - powtab[2].digits_in_base = digits_in_base; - powtab[2].base = base; - n = 1; - pi = 2; - p = &big_base; - for (;;) - { - ++pi; - t = powtab_mem_ptr; - powtab_mem_ptr += 2 * n; - mpn_sqr_n (t, p, n); - n *= 2; n -= t[n - 1] == 0; - digits_in_base *= 2; - p = t; - powtab[pi].p = p; - powtab[pi].n = n; - powtab[pi].digits_in_base = digits_in_base; - powtab[pi].base = base; - - if (2 * n > un) - break; + { + mp_size_t n_pows, xn, pn, exptab[GMP_LIMB_BITS], bexp; + mp_limb_t cy; + n_pows = 0; + xn = 1 + un*(__mp_bases[base].chars_per_bit_exactly*GMP_NUMB_BITS)/__mp_bases[base].chars_per_limb; + for (pn = xn; pn != 1; pn = (pn + 1) >> 1) + { + exptab[n_pows] = pn; + n_pows++; + } + exptab[n_pows] = 1; + + powtab[0].p = &big_base; + powtab[0].n = 1; + powtab[0].digits_in_base = digits_in_base; + powtab[0].base = base; + + powtab[1].p = powtab_mem_ptr; powtab_mem_ptr += 2; + powtab[1].p[0] = big_base; + powtab[1].n = 1; + powtab[1].digits_in_base = digits_in_base; + powtab[1].base = base; + + n = 1; + p = &big_base; + bexp = 1; + for (pi = 2; pi < n_pows; pi++) + { + t = powtab_mem_ptr; + powtab_mem_ptr += 2 * n + 2; + + ASSERT_ALWAYS (powtab_mem_ptr < powtab_mem + POWTAB_ALLOC_SIZE); + + mpn_sqr_n (t, p, n); + + digits_in_base *= 2; + n *= 2; n -= t[n - 1] == 0; + bexp *= 2; + + if (bexp + 1 < exptab[n_pows - pi]) + { + digits_in_base += __mp_bases[base].chars_per_limb; + cy = mpn_mul_1 (t, t, n, big_base); + t[n] = cy; + n += cy != 0; + bexp += 1; + } + p = t; + powtab[pi].p = p; + powtab[pi].n = n; + powtab[pi].digits_in_base = digits_in_base; + powtab[pi].base = base; + } + + for (pi = 1; pi < n_pows; pi++) + { + t = powtab[pi].p; + n = powtab[pi].n; + cy = mpn_mul_1 (t, t, n, big_base); + t[n] = cy; + n += cy != 0; + powtab[pi].n = n; + powtab[pi].digits_in_base += __mp_bases[base].chars_per_limb; + } + +#if 0 + { int i; + printf ("Computed table values for base=%d, un=%d, xn=%d:\n", base, un, xn); + for (i = 0; i < n_pows; i++) + printf ("%2d: %10ld %10ld %11ld\n", i, exptab[n_pows-i], powtab[i].n, powtab[i].digits_in_base); } - ASSERT_ALWAYS (POWTAB_ALLOC_SIZE > powtab_mem_ptr - powtab_mem); +#endif + } /* Using our precomputed powers, now in powtab[], convert our number. */ - tmp = __GMP_ALLOCATE_FUNC_LIMBS (TMP_ALLOC_SIZE); - out_len = mpn_dc_get_str (str, 0, up, un, powtab + pi, tmp) - str; - __GMP_FREE_FUNC_LIMBS (tmp, TMP_ALLOC_SIZE); - - __GMP_FREE_FUNC_LIMBS (powtab_mem, POWTAB_ALLOC_SIZE); + tmp = TMP_BALLOC_LIMBS (TMP_ALLOC_SIZE); + out_len = mpn_dc_get_str (str, 0, up, un, powtab - 1 + pi, tmp) - str; + TMP_FREE; return out_len; } diff --git a/mpn/generic/set_str.c b/mpn/generic/set_str.c index 4166bb1d7..f08965237 100644 --- a/mpn/generic/set_str.c +++ b/mpn/generic/set_str.c @@ -2,72 +2,76 @@ Convert a STR_LEN long base BASE byte string pointed to by STR to a limb vector pointed to by RES_PTR. Return the number of limbs in RES_PTR. -Copyright 1991, 1992, 1993, 1994, 1996, 2000, 2001, 2002 Free Software -Foundation, Inc. +Copyright 1991, 1992, 1993, 1994, 1996, 2000, 2001, 2002, 2004, 2006 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. +under the terms of the GNU Lesser General Public License as published by the +Free Software Foundation; either version 2.1 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 +FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. -You should have received a copy of the GNU Library General Public License along +You should have received a copy of the GNU Lesser General Public License along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free -Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ +/* Originally written 2004, but it was left non-working. + Fixed 2006-03 at LORIA. + + TODO: + + Perhaps do not compute the highest power? + Instead, multiply twice by the 2nd highest power: + + _______ + |_______| hp + |_______| pow + _______________ + |_______________| final result + + + _______ + |_______| hp + |___| pow[-1] + ___________ + |___________| intermediate result + |___| pow[-1] + _______________ + |_______________| final result + + Generalizing that idea, perhaps we should make powtab contain successive + cubes, not squares. +*/ + #include "gmp.h" #include "gmp-impl.h" +#include "longlong.h" -static mp_size_t convert_blocks __GMP_PROTO ((mp_ptr, const unsigned char *, size_t, int)); -/* When to switch to sub-quadratic code. This counts characters/digits in - the input string, not limbs as most other *_THRESHOLD. */ -#ifndef SET_STR_THRESHOLD -#define SET_STR_THRESHOLD 4000 +#ifndef SET_STR_DC_THRESHOLD +#define SET_STR_DC_THRESHOLD 2000 #endif -/* Don't define this to anything but 1 for now. In order to make other values - work well, either the convert_blocks function should be generalized to - handle larger blocks than chars_per_limb, or the basecase code should be - broken out of the main function. Also note that this must be a power of - 2. */ -#ifndef SET_STR_BLOCK_SIZE -#define SET_STR_BLOCK_SIZE 1 /* Must be a power of 2. */ +#ifndef SET_STR_PRECOMPUTE_THRESHOLD +#define SET_STR_PRECOMPUTE_THRESHOLD 5000 #endif +mp_size_t mpn_bc_set_str (mp_ptr, const unsigned char *, size_t, int); +mp_size_t mpn_dc_set_str (mp_ptr, const unsigned char *, size_t, + const powers_t *, mp_ptr); +void mpn_set_str_compute_powtab (powers_t *, mp_ptr, mp_size_t, int); -/* This check interferes with expression based values of SET_STR_THRESHOLD - used for tuning and measuring. -#if SET_STR_BLOCK_SIZE >= SET_STR_THRESHOLD -These values are silly. -The sub-quadratic code would recurse to itself. -#endif -*/ mp_size_t mpn_set_str (mp_ptr rp, const unsigned char *str, size_t str_len, int base) { - mp_size_t size; - mp_limb_t big_base; - int chars_per_limb; - mp_limb_t res_digit; - - ASSERT (base >= 2); - ASSERT (base < numberof (__mp_bases)); - ASSERT (str_len >= 1); - - big_base = __mp_bases[base].big_base; - chars_per_limb = __mp_bases[base].chars_per_limb; - - size = 0; - if (POW2_P (base)) { /* The base is a power of 2. Read the input string from least to most @@ -75,8 +79,11 @@ mpn_set_str (mp_ptr rp, const unsigned char *str, size_t str_len, int base) const unsigned char *s; int next_bitpos; - int bits_per_indigit = big_base; + mp_limb_t res_digit; + mp_size_t size; + int bits_per_indigit = __mp_bases[base].big_base; + size = 0; res_digit = 0; next_bitpos = 0; @@ -98,251 +105,246 @@ mpn_set_str (mp_ptr rp, const unsigned char *str, size_t str_len, int base) rp[size++] = res_digit; return size; } + + if (BELOW_THRESHOLD (str_len, SET_STR_PRECOMPUTE_THRESHOLD)) + return mpn_bc_set_str (rp, str, str_len, base); else { - /* General case. The base is not a power of 2. */ + mp_ptr powtab_mem, tp; + powers_t powtab[GMP_LIMB_BITS]; + int chars_per_limb; + mp_size_t size; + mp_size_t un; + TMP_DECL; - if (str_len < SET_STR_THRESHOLD) - { - size_t i; - int j; - mp_limb_t cy_limb; + TMP_MARK; - for (i = chars_per_limb; i < str_len; i += chars_per_limb) - { - res_digit = *str++; - if (base == 10) - { /* This is a common case. - Help the compiler to avoid multiplication. */ - for (j = MP_BASES_CHARS_PER_LIMB_10 - 1; j != 0; j--) - res_digit = res_digit * 10 + *str++; - } - else - { - for (j = chars_per_limb - 1; j != 0; j--) - res_digit = res_digit * base + *str++; - } - - if (size == 0) - { - if (res_digit != 0) - { - rp[0] = res_digit; - size = 1; - } - } - else - { -#if HAVE_NATIVE_mpn_mul_1c - cy_limb = mpn_mul_1c (rp, rp, size, big_base, res_digit); -#else - cy_limb = mpn_mul_1 (rp, rp, size, big_base); - cy_limb += mpn_add_1 (rp, rp, size, res_digit); -#endif - if (cy_limb != 0) - rp[size++] = cy_limb; - } - } + chars_per_limb = __mp_bases[base].chars_per_limb; - big_base = base; - res_digit = *str++; - if (base == 10) - { /* This is a common case. - Help the compiler to avoid multiplication. */ - for (j = str_len - (i - MP_BASES_CHARS_PER_LIMB_10) - 1; j > 0; j--) - { - res_digit = res_digit * 10 + *str++; - big_base *= 10; - } - } - else - { - for (j = str_len - (i - chars_per_limb) - 1; j > 0; j--) - { - res_digit = res_digit * base + *str++; - big_base *= base; - } - } + un = str_len / chars_per_limb + 1; - if (size == 0) - { - if (res_digit != 0) - { - rp[0] = res_digit; - size = 1; - } - } - else - { -#if HAVE_NATIVE_mpn_mul_1c - cy_limb = mpn_mul_1c (rp, rp, size, big_base, res_digit); + /* Allocate one large block for the powers of big_base. */ + powtab_mem = TMP_BALLOC_LIMBS (mpn_dc_set_str_powtab_alloc (un)); + + mpn_set_str_compute_powtab (powtab, powtab_mem, un, base); + + tp = TMP_BALLOC_LIMBS (mpn_dc_set_str_itch (un)); + size = mpn_dc_set_str (rp, str, str_len, powtab, tp); + + TMP_FREE; + return size; + } +} + +void +mpn_set_str_compute_powtab (powers_t *powtab, mp_ptr powtab_mem, mp_size_t un, int base) +{ + mp_ptr powtab_mem_ptr; + long i, pi; + mp_size_t n; + mp_ptr p, t; + unsigned normalization_steps; + mp_limb_t big_base, big_base_inverted; + int chars_per_limb; + size_t digits_in_base; + + powtab_mem_ptr = powtab_mem; + + chars_per_limb = __mp_bases[base].chars_per_limb; + big_base = __mp_bases[base].big_base; + big_base_inverted = __mp_bases[base].big_base_inverted; + count_leading_zeros (normalization_steps, big_base); + + p = powtab_mem_ptr; + powtab_mem_ptr += 1; + + digits_in_base = chars_per_limb; + + p[0] = big_base; + n = 1; + + count_leading_zeros (i, un - 1); + i = GMP_LIMB_BITS - 1 - i; + + powtab[i].p = p; + powtab[i].n = n; + powtab[i].digits_in_base = digits_in_base; + powtab[i].base = base; + + for (pi = i - 1; pi >= 0; pi--) + { + t = powtab_mem_ptr; + powtab_mem_ptr += 2 * n; + + ASSERT_ALWAYS (powtab_mem_ptr < powtab_mem + mpn_dc_set_str_powtab_alloc (un)); + + mpn_sqr_n (t, p, n); + n = 2 * n - 1; n += t[n] != 0; + digits_in_base *= 2; +#if 1 + if ((((un - 1) >> pi) & 2) == 0) + { + mpn_divexact_1 (t, t, n, big_base); + n -= t[n - 1] == 0; + digits_in_base -= chars_per_limb; + } #else - cy_limb = mpn_mul_1 (rp, rp, size, big_base); - cy_limb += mpn_add_1 (rp, rp, size, res_digit); -#endif - if (cy_limb != 0) - rp[size++] = cy_limb; - } - return size; + if (CLEVER_CONDITION_1 ()) + { + // perform adjustment operation of previous + cy = mpn_mul_1 (p, p, n, big_base); } - else + if (CLEVER_CONDITION_2 ()) { - /* Sub-quadratic code. */ - - mp_ptr dp; - mp_size_t dsize; - mp_ptr xp, tp; - mp_size_t step; - mp_size_t i; - size_t alloc; - mp_size_t n; - mp_ptr pow_mem; - - alloc = (str_len + chars_per_limb - 1) / chars_per_limb; - alloc = 2 * alloc; - dp = __GMP_ALLOCATE_FUNC_LIMBS (alloc); - -#if SET_STR_BLOCK_SIZE == 1 - dsize = convert_blocks (dp, str, str_len, base); -#else - { - const unsigned char *s; - mp_ptr ddp = dp; - - s = str + str_len; - while (s - str > SET_STR_BLOCK_SIZE * chars_per_limb) - { - s -= SET_STR_BLOCK_SIZE * chars_per_limb; - mpn_set_str (ddp, s, SET_STR_BLOCK_SIZE * chars_per_limb, base); - ddp += SET_STR_BLOCK_SIZE; - } - ddp += mpn_set_str (ddp, str, s - str, base); - dsize = ddp - dp; - } + // perform adjustment operation of new + cy = mpn_mul_1 (t, t, n, big_base); + } #endif + p = t; + powtab[pi].p = p; + powtab[pi].n = n; + powtab[pi].digits_in_base = digits_in_base; + powtab[pi].base = base; + } +} - /* Allocate space for powers of big_base. Could trim this in two - ways: - 1. Only really need 2^ceil(log2(dsize)) bits for the largest - power. - 2. Only the variable to get the largest power needs that much - memory. The other variable needs half as much. We need just - figure out which of xp and tp will hold the last one. - Net space savings would be in the range 1/4 to 5/8 of current - allocation, depending on how close dsize is to the next greater - power of 2. */ - pow_mem = __GMP_ALLOCATE_FUNC_LIMBS (2 * alloc); - xp = pow_mem; - tp = pow_mem + alloc; - - xp[0] = big_base; - n = 1; - step = 1; -#if SET_STR_BLOCK_SIZE != 1 - for (i = SET_STR_BLOCK_SIZE; i > 1; i >>= 1) - { - mpn_sqr_n (tp, xp, n); - n = 2 * n; - n -= tp[n - 1] == 0; +mp_size_t +mpn_dc_set_str (mp_ptr rp, const unsigned char *str, size_t str_len, + const powers_t *powtab, mp_ptr tp) +{ + size_t len_lo, len_hi; + mp_limb_t cy; + mp_size_t ln, hn, n; - step = 2 * step; - MP_PTR_SWAP (tp, xp); - } -#endif + len_lo = powtab->digits_in_base; + len_hi = str_len - len_lo; - /* Multiply every second limb block, each `step' limbs large by the - base power currently in xp[], then add this to the adjacent block. - We thereby convert from dsize blocks in base big_base, to dsize/2 - blocks in base big_base^2, then to dsize/4 blocks in base - big_base^4, etc, etc. */ + if (str_len <= len_lo) + { + if (BELOW_THRESHOLD (str_len, SET_STR_DC_THRESHOLD)) + return mpn_bc_set_str (rp, str, str_len, powtab->base); + else + return mpn_dc_set_str (rp, str, str_len, powtab + 1, tp); + } - if (step < dsize) - { - for (;;) - { - for (i = 0; i < dsize - step; i += 2 * step) - { - mp_ptr bp = dp + i; - mp_size_t m = dsize - i - step; - mp_size_t hi; - if (n >= m) - { - mpn_mul (tp, xp, n, bp + step, m); - mpn_add (bp, tp, n + m, bp, n); - hi = i + n + m; - dsize = hi; - dsize -= dp[dsize - 1] == 0; - } - else - { - mpn_mul_n (tp, xp, bp + step, n); - mpn_add (bp, tp, n + n, bp, n); - } - } - - step = 2 * step; - if (! (step < dsize)) - break; - - mpn_sqr_n (tp, xp, n); - n = 2 * n; - n -= tp[n - 1] == 0; - MP_PTR_SWAP (tp, xp); - } - } + ASSERT_ALWAYS (len_lo >= len_hi); - MPN_NORMALIZE (dp, dsize); - MPN_COPY (rp, dp, dsize); - __GMP_FREE_FUNC_LIMBS (pow_mem, 2 * alloc); - __GMP_FREE_FUNC_LIMBS (dp, alloc); - return dsize; - } - } + if (BELOW_THRESHOLD (len_hi, SET_STR_DC_THRESHOLD)) + hn = mpn_bc_set_str (tp, str, len_hi, powtab->base); + else + hn = mpn_dc_set_str (tp, str, len_hi, powtab + 1, rp); + + if (hn == 0) + MPN_ZERO (rp, powtab->n); + else + mpn_mul (rp, powtab->p, powtab->n, tp, hn); + + str = str + str_len - len_lo; + if (BELOW_THRESHOLD (len_lo, SET_STR_DC_THRESHOLD)) + ln = mpn_bc_set_str (tp, str, len_lo, powtab->base); + else + ln = mpn_dc_set_str (tp, str, len_lo, powtab + 1, tp + powtab->n + 1); + + cy = mpn_add_n (rp, rp, tp, ln); + mpn_incr_u (rp + ln, cy); + n = hn + powtab->n; + return n - (rp[n - 1] == 0); } -static mp_size_t -convert_blocks (mp_ptr dp, const unsigned char *str, size_t str_len, int base) +mp_size_t +mpn_bc_set_str (mp_ptr rp, const unsigned char *str, size_t str_len, int base) { + mp_size_t size; + size_t i; + long j; + mp_limb_t cy_limb; + + mp_limb_t big_base; int chars_per_limb; - mp_size_t i; - int j; - int ds; - mp_size_t dsize; mp_limb_t res_digit; - chars_per_limb = __mp_bases[base].chars_per_limb; + ASSERT (base >= 2); + ASSERT (base < numberof (__mp_bases)); + ASSERT (str_len >= 1); - dsize = str_len / chars_per_limb; - ds = str_len % chars_per_limb; + big_base = __mp_bases[base].big_base; + chars_per_limb = __mp_bases[base].chars_per_limb; - if (ds != 0) + size = 0; + for (i = chars_per_limb; i < str_len; i += chars_per_limb) { res_digit = *str++; - for (j = ds - 1; j != 0; j--) - res_digit = res_digit * base + *str++; - dp[dsize] = res_digit; + if (base == 10) + { /* This is a common case. + Help the compiler to avoid multiplication. */ + for (j = MP_BASES_CHARS_PER_LIMB_10 - 1; j != 0; j--) + res_digit = res_digit * 10 + *str++; + } + else + { + for (j = chars_per_limb - 1; j != 0; j--) + res_digit = res_digit * base + *str++; + } + + if (size == 0) + { + if (res_digit != 0) + { + rp[0] = res_digit; + size = 1; + } + } + else + { +#if HAVE_NATIVE_mpn_mul_1c + cy_limb = mpn_mul_1c (rp, rp, size, big_base, res_digit); +#else + cy_limb = mpn_mul_1 (rp, rp, size, big_base); + cy_limb += mpn_add_1 (rp, rp, size, res_digit); +#endif + if (cy_limb != 0) + rp[size++] = cy_limb; + } } + big_base = base; + res_digit = *str++; if (base == 10) - { - for (i = dsize - 1; i >= 0; i--) + { /* This is a common case. + Help the compiler to avoid multiplication. */ + for (j = str_len - (i - MP_BASES_CHARS_PER_LIMB_10) - 1; j > 0; j--) { - res_digit = *str++; - for (j = MP_BASES_CHARS_PER_LIMB_10 - 1; j != 0; j--) - res_digit = res_digit * 10 + *str++; - dp[i] = res_digit; + res_digit = res_digit * 10 + *str++; + big_base *= 10; } } else { - for (i = dsize - 1; i >= 0; i--) + for (j = str_len - (i - chars_per_limb) - 1; j > 0; j--) { - res_digit = *str++; - for (j = chars_per_limb - 1; j != 0; j--) - res_digit = res_digit * base + *str++; - dp[i] = res_digit; + res_digit = res_digit * base + *str++; + big_base *= base; } } - return dsize + (ds != 0); + if (size == 0) + { + if (res_digit != 0) + { + rp[0] = res_digit; + size = 1; + } + } + else + { +#if HAVE_NATIVE_mpn_mul_1c + cy_limb = mpn_mul_1c (rp, rp, size, big_base, res_digit); +#else + cy_limb = mpn_mul_1 (rp, rp, size, big_base); + cy_limb += mpn_add_1 (rp, rp, size, res_digit); +#endif + if (cy_limb != 0) + rp[size++] = cy_limb; + } + return size; } |