> O(N^2) where N = max(abs(exponent of x), precision).
And this patch optimizes the case of huge precision. So that the complexity
decreases from
O(abs(exponent of x)^2 + precision^2)
down to
O(abs(exponent of x)^2 + precision)
2007-05-19 Bruno Haible <[EMAIL PROTECTED]>
* lib/vasnprintf.c (convert_to_decimal): Add an extra_zeroes argument.
(scale10_round_decimal_long_double): Inline scale10_round_long_double.
Instead of multiplying with 10^k, set extra_zeroes to k.
(scale10_round_long_double): Remove function.
*** lib/vasnprintf.c 19 May 2007 00:38:42 -0000 1.47
--- lib/vasnprintf.c 19 May 2007 09:41:59 -0000
***************
*** 656,677 ****
return roomptr;
}
! /* Convert a bignum a >= 0 to decimal representation.
Destroys the contents of a.
Return the allocated memory - containing the decimal digits in low-to-high
order, terminated with a NUL character - in case of success, NULL in case
of memory allocation failure. */
static char *
! convert_to_decimal (mpn_t a)
{
mp_limb_t *a_ptr = a.limbs;
size_t a_len = a.nlimbs;
/* 0.03345 is slightly larger than log(2)/(9*log(10)). */
size_t c_len = 9 * ((size_t)(a_len * (GMP_LIMB_BITS * 0.03345f)) + 1);
! char *c_ptr = (char *) malloc (c_len);
if (c_ptr != NULL)
{
char *d_ptr = c_ptr;
while (a_len > 0)
{
/* Divide a by 10^9, in-place. */
--- 656,680 ----
return roomptr;
}
! /* Convert a bignum a >= 0, multiplied with 10^extra_zeroes, to decimal
! representation.
Destroys the contents of a.
Return the allocated memory - containing the decimal digits in low-to-high
order, terminated with a NUL character - in case of success, NULL in case
of memory allocation failure. */
static char *
! convert_to_decimal (mpn_t a, size_t extra_zeroes)
{
mp_limb_t *a_ptr = a.limbs;
size_t a_len = a.nlimbs;
/* 0.03345 is slightly larger than log(2)/(9*log(10)). */
size_t c_len = 9 * ((size_t)(a_len * (GMP_LIMB_BITS * 0.03345f)) + 1);
! char *c_ptr = (char *) malloc (xsum (c_len, extra_zeroes));
if (c_ptr != NULL)
{
char *d_ptr = c_ptr;
+ for (; extra_zeroes > 0; extra_zeroes--)
+ *d_ptr++ = '0';
while (a_len > 0)
{
/* Divide a by 10^9, in-place. */
***************
*** 789,804 ****
}
/* Assuming x is finite and >= 0, and n is an integer:
! Compute y = round (x * 10^n) as a bignum >= 0.
! Return the allocated memory in case of success, NULL in case of memory
! allocation failure. */
! static void *
! scale10_round_long_double (long double x, int n, mpn_t *yp)
{
int e;
mpn_t m;
void *memory = decode_long_double (x, &e, &m);
int s;
unsigned int abs_n;
unsigned int abs_s;
mp_limb_t *pow5_ptr;
--- 792,809 ----
}
/* Assuming x is finite and >= 0, and n is an integer:
! Returns the decimal representation of round (x * 10^n).
! Return the allocated memory - containing the decimal digits in low-to-high
! order, terminated with a NUL character - in case of success, NULL in case
! of memory allocation failure. */
! static char *
! scale10_round_decimal_long_double (long double x, int n)
{
int e;
mpn_t m;
void *memory = decode_long_double (x, &e, &m);
int s;
+ size_t extra_zeroes;
unsigned int abs_n;
unsigned int abs_s;
mp_limb_t *pow5_ptr;
***************
*** 806,811 ****
--- 811,819 ----
unsigned int s_limbs;
unsigned int s_bits;
mpn_t pow5;
+ mpn_t z;
+ void *z_memory;
+ char *digits;
if (memory == NULL)
return NULL;
***************
*** 813,818 ****
--- 821,837 ----
y = round (2^e * 10^n * m) = round (2^(e+n) * 5^n * m)
= round (2^s * 5^n * m). */
s = e + n;
+ extra_zeroes = 0;
+ /* Factor out a common power of 10 if possible. */
+ if (s > 0 && n > 0)
+ {
+ extra_zeroes = (s < n ? s : n);
+ s -= extra_zeroes;
+ n -= extra_zeroes;
+ }
+ /* Here y = round (2^s * 5^n * m) * 10^extra_zeroes.
+ Before converting to decimal, we need to compute
+ z = round (2^s * 5^n * m). */
/* Compute 5^|n|, possibly shifted by |s| bits if n and s have the same
sign. 2.322 is slightly larger than log(5)/log(2). */
abs_n = (n >= 0 ? n : -n);
***************
*** 895,912 ****
if (n >= 0)
{
/* Multiply m with pow5. No division needed. */
! void *result_memory = multiply (m, pow5, yp);
! free (pow5_ptr);
! free (memory);
! return result_memory;
}
else
{
/* Divide m by pow5 and round. */
! void *result_memory = divide (m, pow5, yp);
! free (pow5_ptr);
! free (memory);
! return result_memory;
}
}
else
--- 914,925 ----
if (n >= 0)
{
/* Multiply m with pow5. No division needed. */
! z_memory = multiply (m, pow5, &z);
}
else
{
/* Divide m by pow5 and round. */
! z_memory = divide (m, pow5, &z);
}
}
else
***************
*** 920,926 ****
mpn_t numerator;
mpn_t denominator;
void *tmp_memory;
- void *result_memory;
tmp_memory = multiply (m, pow5, &numerator);
if (tmp_memory == NULL)
{
--- 933,938 ----
***************
*** 938,948 ****
denominator.limbs = ptr;
denominator.nlimbs = s_limbs + 1;
}
! result_memory = divide (numerator, denominator, yp);
free (tmp_memory);
- free (pow5_ptr);
- free (memory);
- return result_memory;
}
else
{
--- 950,957 ----
denominator.limbs = ptr;
denominator.nlimbs = s_limbs + 1;
}
! z_memory = divide (numerator, denominator, &z);
free (tmp_memory);
}
else
{
***************
*** 950,956 ****
Multiply m with 2^s, then divide by pow5. */
mpn_t numerator;
mp_limb_t *num_ptr;
- void *result_memory;
num_ptr = (mp_limb_t *) malloc ((m.nlimbs + s_limbs + 1)
* sizeof (mp_limb_t));
if (num_ptr == NULL)
--- 959,964 ----
***************
*** 990,1021 ****
numerator.limbs = num_ptr;
numerator.nlimbs = destptr - num_ptr;
}
! result_memory = divide (numerator, pow5, yp);
free (num_ptr);
- free (pow5_ptr);
- free (memory);
- return result_memory;
}
}
! }
! /* Assuming x is finite and >= 0, and n is an integer:
! Returns the decimal representation of round (x * 10^n).
! Return the allocated memory - containing the decimal digits in low-to-high
! order, terminated with a NUL character - in case of success, NULL in case
! of memory allocation failure. */
! static char *
! scale10_round_decimal_long_double (long double x, int n)
! {
! mpn_t y;
! void *memory;
! char *digits;
! memory = scale10_round_long_double (x, n, &y);
! if (memory == NULL)
return NULL;
! digits = convert_to_decimal (y);
! free (memory);
return digits;
}
--- 998,1016 ----
numerator.limbs = num_ptr;
numerator.nlimbs = destptr - num_ptr;
}
! z_memory = divide (numerator, pow5, &z);
free (num_ptr);
}
}
! free (pow5_ptr);
! free (memory);
! /* Here y = round (x * 10^n) = z * 10^extra_zeroes. */
! if (z_memory == NULL)
return NULL;
! digits = convert_to_decimal (z, extra_zeroes);
! free (z_memory);
return digits;
}
