php/archived-php-src
1
0
Fork 0
mirror of https://github.com/php/php-src.git synced 2026-03-24 00:02:20 +01:00
Code Issues Packages Projects Releases Wiki Activity

ext/bcmath: Improving bcpow() performance (#18099)

Browse Source
This commit is contained in:
Saki Takamachi
2025-04-30 08:05:28 +09:00
committed by GitHub
parent 26dbcb7640
commit c5f3281446
5 changed files with 259 additions and 161 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
ext/bcmath/bcmath.c
View File

@@ -179,6 +179,26 @@ static zend_result php_str2num(bc_num *num, const zend_string *str)
}
/* }}} */
static void bc_pow_err(bc_raise_status status, uint32_t arg_num)
{
/* If arg_num is 0, it means it is an op */
switch (status) {
case BC_RAISE_STATUS_DIVIDE_BY_ZERO:
zend_throw_exception_ex(zend_ce_division_by_zero_error, 0, "Negative power of zero");
break;
case BC_RAISE_STATUS_LEN_IS_OVERFLOW:
case BC_RAISE_STATUS_SCALE_IS_OVERFLOW:
case BC_RAISE_STATUS_FULLLEN_IS_OVERFLOW:
if (arg_num == 0) {
zend_value_error("exponent is too large, the number of digits overflowed");
} else {
zend_argument_value_error(arg_num, "exponent is too large, the number of digits overflowed");
}
break;
EMPTY_SWITCH_DEFAULT_CASE();
}
}
/* {{{ Returns the sum of two arbitrary precision numbers */
PHP_FUNCTION(bcadd)
{
@@ -615,11 +635,11 @@ PHP_FUNCTION(bcpow)
goto cleanup;
}
if (!bc_raise(first, exponent, &result, scale)) {
zend_throw_exception_ex(zend_ce_division_by_zero_error, 0, "Negative power of zero");
bc_raise_status ret_status = bc_raise(first, exponent, &result, scale);
if (UNEXPECTED(ret_status != BC_RAISE_STATUS_OK)) {
bc_pow_err(ret_status, 2);
goto cleanup;
}
RETVAL_NEW_STR(bc_num2str_ex(result, scale));
cleanup: {
@@ -1144,8 +1164,9 @@ static zend_result bcmath_number_pow_internal(
}
return FAILURE;
}
if (!bc_raise(n1, exponent, ret, *scale)) {
zend_throw_exception_ex(zend_ce_division_by_zero_error, 0, "Negative power of zero");
bc_raise_status ret_status = bc_raise(n1, exponent, ret, *scale);
if (UNEXPECTED(ret_status != BC_RAISE_STATUS_OK)) {
bc_pow_err(ret_status, is_op ? 0 : 1);
return FAILURE;
}
bc_rm_trailing_zeros(*ret);

12
ext/bcmath/libbcmath/src/bcmath.h
View File

@@ -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);
@@ -159,6 +157,14 @@ bc_num bc_floor_or_ceil(bc_num num, bool is_floor);
size_t bc_round(bc_num num, zend_long places, zend_long mode, bc_num *result);
typedef enum {
BC_RAISE_STATUS_OK,
BC_RAISE_STATUS_LEN_IS_OVERFLOW,
BC_RAISE_STATUS_SCALE_IS_OVERFLOW,
BC_RAISE_STATUS_FULLLEN_IS_OVERFLOW,
BC_RAISE_STATUS_DIVIDE_BY_ZERO,
} bc_raise_status;
typedef enum {
OK,
BASE_HAS_FRACTIONAL,
@@ -170,7 +176,7 @@ typedef enum {
raise_mod_status bc_raisemod(bc_num base, bc_num exponent, bc_num mod, bc_num *result, size_t scale);
bool bc_raise(bc_num base, long exponent, bc_num *result, size_t scale);
bc_raise_status bc_raise(bc_num base, long exponent, bc_num *result, size_t scale);
void bc_raise_bc_exponent(bc_num base, bc_num exponent, bc_num *resul, size_t scale);

3
ext/bcmath/libbcmath/src/private.h
View File

@@ -84,6 +84,9 @@ 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(
const BC_VECTOR *n1_vector, size_t n1_arr_size, const 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

216
ext/bcmath/libbcmath/src/raise.c
View File

@@ -30,31 +30,155 @@
*************************************************************************/
#include "bcmath.h"
#include "convert.h"
#include "private.h"
#include <assert.h>
#include <stdbool.h>
#include <stddef.h>
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_ptr++;
base_len--;
}
size_t base_arr_size = BC_ARR_SIZE_FROM_LEN(base_len);
/* Since it is guaranteed that base_len * exponent does not overflow, there is no possibility of overflow here. */
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, sizeof(BC_VECTOR) * 3, 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. */
memset(pptr, 0, power_leading_zeros);
pptr += power_leading_zeros;
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;
bc_raise_status bc_raise(bc_num base, long exponent, bc_num *result, size_t scale) {
size_t rscale;
size_t pwrscale;
size_t calcscale;
bool is_neg;
/* Special case if exponent is a zero. */
if (exponent == 0) {
bc_free_num (result);
*result = bc_copy_num(BCG(_one_));
return true;
return BC_RAISE_STATUS_OK;
}
/* Other initializations. */
@@ -67,44 +191,66 @@ 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)) {
/* If the exponent is negative, it divides by 0 */
return is_neg ? BC_RAISE_STATUS_DIVIDE_BY_ZERO : BC_RAISE_STATUS_OK;
}
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;
/* check overflow */
if (UNEXPECTED(base->n_len > SIZE_MAX / exponent)) {
return BC_RAISE_STATUS_LEN_IS_OVERFLOW;
}
if (UNEXPECTED(base->n_scale > SIZE_MAX / exponent)) {
return BC_RAISE_STATUS_SCALE_IS_OVERFLOW;
}
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;
/* check overflow */
if (UNEXPECTED(power_len > SIZE_MAX - power_scale)) {
return BC_RAISE_STATUS_FULLLEN_IS_OVERFLOW;
}
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;
return BC_RAISE_STATUS_DIVIDE_BY_ZERO;
}
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;
return BC_RAISE_STATUS_OK;
}
/* This is used internally by BCMath */

158
ext/bcmath/libbcmath/src/recmul.c
View File

@@ -72,43 +72,37 @@ 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(
const BC_VECTOR *n1_vector, size_t n1_arr_size, const 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;
while (pend >= pptr) {
*pend-- = prod_vector % BASE;
prod_vector /= BASE;
for (size_t i = 0; i < prod_arr_size; i++) {
prod_vector[i] = 0;
}
/* 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 +115,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 +140,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 +187,17 @@ 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(
const BC_VECTOR *n1_vector, size_t n1_arr_size, const 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;
}