diff --git a/ext/bcmath/libbcmath/src/bcmath.h b/ext/bcmath/libbcmath/src/bcmath.h index 1f05ad51f7f26..3a0b80dfc6893 100644 --- a/ext/bcmath/libbcmath/src/bcmath.h +++ b/ext/bcmath/libbcmath/src/bcmath.h @@ -147,8 +147,6 @@ bc_num bc_multiply(bc_num n1, bc_num n2, size_t scale); *(result) = mul_ex; \ } while (0) -bc_num bc_square(bc_num n1, size_t scale); - bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, size_t scale); bool bc_modulo(bc_num num1, bc_num num2, bc_num *resul, size_t scale); diff --git a/ext/bcmath/libbcmath/src/private.h b/ext/bcmath/libbcmath/src/private.h index 91facfb2f8b40..327a72c3677af 100644 --- a/ext/bcmath/libbcmath/src/private.h +++ b/ext/bcmath/libbcmath/src/private.h @@ -84,6 +84,8 @@ static const BC_VECTOR BC_POW_10_LUT[9] = { bcmath_compare_result _bc_do_compare (bc_num n1, bc_num n2, size_t scale, bool use_sign); bc_num _bc_do_add (bc_num n1, bc_num n2); bc_num _bc_do_sub (bc_num n1, bc_num n2); +void bc_multiply_vector( + BC_VECTOR *n1_vector, size_t n1_arr_size, BC_VECTOR *n2_vector, size_t n2_arr_size, BC_VECTOR *prod_vector, size_t prod_arr_size); void _bc_rm_leading_zeros (bc_num num); #endif diff --git a/ext/bcmath/libbcmath/src/raise.c b/ext/bcmath/libbcmath/src/raise.c index 1e283864694b6..f9d9dacc35d1f 100644 --- a/ext/bcmath/libbcmath/src/raise.c +++ b/ext/bcmath/libbcmath/src/raise.c @@ -30,24 +30,152 @@ *************************************************************************/ #include "bcmath.h" +#include "convert.h" +#include "private.h" #include #include #include -void bc_square_ex(bc_num n1, bc_num *result, size_t scale_min) { - bc_num square_ex = bc_square(n1, scale_min); - bc_free_num(result); - *(result) = square_ex; +static inline size_t bc_multiply_vector_ex( + BC_VECTOR **n1_vector, size_t n1_arr_size, BC_VECTOR *n2_vector, size_t n2_arr_size, BC_VECTOR **result_vector) +{ + size_t result_arr_size = n1_arr_size + n2_arr_size; + bc_multiply_vector(*n1_vector, n1_arr_size, n2_vector, n2_arr_size, *result_vector, result_arr_size); + + /* Eliminate extra zeros because they increase the number of calculations. */ + while ((*result_vector)[result_arr_size - 1] == 0) { + result_arr_size--; + } + + /* Swap n1_vector and result_vector. */ + BC_VECTOR *tmp = *n1_vector; + *n1_vector = *result_vector; + *result_vector = tmp; + + return result_arr_size; +} + +static inline size_t bc_square_vector_ex(BC_VECTOR **base_vector, size_t base_arr_size, BC_VECTOR **result_vector) +{ + return bc_multiply_vector_ex(base_vector, base_arr_size, *base_vector, base_arr_size, result_vector); +} + +/* Use "exponentiation by squaring". This is the fast path when the results are small. */ +static inline bc_num bc_fast_raise( + const char *base_end, long exponent, size_t base_len, size_t power_len, size_t power_scale, size_t power_full_len) +{ + BC_VECTOR base_vector = 0; + + /* Convert to BC_VECTOR[] */ + bc_convert_to_vector(&base_vector, base_end, base_len); + + while ((exponent & 1) == 0) { + base_vector *= base_vector; + exponent >>= 1; + } + + /* copy base to power */ + BC_VECTOR power_vector = base_vector; + exponent >>= 1; + + while (exponent > 0) { + base_vector *= base_vector; + if ((exponent & 1) == 1) { + power_vector *= base_vector; + } + exponent >>= 1; + } + + bc_num power = bc_new_num_nonzeroed(power_len, power_scale); + char *pptr = power->n_value; + char *pend = pptr + power_full_len - 1; + + while (pend >= pptr) { + *pend-- = power_vector % BASE; + power_vector /= BASE; + } + return power; +} + +/* Use "exponentiation by squaring". This is the standard path. */ +static bc_num bc_standard_raise( + const char *base_ptr, const char *base_end, long exponent, size_t base_len, size_t power_scale) +{ + /* Remove the leading zeros as they will be filled in later. */ + while (*base_ptr++ == 0) { + base_len--; + } + + size_t base_arr_size = BC_ARR_SIZE_FROM_LEN(base_len); + size_t max_power_arr_size = base_arr_size * exponent; + + /* The allocated memory area is reused on a rotational basis, so the same size is required. */ + BC_VECTOR *buf = safe_emalloc(max_power_arr_size * 3, sizeof(BC_VECTOR), 0); + BC_VECTOR *base_vector = buf; + BC_VECTOR *power_vector = base_vector + max_power_arr_size; + BC_VECTOR *tmp_result_vector = power_vector + max_power_arr_size; + + /* Convert to BC_VECTOR[] */ + bc_convert_to_vector(base_vector, base_end, base_len); + + while ((exponent & 1) == 0) { + base_arr_size = bc_square_vector_ex(&base_vector, base_arr_size, &tmp_result_vector); + exponent >>= 1; + } + + /* copy base to power */ + size_t power_arr_size = base_arr_size; + for (size_t i = 0; i < base_arr_size; i++) { + power_vector[i] = base_vector[i]; + } + exponent >>= 1; + + while (exponent > 0) { + base_arr_size = bc_square_vector_ex(&base_vector, base_arr_size, &tmp_result_vector); + if ((exponent & 1) == 1) { + power_arr_size = bc_multiply_vector_ex(&power_vector, power_arr_size, base_vector, base_arr_size, &tmp_result_vector); + } + exponent >>= 1; + } + + /* Convert to bc_num */ + size_t power_leading_zeros = 0; + size_t power_len; + size_t power_full_len = power_arr_size * BC_VECTOR_SIZE; + if (power_full_len > power_scale) { + power_len = power_full_len - power_scale; + } else { + power_len = 1; + power_leading_zeros = power_scale - power_full_len + 1; + power_full_len = power_scale + 1; + } + bc_num power = bc_new_num_nonzeroed(power_len, power_scale); + + char *pptr = power->n_value; + char *pend = pptr + power_full_len - 1; + + /* Pad with leading zeros if necessary. */ + while (power_leading_zeros > sizeof(uint32_t)) { + bc_write_bcd_representation(0, pptr); + pptr += sizeof(uint32_t); + power_leading_zeros -= sizeof(uint32_t); + } + for (size_t i = 0; i < power_leading_zeros; i++) { + *pptr++ = 0; + } + + bc_convert_vector_to_char(power_vector, pptr, pend, power_arr_size); + + efree(buf); + + return power; } /* Raise "base" to the "exponent" power. The result is placed in RESULT. Maximum exponent is LONG_MAX. If a "exponent" is not an integer, only the integer part is used. */ bool bc_raise(bc_num base, long exponent, bc_num *result, size_t scale) { - bc_num temp, power; size_t rscale; - size_t pwrscale; - size_t calcscale; bool is_neg; /* Special case if exponent is a zero. */ @@ -67,43 +195,54 @@ bool bc_raise(bc_num base, long exponent, bc_num *result, size_t scale) { rscale = MIN (base->n_scale * exponent, MAX(scale, base->n_scale)); } - /* Set initial value of temp. */ - power = bc_copy_num(base); - pwrscale = base->n_scale; - while ((exponent & 1) == 0) { - pwrscale = 2 * pwrscale; - bc_square_ex(power, &power, pwrscale); - exponent = exponent >> 1; + if (bc_is_zero(base)) { + bc_free_num(result); + *result = bc_copy_num(BCG(_zero_)); + /* If the exponent is negative, it divides by 0, so it is false. */ + return !is_neg; } - temp = bc_copy_num(power); - calcscale = pwrscale; - exponent = exponent >> 1; - /* Do the calculation. */ - while (exponent > 0) { - pwrscale = 2 * pwrscale; - bc_square_ex(power, &power, pwrscale); - if ((exponent & 1) == 1) { - calcscale = pwrscale + calcscale; - bc_multiply_ex(temp, power, &temp, calcscale); - } - exponent = exponent >> 1; + size_t base_len = base->n_len + base->n_scale; + size_t power_len = base->n_len * exponent; + size_t power_scale = base->n_scale * exponent; + size_t power_full_len = power_len + power_scale; + + sign power_sign; + if (base->n_sign == MINUS && (exponent & 1) == 1) { + power_sign = MINUS; + } else { + power_sign = PLUS; + } + + const char *base_end = base->n_value + base_len - 1; + + bc_num power; + if (base_len <= BC_VECTOR_SIZE && power_full_len <= BC_VECTOR_SIZE * 2) { + power = bc_fast_raise(base_end, exponent, base_len, power_len, power_scale, power_full_len); + } else { + power = bc_standard_raise(base->n_value, base_end, exponent, base_len, power_scale); + } + + _bc_rm_leading_zeros(power); + if (bc_is_zero(power)) { + power->n_sign = PLUS; + power->n_scale = 0; + } else { + power->n_sign = power_sign; } /* Assign the value. */ if (is_neg) { - if (bc_divide(BCG(_one_), temp, result, rscale) == false) { - bc_free_num (&temp); + if (bc_divide(BCG(_one_), power, result, rscale) == false) { bc_free_num (&power); return false; } - bc_free_num (&temp); + bc_free_num (&power); } else { bc_free_num (result); - *result = temp; + *result = power; (*result)->n_scale = MIN(scale, (*result)->n_scale); } - bc_free_num (&power); return true; } diff --git a/ext/bcmath/libbcmath/src/recmul.c b/ext/bcmath/libbcmath/src/recmul.c index 26ce1641db410..64934cd4c3217 100644 --- a/ext/bcmath/libbcmath/src/recmul.c +++ b/ext/bcmath/libbcmath/src/recmul.c @@ -72,43 +72,36 @@ static inline void bc_fast_mul(bc_num n1, size_t n1len, bc_num n2, size_t n2len, } } -/* - * Equivalent of bc_fast_mul for small numbers to perform computations - * without using array. - */ -static inline void bc_fast_square(bc_num n1, size_t n1len, bc_num *prod) +static inline void bc_standard_vector_mul( + BC_VECTOR *n1_vector, size_t n1_arr_size, BC_VECTOR *n2_vector, size_t n2_arr_size, BC_VECTOR *prod_vector, size_t prod_arr_size) { - const char *n1end = n1->n_value + n1len - 1; - - BC_VECTOR n1_vector = bc_partial_convert_to_vector(n1end, n1len); - BC_VECTOR prod_vector = n1_vector * n1_vector; - - size_t prodlen = n1len + n1len; - *prod = bc_new_num_nonzeroed(prodlen, 0); - char *pptr = (*prod)->n_value; - char *pend = pptr + prodlen - 1; + for (size_t i = 0; i < prod_arr_size; i++) { + prod_vector[i] = 0; + } - while (pend >= pptr) { - *pend-- = prod_vector % BASE; - prod_vector /= BASE; + /* Multiplication and addition */ + size_t count = 0; + for (size_t i = 0; i < n1_arr_size; i++) { + /* + * This calculation adds the result multiple times to the array entries. + * When multiplying large numbers of digits, there is a possibility of + * overflow, so digit adjustment is performed beforehand. + */ + if (UNEXPECTED(count >= BC_VECTOR_NO_OVERFLOW_ADD_COUNT)) { + bc_mul_carry_calc(prod_vector, prod_arr_size); + count = 0; + } + count++; + for (size_t j = 0; j < n2_arr_size; j++) { + prod_vector[i + j] += n1_vector[i] * n2_vector[j]; + } } -} -/* Common part of functions bc_standard_mul and bc_standard_square - * that takes a vector and converts it to a bc_num */ -static inline void bc_mul_finish_from_vector(BC_VECTOR *prod_vector, size_t prod_arr_size, size_t prodlen, bc_num *prod) { /* * Move a value exceeding 4/8 digits by carrying to the next digit. * However, the last digit does nothing. */ bc_mul_carry_calc(prod_vector, prod_arr_size); - - /* Convert to bc_num */ - *prod = bc_new_num_nonzeroed(prodlen, 0); - char *pptr = (*prod)->n_value; - char *pend = pptr + prodlen - 1; - - bc_convert_vector_to_char(prod_vector, pptr, pend, prod_arr_size); } /* @@ -121,7 +114,6 @@ static inline void bc_mul_finish_from_vector(BC_VECTOR *prod_vector, size_t prod */ static void bc_standard_mul(bc_num n1, size_t n1len, bc_num n2, size_t n2len, bc_num *prod) { - size_t i; const char *n1end = n1->n_value + n1len - 1; const char *n2end = n2->n_value + n2len - 1; size_t prodlen = n1len + n2len; @@ -147,88 +139,24 @@ static void bc_standard_mul(bc_num n1, size_t n1len, bc_num n2, size_t n2len, bc BC_VECTOR *n2_vector = n1_vector + n1_arr_size; BC_VECTOR *prod_vector = n2_vector + n2_arr_size; - for (i = 0; i < prod_arr_size; i++) { - prod_vector[i] = 0; - } - /* Convert to BC_VECTOR[] */ bc_convert_to_vector(n1_vector, n1end, n1len); bc_convert_to_vector(n2_vector, n2end, n2len); - /* Multiplication and addition */ - size_t count = 0; - for (i = 0; i < n1_arr_size; i++) { - /* - * This calculation adds the result multiple times to the array entries. - * When multiplying large numbers of digits, there is a possibility of - * overflow, so digit adjustment is performed beforehand. - */ - if (UNEXPECTED(count >= BC_VECTOR_NO_OVERFLOW_ADD_COUNT)) { - bc_mul_carry_calc(prod_vector, prod_arr_size); - count = 0; - } - count++; - for (size_t j = 0; j < n2_arr_size; j++) { - prod_vector[i + j] += n1_vector[i] * n2_vector[j]; - } - } + /* Do multiply */ + bc_standard_vector_mul(n1_vector, n1_arr_size, n2_vector, n2_arr_size, prod_vector, prod_arr_size); - bc_mul_finish_from_vector(prod_vector, prod_arr_size, prodlen, prod); + /* Convert to bc_num */ + *prod = bc_new_num_nonzeroed(prodlen, 0); + char *pptr = (*prod)->n_value; + char *pend = pptr + prodlen - 1; + bc_convert_vector_to_char(prod_vector, pptr, pend, prod_arr_size); if (allocation_arr_size > BC_STACK_VECTOR_SIZE) { efree(n1_vector); } } -/** This is bc_standard_mul implementation for square */ -static void bc_standard_square(bc_num n1, size_t n1len, bc_num *prod) -{ - size_t i; - const char *n1end = n1->n_value + n1len - 1; - size_t prodlen = n1len + n1len; - - size_t n1_arr_size = BC_ARR_SIZE_FROM_LEN(n1len); - size_t prod_arr_size = BC_ARR_SIZE_FROM_LEN(prodlen); - - BC_VECTOR *buf = safe_emalloc(n1_arr_size + n1_arr_size + prod_arr_size, sizeof(BC_VECTOR), 0); - - BC_VECTOR *n1_vector = buf; - BC_VECTOR *prod_vector = n1_vector + n1_arr_size + n1_arr_size; - - for (i = 0; i < prod_arr_size; i++) { - prod_vector[i] = 0; - } - - /* Convert to BC_VECTOR[] */ - bc_convert_to_vector(n1_vector, n1end, n1len); - - /* Multiplication and addition */ - size_t count = 0; - for (i = 0; i < n1_arr_size; i++) { - /* - * This calculation adds the result multiple times to the array entries. - * When multiplying large numbers of digits, there is a possibility of - * overflow, so digit adjustment is performed beforehand. - */ - if (UNEXPECTED(count >= BC_VECTOR_NO_OVERFLOW_ADD_COUNT)) { - bc_mul_carry_calc(prod_vector, prod_arr_size); - count = 0; - } - count++; - for (size_t j = 0; j < n1_arr_size; j++) { - prod_vector[i + j] += n1_vector[i] * n1_vector[j]; - } - } - - bc_mul_finish_from_vector(prod_vector, prod_arr_size, prodlen, prod); - - efree(buf); -} - -/* The multiply routine. N2 times N1 is put int PROD with the scale of - the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). - */ - bc_num bc_multiply(bc_num n1, bc_num n2, size_t scale) { bc_num prod; @@ -258,24 +186,16 @@ bc_num bc_multiply(bc_num n1, bc_num n2, size_t scale) return prod; } -bc_num bc_square(bc_num n1, size_t scale) +void bc_multiply_vector( + BC_VECTOR *n1_vector, size_t n1_arr_size, BC_VECTOR *n2_vector, size_t n2_arr_size, BC_VECTOR *prod_vector, size_t prod_arr_size) { - bc_num prod; - - size_t len1 = n1->n_len + n1->n_scale; - size_t full_scale = n1->n_scale + n1->n_scale; - size_t prod_scale = MIN(full_scale, MAX(scale, n1->n_scale)); - - if (len1 <= BC_VECTOR_SIZE) { - bc_fast_square(n1, len1, &prod); + if (n1_arr_size == 1 && n2_arr_size == 1) { + prod_vector[0] = *n1_vector * *n2_vector; + if (prod_arr_size == 2) { + prod_vector[1] = prod_vector[0] / BC_VECTOR_BOUNDARY_NUM; + prod_vector[0] %= BC_VECTOR_BOUNDARY_NUM; + } } else { - bc_standard_square(n1, len1, &prod); + bc_standard_vector_mul(n1_vector, n1_arr_size, n2_vector, n2_arr_size, prod_vector, prod_arr_size); } - - prod->n_sign = PLUS; - prod->n_len -= full_scale; - prod->n_scale = prod_scale; - _bc_rm_leading_zeros(prod); - - return prod; }