Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions ext/bcmath/libbcmath/src/convert.h
Original file line number Diff line number Diff line change
Expand Up @@ -85,16 +85,16 @@ static inline void bc_convert_to_vector_with_zero_pad(BC_VECTOR *n_vector, const
bc_convert_to_vector(n_vector, nend, nlen);
}

static inline void bc_convert_vector_to_char(const BC_VECTOR *vector, char *nptr, char *nend, size_t arr_size)
static inline void bc_convert_vector_to_char(const bc_vector_arr *vector, char *nptr, char *nend)
{
size_t i = 0;
while (i < arr_size - 1) {
while (i < vector->size - 1) {
#if BC_VECTOR_SIZE == 4
bc_write_bcd_representation(vector[i], nend - 3);
nend -= 4;
#else
bc_write_bcd_representation(vector[i] / 10000, nend - 7);
bc_write_bcd_representation(vector[i] % 10000, nend - 3);
bc_write_bcd_representation(vector->values[i] / 10000, nend - 7);
bc_write_bcd_representation(vector->values[i] % 10000, nend - 3);
nend -= 8;
#endif
i++;
Expand All @@ -104,7 +104,7 @@ static inline void bc_convert_vector_to_char(const BC_VECTOR *vector, char *nptr
* The last digit may carry over.
* Also need to fill it to the end with zeros, so loop until the end of the string.
*/
BC_VECTOR last = vector[i];
BC_VECTOR last = vector->values[i];
while (nend >= nptr) {
*nend-- = last % BASE;
last /= BASE;
Expand Down
117 changes: 60 additions & 57 deletions ext/bcmath/libbcmath/src/div.c
Original file line number Diff line number Diff line change
Expand Up @@ -42,25 +42,24 @@
* This is because the algorithm can be simplified.
* This function is therefore an optimized version of bc_standard_div().
*/
static inline void bc_fast_div(
BC_VECTOR *numerator_vectors, size_t numerator_arr_size, BC_VECTOR divisor_vector, BC_VECTOR *quot_vectors, size_t quot_arr_size)
static inline void bc_fast_div(bc_vector_arr *numerator_vectors, BC_VECTOR divisor_vector, bc_vector_arr *quot_vectors)
{
size_t numerator_top_index = numerator_arr_size - 1;
size_t quot_top_index = quot_arr_size - 1;
size_t numerator_top_index = numerator_vectors->size - 1;
size_t quot_top_index = quot_vectors->size - 1;
for (size_t i = 0; i < quot_top_index; i++) {
if (numerator_vectors[numerator_top_index - i] < divisor_vector) {
quot_vectors[quot_top_index - i] = 0;
/* numerator_vectors[numerator_top_index - i] < divisor_vector, so there will be no overflow. */
numerator_vectors[numerator_top_index - i - 1] += numerator_vectors[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM;
if (numerator_vectors->values[numerator_top_index - i] < divisor_vector) {
quot_vectors->values[quot_top_index - i] = 0;
/* numerator_vectors->values[numerator_top_index - i] < divisor_vector, so there will be no overflow. */
numerator_vectors->values[numerator_top_index - i - 1] += numerator_vectors->values[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM;
continue;
}
quot_vectors[quot_top_index - i] = numerator_vectors[numerator_top_index - i] / divisor_vector;
numerator_vectors[numerator_top_index - i] -= divisor_vector * quot_vectors[quot_top_index - i];
numerator_vectors[numerator_top_index - i - 1] += numerator_vectors[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM;
numerator_vectors[numerator_top_index - i] = 0;
quot_vectors->values[quot_top_index - i] = numerator_vectors->values[numerator_top_index - i] / divisor_vector;
numerator_vectors->values[numerator_top_index - i] -= divisor_vector * quot_vectors->values[quot_top_index - i];
numerator_vectors->values[numerator_top_index - i - 1] += numerator_vectors->values[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM;
numerator_vectors->values[numerator_top_index - i] = 0;
}
/* last */
quot_vectors[0] = numerator_vectors[0] / divisor_vector;
quot_vectors->values[0] = numerator_vectors->values[0] / divisor_vector;
}

/*
Expand All @@ -69,13 +68,11 @@ static inline void bc_fast_div(
* https://en.wikipedia.org/wiki/Division_algorithm#Restoring_division
*/
static inline void bc_standard_div(
BC_VECTOR *numerator_vectors, size_t numerator_arr_size,
const BC_VECTOR *divisor_vectors, size_t divisor_arr_size, size_t divisor_len,
BC_VECTOR *quot_vectors, size_t quot_arr_size
) {
size_t numerator_top_index = numerator_arr_size - 1;
size_t divisor_top_index = divisor_arr_size - 1;
size_t quot_top_index = quot_arr_size - 1;
bc_vector_arr *numerator_vectors, const bc_vector_arr *divisor_vectors, size_t divisor_len, bc_vector_arr *quot_vectors)
{
size_t numerator_top_index = numerator_vectors->size - 1;
size_t divisor_top_index = divisor_vectors->size - 1;
size_t quot_top_index = quot_vectors->size - 1;

BC_VECTOR div_carry = 0;

Expand Down Expand Up @@ -172,9 +169,11 @@ static inline void bc_standard_div(

size_t high_part_shift = BC_POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1];
size_t low_part_shift = BC_POW_10_LUT[divisor_top_digits - 1];
BC_VECTOR divisor_high_part = divisor_vectors[divisor_top_index] * high_part_shift + divisor_vectors[divisor_top_index - 1] / low_part_shift;
for (size_t i = 0; i < quot_arr_size; i++) {
BC_VECTOR numerator_high_part = numerator_vectors[numerator_top_index - i] * high_part_shift + numerator_vectors[numerator_top_index - i - 1] / low_part_shift;
BC_VECTOR divisor_high_part = divisor_vectors->values[divisor_top_index] * high_part_shift
+ divisor_vectors->values[divisor_top_index - 1] / low_part_shift;
for (size_t i = 0; i < quot_vectors->size; i++) {
BC_VECTOR numerator_high_part = numerator_vectors->values[numerator_top_index - i] * high_part_shift
+ numerator_vectors->values[numerator_top_index - i - 1] / low_part_shift;

/*
* Determine the temporary quotient.
Expand All @@ -187,35 +186,35 @@ static inline void bc_standard_div(

/* numerator_high_part < divisor_high_part => quotient == 0 */
if (numerator_high_part < divisor_high_part) {
quot_vectors[quot_top_index - i] = 0;
div_carry = numerator_vectors[numerator_top_index - i];
numerator_vectors[numerator_top_index - i] = 0;
quot_vectors->values[quot_top_index - i] = 0;
div_carry = numerator_vectors->values[numerator_top_index - i];
numerator_vectors->values[numerator_top_index - i] = 0;
continue;
}

BC_VECTOR quot_guess = numerator_high_part / divisor_high_part;

/* For sub, add the remainder to the high-order digit */
numerator_vectors[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM;
numerator_vectors->values[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM;

/*
* It is impossible for the temporary quotient to underestimate the true quotient.
* Therefore a temporary quotient of 0 implies the true quotient is also 0.
*/
if (quot_guess == 0) {
quot_vectors[quot_top_index - i] = 0;
div_carry = numerator_vectors[numerator_top_index - i];
numerator_vectors[numerator_top_index - i] = 0;
quot_vectors->values[quot_top_index - i] = 0;
div_carry = numerator_vectors->values[numerator_top_index - i];
numerator_vectors->values[numerator_top_index - i] = 0;
continue;
}

/* sub */
BC_VECTOR sub;
BC_VECTOR borrow = 0;
BC_VECTOR *numerator_calc_bottom = numerator_vectors + numerator_arr_size - divisor_arr_size - i;
BC_VECTOR *numerator_calc_bottom = numerator_vectors->values + numerator_vectors->size - divisor_vectors->size - i;
size_t j;
for (j = 0; j < divisor_arr_size - 1; j++) {
sub = divisor_vectors[j] * quot_guess + borrow;
for (j = 0; j < divisor_vectors->size - 1; j++) {
sub = divisor_vectors->values[j] * quot_guess + borrow;
BC_VECTOR sub_low = sub % BC_VECTOR_BOUNDARY_NUM;
borrow = sub / BC_VECTOR_BOUNDARY_NUM;

Expand All @@ -227,26 +226,26 @@ static inline void bc_standard_div(
}
}
/* last digit sub */
sub = divisor_vectors[j] * quot_guess + borrow;
sub = divisor_vectors->values[j] * quot_guess + borrow;
bool neg_remainder = numerator_calc_bottom[j] < sub;
numerator_calc_bottom[j] -= sub;

/* If the temporary quotient is too large, add back divisor */
BC_VECTOR carry = 0;
if (neg_remainder) {
quot_guess--;
for (j = 0; j < divisor_arr_size - 1; j++) {
numerator_calc_bottom[j] += divisor_vectors[j] + carry;
for (j = 0; j < divisor_vectors->size - 1; j++) {
numerator_calc_bottom[j] += divisor_vectors->values[j] + carry;
carry = numerator_calc_bottom[j] / BC_VECTOR_BOUNDARY_NUM;
numerator_calc_bottom[j] %= BC_VECTOR_BOUNDARY_NUM;
}
/* last add */
numerator_calc_bottom[j] += divisor_vectors[j] + carry;
numerator_calc_bottom[j] += divisor_vectors->values[j] + carry;
}

quot_vectors[quot_top_index - i] = quot_guess;
div_carry = numerator_vectors[numerator_top_index - i];
numerator_vectors[numerator_top_index - i] = 0;
quot_vectors->values[quot_top_index - i] = quot_guess;
div_carry = numerator_vectors->values[numerator_top_index - i];
numerator_vectors->values[numerator_top_index - i] = 0;
}
}

Expand All @@ -255,48 +254,52 @@ static void bc_do_div(
const char *divisor, size_t divisor_size,
bc_num *quot, size_t quot_size
) {
size_t numerator_arr_size = BC_ARR_SIZE_FROM_LEN(numerator_size);
size_t divisor_arr_size = BC_ARR_SIZE_FROM_LEN(divisor_size);
size_t quot_arr_size = numerator_arr_size - divisor_arr_size + 1;
size_t quot_real_arr_size = MIN(quot_arr_size, BC_ARR_SIZE_FROM_LEN(quot_size));
bc_vector_arr numerator_vectors;
bc_vector_arr divisor_vectors;
bc_vector_arr quot_vectors;

numerator_vectors.size = BC_ARR_SIZE_FROM_LEN(numerator_size);
divisor_vectors.size = BC_ARR_SIZE_FROM_LEN(divisor_size);
quot_vectors.size = numerator_vectors.size - divisor_vectors.size + 1;
size_t quot_real_arr_size = MIN(quot_vectors.size, BC_ARR_SIZE_FROM_LEN(quot_size));

BC_VECTOR stack_vectors[BC_STACK_VECTOR_SIZE];
size_t allocation_arr_size = numerator_arr_size + divisor_arr_size + quot_arr_size;
size_t allocation_arr_size = numerator_vectors.size + divisor_vectors.size + quot_vectors.size;

BC_VECTOR *numerator_vectors;
if (allocation_arr_size <= BC_STACK_VECTOR_SIZE) {
numerator_vectors = stack_vectors;
numerator_vectors.values = stack_vectors;
} else {
numerator_vectors = safe_emalloc(allocation_arr_size, sizeof(BC_VECTOR), 0);
numerator_vectors.values = safe_emalloc(allocation_arr_size, sizeof(BC_VECTOR), 0);
}
BC_VECTOR *divisor_vectors = numerator_vectors + numerator_arr_size;
BC_VECTOR *quot_vectors = divisor_vectors + divisor_arr_size;
divisor_vectors.values = numerator_vectors.values + numerator_vectors.size;
quot_vectors.values = divisor_vectors.values + divisor_vectors.size;

size_t numerator_extension = numerator_size > numerator_readable_size ? numerator_size - numerator_readable_size : 0;

/* Bulk convert numerator and divisor to vectors */
size_t numerator_use_size = numerator_size - numerator_extension;
const char *numerator_end = numerator + numerator_use_size - 1;
bc_convert_to_vector_with_zero_pad(numerator_vectors, numerator_end, numerator_use_size, numerator_extension);
bc_convert_to_vector_with_zero_pad(numerator_vectors.values, numerator_end, numerator_use_size, numerator_extension);

const char *divisor_end = divisor + divisor_size - 1;
bc_convert_to_vector(divisor_vectors, divisor_end, divisor_size);
bc_convert_to_vector(divisor_vectors.values, divisor_end, divisor_size);

/* Do the division */
if (divisor_arr_size == 1) {
bc_fast_div(numerator_vectors, numerator_arr_size, divisor_vectors[0], quot_vectors, quot_arr_size);
if (divisor_vectors.size == 1) {
bc_fast_div(&numerator_vectors, divisor_vectors.values[0], &quot_vectors);
} else {
bc_standard_div(numerator_vectors, numerator_arr_size, divisor_vectors, divisor_arr_size, divisor_size, quot_vectors, quot_arr_size);
bc_standard_div(&numerator_vectors, &divisor_vectors, divisor_size, &quot_vectors);
}

/* Convert to bc_num */
char *qptr = (*quot)->n_value;
char *qend = qptr + (*quot)->n_len + (*quot)->n_scale - 1;

bc_convert_vector_to_char(quot_vectors, qptr, qend, quot_real_arr_size);
quot_vectors.size = quot_real_arr_size;
bc_convert_vector_to_char(&quot_vectors, qptr, qend);

if (allocation_arr_size > BC_STACK_VECTOR_SIZE) {
efree(numerator_vectors);
efree(numerator_vectors.values);
}
}

Expand Down
9 changes: 6 additions & 3 deletions ext/bcmath/libbcmath/src/private.h
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,11 @@
# define BC_VECTOR_BOUNDARY_NUM (BC_VECTOR) 10000
#endif

typedef struct bc_vector_arr {
BC_VECTOR *values;
size_t size;
} bc_vector_arr;

#ifdef WORDS_BIGENDIAN
# define BC_LITTLE_ENDIAN 0
#else
Expand All @@ -84,9 +89,7 @@ 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_multiply_vector(const bc_vector_arr *n1_vector, const bc_vector_arr *n2_vector, bc_vector_arr *prod_vector);
void _bc_rm_leading_zeros (bc_num num);

#endif
57 changes: 30 additions & 27 deletions ext/bcmath/libbcmath/src/raise.c
Original file line number Diff line number Diff line change
Expand Up @@ -36,28 +36,27 @@
#include <stdbool.h>
#include <stddef.h>

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)
static inline size_t bc_multiply_vector_ex(bc_vector_arr *n1_vector, bc_vector_arr *n2_vector, bc_vector_arr *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);
result_vector->size = n1_vector->size + n2_vector->size;
bc_multiply_vector(n1_vector, n2_vector, result_vector);

/* Eliminate extra zeros because they increase the number of calculations. */
while ((*result_vector)[result_arr_size - 1] == 0) {
result_arr_size--;
while (result_vector->values[result_vector->size - 1] == 0) {
result_vector->size--;
}

/* Swap n1_vector and result_vector. */
BC_VECTOR *tmp = *n1_vector;
*n1_vector = *result_vector;
*result_vector = tmp;
/* Swap n1_vector->values and result_vector->values. */
BC_VECTOR *tmp = n1_vector->values;
n1_vector->values = result_vector->values;
result_vector->values = tmp;

return result_arr_size;
return result_vector->size;
}

static inline size_t bc_square_vector_ex(BC_VECTOR **base_vector, size_t base_arr_size, BC_VECTOR **result_vector)
static inline size_t bc_square_vector_ex(bc_vector_arr *base_vector, bc_vector_arr *result_vector)
{
return bc_multiply_vector_ex(base_vector, base_arr_size, *base_vector, base_arr_size, result_vector);
return bc_multiply_vector_ex(base_vector, base_vector, result_vector);
}

/* Use "exponentiation by squaring". This is the fast path when the results are small. */
Expand Down Expand Up @@ -107,43 +106,47 @@ static bc_num bc_standard_raise(
base_len--;
}

size_t base_arr_size = BC_ARR_SIZE_FROM_LEN(base_len);
bc_vector_arr base_vector;
bc_vector_arr power_vector;
bc_vector_arr tmp_result_vector;

base_vector.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;
size_t max_power_arr_size = base_vector.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;
base_vector.values = buf;
power_vector.values = base_vector.values + max_power_arr_size;
tmp_result_vector.values = power_vector.values + max_power_arr_size;

/* Convert to BC_VECTOR[] */
bc_convert_to_vector(base_vector, base_end, base_len);
bc_convert_to_vector(base_vector.values, base_end, base_len);

while ((exponent & 1) == 0) {
base_arr_size = bc_square_vector_ex(&base_vector, base_arr_size, &tmp_result_vector);
base_vector.size = bc_square_vector_ex(&base_vector, &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];
power_vector.size = base_vector.size;
for (size_t i = 0; i < base_vector.size; i++) {
power_vector.values[i] = base_vector.values[i];
}
exponent >>= 1;

while (exponent > 0) {
base_arr_size = bc_square_vector_ex(&base_vector, base_arr_size, &tmp_result_vector);
base_vector.size = bc_square_vector_ex(&base_vector, &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);
power_vector.size = bc_multiply_vector_ex(&power_vector, &base_vector, &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;
size_t power_full_len = power_vector.size * BC_VECTOR_SIZE;
if (power_full_len > power_scale) {
power_len = power_full_len - power_scale;
} else {
Expand All @@ -160,7 +163,7 @@ static bc_num bc_standard_raise(
memset(pptr, 0, power_leading_zeros);
pptr += power_leading_zeros;

bc_convert_vector_to_char(power_vector, pptr, pend, power_arr_size);
bc_convert_vector_to_char(&power_vector, pptr, pend);

efree(buf);

Expand Down
Loading
Loading