From 5992fefe221fdfa9b92a843d597e27f0e048b76a Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Wed, 17 Jul 2024 21:53:23 +0900 Subject: [PATCH 01/34] Used BC_VECTOR for division to speed up calculations Co-authored-by: Niels Dossche <7771979+nielsdos@users.noreply.github.com> --- ext/bcmath/libbcmath/src/bcmath.h | 2 +- ext/bcmath/libbcmath/src/div.c | 549 +++++++++++++++++++++--------- 2 files changed, 380 insertions(+), 171 deletions(-) diff --git a/ext/bcmath/libbcmath/src/bcmath.h b/ext/bcmath/libbcmath/src/bcmath.h index 3ca5db6d408b1..75a17604f379f 100644 --- a/ext/bcmath/libbcmath/src/bcmath.h +++ b/ext/bcmath/libbcmath/src/bcmath.h @@ -146,7 +146,7 @@ bc_num bc_multiply(bc_num n1, bc_num n2, size_t scale); *(result) = mul_ex; \ } while (0) -bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int 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/div.c b/ext/bcmath/libbcmath/src/div.c index 5d3f0439cc0b7..1638520affaa6 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -31,219 +31,428 @@ #include "bcmath.h" #include "private.h" +#include "convert.h" #include #include #include #include "zend_alloc.h" +static const BC_VECTOR POW_10_LUT[9] = { + 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000 +}; -/* Some utility routines for the divide: First a one digit multiply. - NUM (with SIZE digits) is multiplied by DIGIT and the result is - placed into RESULT. It is written so that NUM and RESULT can be - the same pointers. */ - -static void _one_mult(unsigned char *num, size_t size, int digit, unsigned char *result) +/* + * Use this function when the number of array elements in n2 is 1, that is, when n2 is not divided into multiple array elements. + * The algorithm can be simplified if n2 is not divided, and this function is an optimization of bc_standard_div. + */ +static inline void bc_fast_div(BC_VECTOR *n1_vector, size_t n1_arr_size, BC_VECTOR n2_vector, BC_VECTOR *quot_vector, size_t quot_arr_size) { - size_t carry, value; - unsigned char *nptr, *rptr; + size_t n1_top_index = n1_arr_size - 1; + size_t quot_top_index = quot_arr_size - 1; + for (size_t i = 0; i < quot_arr_size - 1; i++) { + if (n1_vector[n1_top_index - i] < n2_vector) { + quot_vector[quot_top_index - i] = 0; + /* n1_vector[n1_top_index - i] is always a smaller number than n2_vector. Therefore, there will be no overflow. */ + n1_vector[n1_top_index - i - 1] += n1_vector[n1_top_index - i] * BC_VECTOR_BOUNDARY_NUM; + continue; + } + quot_vector[quot_top_index - i] = n1_vector[n1_top_index - i] / n2_vector; + n1_vector[n1_top_index - i] -= n2_vector * quot_vector[quot_top_index - i]; + n1_vector[n1_top_index - i - 1] += n1_vector[n1_top_index - i] * BC_VECTOR_BOUNDARY_NUM; + n1_vector[n1_top_index - i] = 0; + } + /* last */ + quot_vector[0] = n1_vector[0] / n2_vector; +} - if (digit == 0) { - memset(result, 0, size); - } else { - if (digit == 1) { - memcpy(result, num, size); - } else { - /* Initialize */ - nptr = (unsigned char *) (num + size - 1); - rptr = (unsigned char *) (result + size - 1); - carry = 0; - - while (size-- > 0) { - value = *nptr-- * digit + carry; - *rptr-- = value % BASE; - carry = value / BASE; - } +/* + * Used when the number of elements in the n2 vector is 2 or more. + */ +static inline void bc_standard_div( + BC_VECTOR *n1_vector, size_t n1_arr_size, + BC_VECTOR *n2_vector, size_t n2_arr_size, size_t n2_len, + BC_VECTOR *quot_vector, size_t quot_arr_size +) { + size_t i, j; + size_t n1_top_index = n1_arr_size - 1; + size_t n2_top_index = n2_arr_size - 1; + size_t quot_top_index = quot_arr_size - 1; + + BC_VECTOR div_carry = 0; + + /* + * Consider the error E between the actual quotient and the temporary quotient. + * For example, the quotient of 240 / 121 is 1, but if calculate the quotient only using the high-order digits (24 / 12), + * it will be 2. In this case, 2 is the temporary quotient. Here, the error E is 1. + * Also note that for example 2400000 / 120, there will be 5 divisions. + * + * Another example: Calculating 99900000000 / 10009999 with radix 10000. + * The actual quotient is 9980, but if calculate it using only the high-order digits: + * 999 / 1000 = 0, Numbers that are not divisible are carried forward to the next division. + * 9990000 / 1000 = 9990 (Here, the temporary quotient is 9990 and the error E is 10.) + * + * If calculate the temporary quotient using only one array of n1 and n2, the error E can be larger than 1. + * In other words, in the restoring division, the count of additions for restore increases significantly. + * + * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to one + * at most, adjust the number of high-order digits used to calculate the temporary quotient as follows. + * - Adjust the number of digits of n2 used in the calculation to BC_VECTOR_SIZE + 1 digit. The missing digits are + * filled in from the next array element. + * - Add digits to n1 in the same way as the number of digits adjusted by n2. + * + * e.g. n1 = 123456780000, n2 = 123456789, radix = 10000 + * n1_arr = [0, 5678, 1234], n2_arr = [6789, 2345, 1] + * n1_top = 1234, n2_top = 1 + * n2_top_tmp = 12345 (+4 digits), n1_top_tmp = 12345678 (+4 digits) + * tmp_quot = n1_top_tmp / n2_top_tmp = 1000, tmp_rem = -9000 + * tmp_quot is too large, so tmp_quot--, tmp_rem += n2 (restoring) + * quot = 999, rem = 123447789 + * + * Details: + * Suppose that when calculating the temporary quotient, only the high-order elements of the BC_VECTOR array are used. + * At this time, n1 and n2 can be considered to be decomposed as follows. (Here, B = 10^b, any b ∈ ℕ) + * n1 = n1_high * B^k + n1_low + * n2 = n2_high * B^k + n2_low + * + * The error E can be expressed by the following formula. + * + * E = (n1_high * B^k) / (n2_high * B^k) - (n1_high * B^k + n1_low) / (n2_high * B^k + n2_low) + * E = n1_high / n2_high - (n1_high * B^k + n1_low) / (n2_high * B^k + n2_low) + * E = (n1_high * (n2_high * B^k + n2_low) - (n1_high * B^k + n1_low) * n2_high) / (n2_high * (n2_high * B^k + n2_low)) + * E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) + * + * Find the error MAX_E when the error E is maximum. + * First, n1_high, which only exists in the numerator, uses its maximum value. Considering carry-back, + * n1_high can be expressed as follows. + * n1_high = n2_high * B + * And n1_low is also present only in the numerator, but since it is a negative number, use the largest possible value, i.e. 0. + * n1_low = 0 + * + * MAX_E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) + * MAX_E = (n2_high * B * n2_low) / (n2_high * (n2_high * B^k + n2_low)) + * + * n2_low uses the maximum value. + * n2_low = B^k - 1 + * MAX_E = (n2_high * B * n2_low) / (n2_high * (n2_high * B^k + n2_low)) + * MAX_E = (n2_high * B * (B^k - 1)) / (n2_high * (n2_high * B^k + B^k - 1)) + * + * n2_high uses the minimum value, but want to see how the number of digits affects the error, so represent + * the minimum value as: + * n2_high = 10^x (any x ∈ ℕ) + * Since B = 10^b, the formula becomes: + * MAX_E = (n2_high * B * (B^k - 1)) / (n2_high * (n2_high * B^k + B^k - 1)) + * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) + * + * Now let's make an approximation. Remove -1 from the numerator. Approximate the numerator to be + * large and the denominator to be small, such that MAX_E is less than this expression. + * Also, the denominator "(10^b)^k - 1" is never a negative value, so delete it. + * + * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) + * MAX_E < (10^x * 10^b * (10^b)^k) / (10^x * 10^x * (10^b)^k) + * MAX_E < 10^b / 10^x + * MAX_E < 10^(b - x) + * + * Therefore, found that if the number of digits in radix B and n2_high are equal, the error will be less + * than 1 regardless of the number of digits in the value of k. + * + * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in divider to be BC_VECTOR_SIZE + 1. + */ + size_t n2_top_digits = n2_len % BC_VECTOR_SIZE; + if (n2_top_digits == 0) { + n2_top_digits = BC_VECTOR_SIZE; + } + size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - n2_top_digits + 1]; + size_t low_part_shift = POW_10_LUT[n2_top_digits - 1]; + BC_VECTOR divider = n2_vector[n2_top_index] * high_part_shift + n2_vector[n2_top_index - 1] / low_part_shift; + for (i = 0; i < quot_arr_size; i++) { + BC_VECTOR divided = n1_vector[n1_top_index - i] * high_part_shift + n1_vector[n1_top_index - i - 1] / low_part_shift; + + /* If it is clear that n2 is greater in this loop, then the quotient is 0. */ + if (div_carry == 0 && divided < divider) { + quot_vector[quot_top_index - i] = 0; + div_carry = n1_vector[n1_top_index - i]; + n1_vector[n1_top_index - i] = 0; + continue; + } - if (carry != 0) { - *rptr = carry; + /* + * Determine the temporary quotient. + * "divided" is n1_high in the previous example. The maximum value of n1_high is n2_high * B, + * so here it is divider * BC_VECTOR_BOUNDARY_NUM. + * The number of digits in divider is BC_VECTOR_SIZE + 1, so even if divider is at its maximum value, + * it will never overflow here. + */ + divided += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; + BC_VECTOR quot_guess = divided / divider; + + /* For sub, add the remainder to the high-order digit */ + n1_vector[n1_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; + + /* + * The temporary quotient can only have errors in the "too large" direction. + * So if the temporary quotient is 0 here, the quotient is 0. + */ + if (quot_guess == 0) { + quot_vector[quot_top_index - i] = 0; + div_carry = n1_vector[n1_top_index - i]; + n1_vector[n1_top_index - i] = 0; + continue; + } + + /* sub */ + BC_VECTOR sub; + BC_VECTOR borrow = 0; + BC_VECTOR *n1_calc_bottom = n1_vector + n1_arr_size - n2_arr_size - i; + for (j = 0; j < n2_arr_size - 1; j++) { + sub = n2_vector[j] * quot_guess + borrow; + BC_VECTOR sub_low = sub % BC_VECTOR_BOUNDARY_NUM; + borrow = sub / BC_VECTOR_BOUNDARY_NUM; + + if (n1_calc_bottom[j] >= sub_low) { + n1_calc_bottom[j] -= sub_low; + } else { + n1_calc_bottom[j] += BC_VECTOR_BOUNDARY_NUM - sub_low; + borrow++; } } + /* last digit sub */ + sub = n2_vector[j] * quot_guess + borrow; + bool neg_remainder = n1_calc_bottom[j] < sub; + n1_calc_bottom[j] -= sub; + + /* If the temporary quotient is too large, add back n2 */ + BC_VECTOR carry = 0; + if (neg_remainder) { + quot_guess--; + for (j = 0; j < n2_arr_size - 1; j++) { + n1_calc_bottom[j] += n2_vector[j] + carry; + carry = n1_calc_bottom[j] / BC_VECTOR_BOUNDARY_NUM; + n1_calc_bottom[j] %= BC_VECTOR_BOUNDARY_NUM; + } + /* last add */ + n1_calc_bottom[j] += n2_vector[j] + carry; + } + + quot_vector[quot_top_index - i] = quot_guess; + div_carry = n1_vector[n1_top_index - i]; + n1_vector[n1_top_index - i] = 0; } } +static void bc_do_div(char *n1, size_t n1_readable_len, size_t n1_bottom_extension, char *n2, size_t n2_len, bc_num *quot, size_t quot_len) +{ + size_t i; + + size_t n2_arr_size = (n2_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; + size_t n1_arr_size = (n1_readable_len + n1_bottom_extension + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; + size_t quot_arr_size = n1_arr_size - n2_arr_size + 1; + size_t quot_real_arr_size = MIN(quot_arr_size, (quot_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE); + + BC_VECTOR *n1_vector = safe_emalloc(n1_arr_size + n2_arr_size + quot_arr_size, sizeof(BC_VECTOR), 0); + BC_VECTOR *n2_vector = n1_vector + n1_arr_size; + BC_VECTOR *quot_vector = n2_vector + n2_arr_size; + + /* Fill with zeros and convert as many vector elements as needed */ + size_t n1_vector_count = 0; + while (n1_bottom_extension >= BC_VECTOR_SIZE) { + n1_vector[n1_vector_count] = 0; + n1_bottom_extension -= BC_VECTOR_SIZE; + n1_vector_count++; + } + + size_t n1_bottom_read_len = BC_VECTOR_SIZE - n1_bottom_extension; + + size_t base; + size_t n1_read = 0; + if (n1_bottom_read_len < BC_VECTOR_SIZE) { + n1_read = MIN(n1_bottom_read_len, n1_readable_len); + base = POW_10_LUT[n1_bottom_extension]; + n1_vector[n1_vector_count] = 0; + for (i = 0; i < n1_read; i++) { + n1_vector[n1_vector_count] += *n1 * base; + base *= BASE; + n1--; + } + n1_vector_count++; + } + + /* Bulk convert n1 and n2 to vectors */ + if (n1_readable_len > n1_read) { + bc_convert_to_vector(n1_vector + n1_vector_count, n1, n1_readable_len - n1_read); + } + bc_convert_to_vector(n2_vector, n2, n2_len); + + /* Do the divide */ + if (n2_arr_size == 1) { + bc_fast_div(n1_vector, n1_arr_size, n2_vector[0], quot_vector, quot_arr_size); + } else { + bc_standard_div(n1_vector, n1_arr_size, n2_vector, n2_arr_size, n2_len, quot_vector, quot_arr_size); + } + + /* Convert to bc_num */ + char *qptr = (*quot)->n_value; + char *qend = qptr + quot_len - 1; + + for (i = 0; i < quot_real_arr_size - 1; i++) { +#if BC_VECTOR_SIZE == 4 + bc_write_bcd_representation(quot_vector[i], qend - 3); + qend -= 4; +#else + bc_write_bcd_representation(quot_vector[i] / 10000, qend - 7); + bc_write_bcd_representation(quot_vector[i] % 10000, qend - 3); + qend -= 8; +#endif + } -/* The full division routine. This computes N1 / N2. It returns - true if the division is ok and the result is in QUOT. The number of - digits after the decimal point is SCALE. It returns false if division - by zero is tried. The algorithm is found in Knuth Vol 2. p237. */ -bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) + while (qend >= qptr) { + *qend-- = quot_vector[i] % BASE; + quot_vector[i] /= BASE; + } + + efree(n1_vector); +} + +bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, size_t scale) { - bc_num qval; - unsigned char *num1, *num2; - unsigned char *ptr1, *ptr2, *n2ptr, *qptr; - int scale1, val; - unsigned int len1, len2, scale2, qdigits, extra, count; - unsigned int qdig, qguess, borrow, carry; - unsigned char *mval; - unsigned int norm; - bool zero; - - /* Test for divide by zero. */ + /* divide by zero */ if (bc_is_zero(n2)) { return false; } - /* Test for divide by 1. If it is we must truncate. */ - if (n2->n_scale == 0 && n2->n_len == 1 && *n2->n_value == 1) { - qval = bc_new_num (n1->n_len, scale); - qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); - memcpy(qval->n_value, n1->n_value, n1->n_len + MIN(n1->n_scale, scale)); - bc_free_num (quot); - *quot = qval; + bc_free_num(quot); + + /* If n1 is zero, the quotient is always zero. */ + if (bc_is_zero(n1)) { + *quot = bc_copy_num(BCG(_zero_)); + return true; } - /* Set up the divide. Move the decimal point on n1 by n2's scale. - Remember, zeros on the end of num2 are wasted effort for dividing. */ - scale2 = n2->n_scale; - n2ptr = (unsigned char *) n2->n_value + n2->n_len + scale2 - 1; - while ((scale2 > 0) && (*n2ptr == 0)) { - scale2--; - n2ptr--; + char *n1ptr = n1->n_value; + char *n1end = n1ptr + n1->n_len + n1->n_scale - 1; + size_t n1_len = n1->n_len; + size_t n1_scale = n1->n_scale; + + char *n2ptr = n2->n_value; + char *n2end = n2ptr + n2->n_len + n2->n_scale - 1; + size_t n2_len = n2->n_len; + size_t n2_scale = n2->n_scale; + size_t n2_int_right_zeros = 0; + + /* remove n2 trailing zeros */ + while (*n2end == 0 && n2_scale > 0) { + n2end--; + n2_scale--; + } + while (*n2end == 0) { + n2end--; + n2_int_right_zeros++; } - len1 = n1->n_len + scale2; - scale1 = n1->n_scale - scale2; - extra = MAX(scale - scale1, 0); + if (*n1ptr == 0 && n1_len == 1) { + n1ptr++; + n1_len = 0; + } - num1 = (unsigned char *) safe_emalloc(1, n1->n_len + n1->n_scale, extra + 2); - memset(num1, 0, n1->n_len + n1->n_scale + extra + 2); - memcpy(num1 + 1, n1->n_value, n1->n_len + n1->n_scale); + size_t n1_top_extension = 0; + size_t n1_bottom_extension = 0; + if (n2_scale > 0) { + /* + * e.g. n2_scale = 4 + * n2 = .0002, to be 2 or n2 = 200.001, to be 200001 + * n1 = .03, to be 300 or n1 = .000003, to be .03 + * n1 may become longer than the original data length due to the addition of + * trailing zeros in the integer part. + */ + n1_len += n2_scale; + n1_bottom_extension = n1_scale < n2_scale ? n2_scale - n1_scale : 0; + n1_scale = n1_scale > n2_scale ? n1_scale - n2_scale : 0; + n2_len += n2_scale; + n2_scale = 0; + } else if (n2_int_right_zeros > 0) { + /* + * e.g. n2_int_right_zeros = 4 + * n2 = 2000, to be 2 + * n1 = 30, to be .03 or n1 = 30000, to be 30 + * Also, n1 may become longer than the original data length due to the addition of + * leading zeros in the fractional part. + */ + n1_top_extension = n1_len < n2_int_right_zeros ? n2_int_right_zeros - n1_len : 0; + n1_len = n1_len > n2_int_right_zeros ? n1_len - n2_int_right_zeros : 0; + n1_scale += n2_int_right_zeros; + n2_len -= n2_int_right_zeros; + n2_scale = 0; + } - len2 = n2->n_len + scale2; - num2 = (unsigned char *) safe_emalloc(1, len2, 1); - memcpy(num2, n2->n_value, len2); - *(num2 + len2) = 0; - n2ptr = num2; + /* remove n1 leading zeros */ + while (*n1ptr == 0 && n1_len > 0) { + n1ptr++; + n1_len--; + } + /* remove n2 leading zeros */ while (*n2ptr == 0) { n2ptr++; - len2--; + n2_len--; } - /* Calculate the number of quotient digits. */ - if (len2 > len1 + scale) { - qdigits = scale + 1; - zero = true; - } else { - zero = false; - if (len2 > len1) { - /* One for the zero integer part. */ - qdigits = scale + 1; + /* Considering the scale specification, the quotient is always 0 if this condition is met */ + if (n2_len > n1_len + scale) { + *quot = bc_copy_num(BCG(_zero_)); + return true; + } + + /* set scale to n1 */ + if (n1_scale > scale) { + size_t scale_diff = n1_scale - scale; + if (n1_bottom_extension > scale_diff) { + n1_bottom_extension -= scale_diff; } else { - qdigits = len1 - len2 + scale + 1; + n1_bottom_extension = 0; + n1end -= scale_diff - n1_bottom_extension; } + } else { + n1_bottom_extension += scale - n1_scale; } + n1_scale = scale; - /* Allocate and zero the storage for the quotient. */ - qval = bc_new_num (qdigits - scale, scale); + /* Length of n1 data that can be read */ + size_t n1_readable_len = n1end - n1ptr + 1; - /* Allocate storage for the temporary storage mval. */ - mval = (unsigned char *) safe_emalloc(1, len2, 1); - - /* Now for the full divide algorithm. */ - if (!zero) { - /* Normalize */ - norm = 10 / ((int) *n2ptr + 1); - if (norm != 1) { - _one_mult(num1, len1 + scale1 + extra + 1, norm, num1); - _one_mult(n2ptr, len2, norm, n2ptr); + /* If n2 is 1 here, return the result of adjusting the decimal point position of n1. */ + if (n2_len == 1 && *n2ptr == 1) { + if (n1_len == 0) { + n1_len = 1; + n1_top_extension++; } + size_t quot_scale = n1_scale > n1_bottom_extension ? n1_scale - n1_bottom_extension : 0; + n1_bottom_extension = n1_scale < n1_bottom_extension ? n1_bottom_extension - n1_scale : 0; - /* Initialize divide loop. */ - qdig = 0; - if (len2 > len1) { - qptr = (unsigned char *) qval->n_value + len2 - len1; - } else { - qptr = (unsigned char *) qval->n_value; + *quot = bc_new_num_nonzeroed(n1_len, quot_scale); + char *qptr = (*quot)->n_value; + for (size_t i = 0; i < n1_top_extension; i++) { + *qptr++ = 0; } - - /* Loop */ - while (qdig <= len1 + scale - len2) { - /* Calculate the quotient digit guess. */ - if (*n2ptr == num1[qdig]) { - qguess = 9; - } else { - qguess = (num1[qdig] * 10 + num1[qdig + 1]) / *n2ptr; - } - - /* Test qguess. */ - if (n2ptr[1] * qguess > (num1[qdig] * 10 + num1[qdig + 1] - *n2ptr * qguess) * 10 + num1[qdig + 2]) { - qguess--; - /* And again. */ - if (n2ptr[1] * qguess > (num1[qdig] * 10 + num1[qdig + 1] - *n2ptr * qguess) * 10 + num1[qdig + 2]) { - qguess--; - } - } - - /* Multiply and subtract. */ - borrow = 0; - if (qguess != 0) { - *mval = 0; - _one_mult(n2ptr, len2, qguess, mval + 1); - ptr1 = (unsigned char *) num1 + qdig + len2; - ptr2 = (unsigned char *) mval + len2; - for (count = 0; count < len2 + 1; count++) { - val = (int) *ptr1 - (int) *ptr2-- - borrow; - if (val < 0) { - val += 10; - borrow = 1; - } else { - borrow = 0; - } - *ptr1-- = val; - } - } - - /* Test for negative result. */ - if (borrow == 1) { - qguess--; - ptr1 = (unsigned char *) num1 + qdig + len2; - ptr2 = (unsigned char *) n2ptr + len2 - 1; - carry = 0; - for (count = 0; count < len2; count++) { - val = (int) *ptr1 + (int) *ptr2-- + carry; - if (val > 9) { - val -= 10; - carry = 1; - } else { - carry = 0; - } - *ptr1-- = val; - } - if (carry == 1) { - *ptr1 = (*ptr1 + 1) % 10; - } - } - - /* We now know the quotient digit. */ - *qptr++ = qguess; - qdig++; + memcpy(qptr, n1ptr, n1_readable_len); + qptr += n1_readable_len; + for (size_t i = 0; i < n1_bottom_extension; i++) { + *qptr++ = 0; } + (*quot)->n_sign = n1->n_sign == n2->n_sign ? PLUS : MINUS; + return true; } - /* Clean up and return the number. */ - qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); - if (bc_is_zero(qval)) { - qval->n_sign = PLUS; + size_t quot_full_len; + if (n2_len > n1_len) { + *quot = bc_new_num_nonzeroed(1, scale); + quot_full_len = 1 + scale; + } else { + *quot = bc_new_num_nonzeroed(n1_len - n2_len + 1, scale); + quot_full_len = n1_len - n2_len + 1 + scale; } - _bc_rm_leading_zeros(qval); - bc_free_num(quot); - *quot = qval; - /* Clean up temporary storage. */ - efree(mval); - efree(num1); - efree(num2); + /* do divide */ + bc_do_div(n1end, n1_readable_len, n1_bottom_extension, n2end, n2_len, quot, quot_full_len); + (*quot)->n_sign = n1->n_sign == n2->n_sign ? PLUS : MINUS; + _bc_rm_leading_zeros(*quot); - /* Everything is OK. */ return true; } From 3e3750e3e211a7090f0fc5049f8ac5f9ed84c8d8 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Wed, 17 Jul 2024 21:57:44 +0900 Subject: [PATCH 02/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 1638520affaa6..16f8c3a6be394 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -94,8 +94,8 @@ static inline void bc_standard_div( * If calculate the temporary quotient using only one array of n1 and n2, the error E can be larger than 1. * In other words, in the restoring division, the count of additions for restore increases significantly. * - * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to one - * at most, adjust the number of high-order digits used to calculate the temporary quotient as follows. + * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to + * at most one, adjust the number of high-order digits used to calculate the temporary quotient as follows. * - Adjust the number of digits of n2 used in the calculation to BC_VECTOR_SIZE + 1 digit. The missing digits are * filled in from the next array element. * - Add digits to n1 in the same way as the number of digits adjusted by n2. From d34e9df04fbc54e37e5f010822fe368ddf608cbd Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Wed, 17 Jul 2024 22:50:03 +0900 Subject: [PATCH 03/34] Fixed comments --- ext/bcmath/libbcmath/src/div.c | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 16f8c3a6be394..64e832c5e846f 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -86,7 +86,7 @@ static inline void bc_standard_div( * it will be 2. In this case, 2 is the temporary quotient. Here, the error E is 1. * Also note that for example 2400000 / 120, there will be 5 divisions. * - * Another example: Calculating 99900000000 / 10009999 with radix 10000. + * Another example: Calculating 99900000000 / 10009999 with base 10000. * The actual quotient is 9980, but if calculate it using only the high-order digits: * 999 / 1000 = 0, Numbers that are not divisible are carried forward to the next division. * 9990000 / 1000 = 9990 (Here, the temporary quotient is 9990 and the error E is 10.) @@ -100,17 +100,25 @@ static inline void bc_standard_div( * filled in from the next array element. * - Add digits to n1 in the same way as the number of digits adjusted by n2. * - * e.g. n1 = 123456780000, n2 = 123456789, radix = 10000 - * n1_arr = [0, 5678, 1234], n2_arr = [6789, 2345, 1] - * n1_top = 1234, n2_top = 1 - * n2_top_tmp = 12345 (+4 digits), n1_top_tmp = 12345678 (+4 digits) - * tmp_quot = n1_top_tmp / n2_top_tmp = 1000, tmp_rem = -9000 + * e.g. + * n1 = 123456780000 + * n2 = 123456789 + * base = 10000 + * n1_arr = [0, 5678, 1234] + * n2_arr = [6789, 2345, 1] + * n1_top = 1234 + * n2_top = 1 + * n2_top_tmp = 12345 (+4 digits) + * n1_top_tmp = 12345678 (+4 digits) + * tmp_quot = n1_top_tmp / n2_top_tmp = 1000 + * tmp_rem = -9000 * tmp_quot is too large, so tmp_quot--, tmp_rem += n2 (restoring) - * quot = 999, rem = 123447789 + * quot = 999 + * rem = 123447789 * * Details: * Suppose that when calculating the temporary quotient, only the high-order elements of the BC_VECTOR array are used. - * At this time, n1 and n2 can be considered to be decomposed as follows. (Here, B = 10^b, any b ∈ ℕ) + * At this time, n1 and n2 can be considered to be decomposed as follows. (Here, B = 10^b, any b ∈ ℕ, any k ∈ ℕ) * n1 = n1_high * B^k + n1_low * n2 = n2_high * B^k + n2_low * @@ -125,7 +133,7 @@ static inline void bc_standard_div( * First, n1_high, which only exists in the numerator, uses its maximum value. Considering carry-back, * n1_high can be expressed as follows. * n1_high = n2_high * B - * And n1_low is also present only in the numerator, but since it is a negative number, use the largest possible value, i.e. 0. + * Also, n1_low is only present in the numerator, but since this is a subtraction, use the smallest possible value here, 0. * n1_low = 0 * * MAX_E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) @@ -152,7 +160,7 @@ static inline void bc_standard_div( * MAX_E < 10^b / 10^x * MAX_E < 10^(b - x) * - * Therefore, found that if the number of digits in radix B and n2_high are equal, the error will be less + * Therefore, found that if the number of digits in base B and n2_high are equal, the error will be less * than 1 regardless of the number of digits in the value of k. * * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in divider to be BC_VECTOR_SIZE + 1. From c56d0f3ccd2e348f5cfadfd02d0f2dc56842086c Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Wed, 17 Jul 2024 23:03:45 +0900 Subject: [PATCH 04/34] Correction of line break position --- ext/bcmath/libbcmath/src/div.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 64e832c5e846f..a06708536f6ac 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -130,8 +130,8 @@ static inline void bc_standard_div( * E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) * * Find the error MAX_E when the error E is maximum. - * First, n1_high, which only exists in the numerator, uses its maximum value. Considering carry-back, - * n1_high can be expressed as follows. + * First, n1_high, which only exists in the numerator, uses its maximum value. + * Considering carry-back, n1_high can be expressed as follows. * n1_high = n2_high * B * Also, n1_low is only present in the numerator, but since this is a subtraction, use the smallest possible value here, 0. * n1_low = 0 From 9177d4c16aec0e5fddc6cf9c5307a42f66294946 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 08:53:48 +0900 Subject: [PATCH 05/34] Changed n1 to numerator and n2 to divisor --- ext/bcmath/libbcmath/src/div.c | 167 +++++++++++++++++---------------- 1 file changed, 84 insertions(+), 83 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index a06708536f6ac..c53bc180ecdf8 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -42,40 +42,41 @@ static const BC_VECTOR POW_10_LUT[9] = { }; /* - * Use this function when the number of array elements in n2 is 1, that is, when n2 is not divided into multiple array elements. - * The algorithm can be simplified if n2 is not divided, and this function is an optimization of bc_standard_div. + * Use this function when the number of array elements in divisor is 1, that is, when divisor is not divided into multiple array elements. + * The algorithm can be simplified if divisor is not divided, and this function is an optimization of bc_standard_div. */ -static inline void bc_fast_div(BC_VECTOR *n1_vector, size_t n1_arr_size, BC_VECTOR n2_vector, BC_VECTOR *quot_vector, size_t quot_arr_size) +static inline void bc_fast_div( + BC_VECTOR *numerator_vector, size_t numerator_arr_size, BC_VECTOR divisor_vector, BC_VECTOR *quot_vector, size_t quot_arr_size) { - size_t n1_top_index = n1_arr_size - 1; + size_t numerator_top_index = numerator_arr_size - 1; size_t quot_top_index = quot_arr_size - 1; for (size_t i = 0; i < quot_arr_size - 1; i++) { - if (n1_vector[n1_top_index - i] < n2_vector) { + if (numerator_vector[numerator_top_index - i] < divisor_vector) { quot_vector[quot_top_index - i] = 0; - /* n1_vector[n1_top_index - i] is always a smaller number than n2_vector. Therefore, there will be no overflow. */ - n1_vector[n1_top_index - i - 1] += n1_vector[n1_top_index - i] * BC_VECTOR_BOUNDARY_NUM; + /* numerator_vector[numerator_top_index - i] is always a smaller number than divisor_vector. Therefore, there will be no overflow. */ + numerator_vector[numerator_top_index - i - 1] += numerator_vector[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM; continue; } - quot_vector[quot_top_index - i] = n1_vector[n1_top_index - i] / n2_vector; - n1_vector[n1_top_index - i] -= n2_vector * quot_vector[quot_top_index - i]; - n1_vector[n1_top_index - i - 1] += n1_vector[n1_top_index - i] * BC_VECTOR_BOUNDARY_NUM; - n1_vector[n1_top_index - i] = 0; + quot_vector[quot_top_index - i] = numerator_vector[numerator_top_index - i] / divisor_vector; + numerator_vector[numerator_top_index - i] -= divisor_vector * quot_vector[quot_top_index - i]; + numerator_vector[numerator_top_index - i - 1] += numerator_vector[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM; + numerator_vector[numerator_top_index - i] = 0; } /* last */ - quot_vector[0] = n1_vector[0] / n2_vector; + quot_vector[0] = numerator_vector[0] / divisor_vector; } /* - * Used when the number of elements in the n2 vector is 2 or more. + * Used when the number of elements in the divisor vector is 2 or more. */ static inline void bc_standard_div( - BC_VECTOR *n1_vector, size_t n1_arr_size, - BC_VECTOR *n2_vector, size_t n2_arr_size, size_t n2_len, + BC_VECTOR *numerator_vector, size_t numerator_arr_size, + BC_VECTOR *divisor_vector, size_t divisor_arr_size, size_t divisor_len, BC_VECTOR *quot_vector, size_t quot_arr_size ) { size_t i, j; - size_t n1_top_index = n1_arr_size - 1; - size_t n2_top_index = n2_arr_size - 1; + 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 div_carry = 0; @@ -91,64 +92,64 @@ static inline void bc_standard_div( * 999 / 1000 = 0, Numbers that are not divisible are carried forward to the next division. * 9990000 / 1000 = 9990 (Here, the temporary quotient is 9990 and the error E is 10.) * - * If calculate the temporary quotient using only one array of n1 and n2, the error E can be larger than 1. + * If calculate the temporary quotient using only one array of numerator and divisor, the error E can be larger than 1. * In other words, in the restoring division, the count of additions for restore increases significantly. * * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to * at most one, adjust the number of high-order digits used to calculate the temporary quotient as follows. - * - Adjust the number of digits of n2 used in the calculation to BC_VECTOR_SIZE + 1 digit. The missing digits are + * - Adjust the number of digits of divisor used in the calculation to BC_VECTOR_SIZE + 1 digit. The missing digits are * filled in from the next array element. - * - Add digits to n1 in the same way as the number of digits adjusted by n2. + * - Add digits to numerator in the same way as the number of digits adjusted by divisor. * * e.g. - * n1 = 123456780000 - * n2 = 123456789 + * numerator = 123456780000 + * divisor = 123456789 * base = 10000 - * n1_arr = [0, 5678, 1234] - * n2_arr = [6789, 2345, 1] - * n1_top = 1234 - * n2_top = 1 - * n2_top_tmp = 12345 (+4 digits) - * n1_top_tmp = 12345678 (+4 digits) - * tmp_quot = n1_top_tmp / n2_top_tmp = 1000 + * numerator_arr = [0, 5678, 1234] + * divisor_arr = [6789, 2345, 1] + * numerator_top = 1234 + * divisor_top = 1 + * divisor_top_tmp = 12345 (+4 digits) + * numerator_top_tmp = 12345678 (+4 digits) + * tmp_quot = numerator_top_tmp / divisor_top_tmp = 1000 * tmp_rem = -9000 - * tmp_quot is too large, so tmp_quot--, tmp_rem += n2 (restoring) + * tmp_quot is too large, so tmp_quot--, tmp_rem += divisor (restoring) * quot = 999 * rem = 123447789 * * Details: * Suppose that when calculating the temporary quotient, only the high-order elements of the BC_VECTOR array are used. - * At this time, n1 and n2 can be considered to be decomposed as follows. (Here, B = 10^b, any b ∈ ℕ, any k ∈ ℕ) - * n1 = n1_high * B^k + n1_low - * n2 = n2_high * B^k + n2_low + * At this time, numerator and divisor can be considered to be decomposed as follows. (Here, B = 10^b, any b ∈ ℕ, any k ∈ ℕ) + * numerator = numerator_high * B^k + numerator_low + * divisor = divisor_high * B^k + divisor_low * * The error E can be expressed by the following formula. * - * E = (n1_high * B^k) / (n2_high * B^k) - (n1_high * B^k + n1_low) / (n2_high * B^k + n2_low) - * E = n1_high / n2_high - (n1_high * B^k + n1_low) / (n2_high * B^k + n2_low) - * E = (n1_high * (n2_high * B^k + n2_low) - (n1_high * B^k + n1_low) * n2_high) / (n2_high * (n2_high * B^k + n2_low)) - * E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) + * E = (numerator_high * B^k) / (divisor_high * B^k) - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) + * E = numerator_high / divisor_high - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) + * E = (numerator_high * (divisor_high * B^k + divisor_low) - (numerator_high * B^k + numerator_low) * divisor_high) / (divisor_high * (divisor_high * B^k + divisor_low)) + * E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * * Find the error MAX_E when the error E is maximum. - * First, n1_high, which only exists in the numerator, uses its maximum value. - * Considering carry-back, n1_high can be expressed as follows. - * n1_high = n2_high * B - * Also, n1_low is only present in the numerator, but since this is a subtraction, use the smallest possible value here, 0. - * n1_low = 0 + * First, numerator_high, which only exists in the numerator, uses its maximum value. + * Considering carry-back, numerator_high can be expressed as follows. + * numerator_high = divisor_high * B + * Also, numerator_low is only present in the numerator, but since this is a subtraction, use the smallest possible value here, 0. + * numerator_low = 0 * - * MAX_E = (n1_high * n2_low - n2_high * n1_low) / (n2_high * (n2_high * B^k + n2_low)) - * MAX_E = (n2_high * B * n2_low) / (n2_high * (n2_high * B^k + n2_low)) + * MAX_E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * - * n2_low uses the maximum value. - * n2_low = B^k - 1 - * MAX_E = (n2_high * B * n2_low) / (n2_high * (n2_high * B^k + n2_low)) - * MAX_E = (n2_high * B * (B^k - 1)) / (n2_high * (n2_high * B^k + B^k - 1)) + * divisor_low uses the maximum value. + * divisor_low = B^k - 1 + * MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) * - * n2_high uses the minimum value, but want to see how the number of digits affects the error, so represent + * divisor_high uses the minimum value, but want to see how the number of digits affects the error, so represent * the minimum value as: - * n2_high = 10^x (any x ∈ ℕ) + * divisor_high = 10^x (any x ∈ ℕ) * Since B = 10^b, the formula becomes: - * MAX_E = (n2_high * B * (B^k - 1)) / (n2_high * (n2_high * B^k + B^k - 1)) + * MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) * * Now let's make an approximation. Remove -1 from the numerator. Approximate the numerator to be @@ -160,32 +161,32 @@ static inline void bc_standard_div( * MAX_E < 10^b / 10^x * MAX_E < 10^(b - x) * - * Therefore, found that if the number of digits in base B and n2_high are equal, the error will be less + * Therefore, found that if the number of digits in base B and divisor_high are equal, the error will be less * than 1 regardless of the number of digits in the value of k. * * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in divider to be BC_VECTOR_SIZE + 1. */ - size_t n2_top_digits = n2_len % BC_VECTOR_SIZE; - if (n2_top_digits == 0) { - n2_top_digits = BC_VECTOR_SIZE; + size_t divisor_top_digits = divisor_len % BC_VECTOR_SIZE; + if (divisor_top_digits == 0) { + divisor_top_digits = BC_VECTOR_SIZE; } - size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - n2_top_digits + 1]; - size_t low_part_shift = POW_10_LUT[n2_top_digits - 1]; - BC_VECTOR divider = n2_vector[n2_top_index] * high_part_shift + n2_vector[n2_top_index - 1] / low_part_shift; + size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; + size_t low_part_shift = POW_10_LUT[divisor_top_digits - 1]; + BC_VECTOR divider = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; for (i = 0; i < quot_arr_size; i++) { - BC_VECTOR divided = n1_vector[n1_top_index - i] * high_part_shift + n1_vector[n1_top_index - i - 1] / low_part_shift; + BC_VECTOR divided = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; - /* If it is clear that n2 is greater in this loop, then the quotient is 0. */ + /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ if (div_carry == 0 && divided < divider) { quot_vector[quot_top_index - i] = 0; - div_carry = n1_vector[n1_top_index - i]; - n1_vector[n1_top_index - i] = 0; + div_carry = numerator_vector[numerator_top_index - i]; + numerator_vector[numerator_top_index - i] = 0; continue; } /* * Determine the temporary quotient. - * "divided" is n1_high in the previous example. The maximum value of n1_high is n2_high * B, + * "divided" is numerator_high in the previous example. The maximum value of numerator_high is divisor_high * B, * so here it is divider * BC_VECTOR_BOUNDARY_NUM. * The number of digits in divider is BC_VECTOR_SIZE + 1, so even if divider is at its maximum value, * it will never overflow here. @@ -194,7 +195,7 @@ static inline void bc_standard_div( BC_VECTOR quot_guess = divided / divider; /* For sub, add the remainder to the high-order digit */ - n1_vector[n1_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; + numerator_vector[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; /* * The temporary quotient can only have errors in the "too large" direction. @@ -202,48 +203,48 @@ static inline void bc_standard_div( */ if (quot_guess == 0) { quot_vector[quot_top_index - i] = 0; - div_carry = n1_vector[n1_top_index - i]; - n1_vector[n1_top_index - i] = 0; + div_carry = numerator_vector[numerator_top_index - i]; + numerator_vector[numerator_top_index - i] = 0; continue; } /* sub */ BC_VECTOR sub; BC_VECTOR borrow = 0; - BC_VECTOR *n1_calc_bottom = n1_vector + n1_arr_size - n2_arr_size - i; - for (j = 0; j < n2_arr_size - 1; j++) { - sub = n2_vector[j] * quot_guess + borrow; + BC_VECTOR *numerator_calc_bottom = numerator_vector + numerator_arr_size - divisor_arr_size - i; + for (j = 0; j < divisor_arr_size - 1; j++) { + sub = divisor_vector[j] * quot_guess + borrow; BC_VECTOR sub_low = sub % BC_VECTOR_BOUNDARY_NUM; borrow = sub / BC_VECTOR_BOUNDARY_NUM; - if (n1_calc_bottom[j] >= sub_low) { - n1_calc_bottom[j] -= sub_low; + if (numerator_calc_bottom[j] >= sub_low) { + numerator_calc_bottom[j] -= sub_low; } else { - n1_calc_bottom[j] += BC_VECTOR_BOUNDARY_NUM - sub_low; + numerator_calc_bottom[j] += BC_VECTOR_BOUNDARY_NUM - sub_low; borrow++; } } /* last digit sub */ - sub = n2_vector[j] * quot_guess + borrow; - bool neg_remainder = n1_calc_bottom[j] < sub; - n1_calc_bottom[j] -= sub; + sub = divisor_vector[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 n2 */ + /* If the temporary quotient is too large, add back divisor */ BC_VECTOR carry = 0; if (neg_remainder) { quot_guess--; - for (j = 0; j < n2_arr_size - 1; j++) { - n1_calc_bottom[j] += n2_vector[j] + carry; - carry = n1_calc_bottom[j] / BC_VECTOR_BOUNDARY_NUM; - n1_calc_bottom[j] %= BC_VECTOR_BOUNDARY_NUM; + for (j = 0; j < divisor_arr_size - 1; j++) { + numerator_calc_bottom[j] += divisor_vector[j] + carry; + carry = numerator_calc_bottom[j] / BC_VECTOR_BOUNDARY_NUM; + numerator_calc_bottom[j] %= BC_VECTOR_BOUNDARY_NUM; } /* last add */ - n1_calc_bottom[j] += n2_vector[j] + carry; + numerator_calc_bottom[j] += divisor_vector[j] + carry; } quot_vector[quot_top_index - i] = quot_guess; - div_carry = n1_vector[n1_top_index - i]; - n1_vector[n1_top_index - i] = 0; + div_carry = numerator_vector[numerator_top_index - i]; + numerator_vector[numerator_top_index - i] = 0; } } From bf0f88df7609b9c506ec007b9b203901023ea3c9 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 08:57:10 +0900 Subject: [PATCH 06/34] Changed the naming of divided and divider --- ext/bcmath/libbcmath/src/div.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index c53bc180ecdf8..f9b4d1f1924d9 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -164,7 +164,7 @@ static inline void bc_standard_div( * Therefore, found that if the number of digits in base B and divisor_high are equal, the error will be less * than 1 regardless of the number of digits in the value of k. * - * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in divider to be BC_VECTOR_SIZE + 1. + * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in high part of divisor to be BC_VECTOR_SIZE + 1. */ size_t divisor_top_digits = divisor_len % BC_VECTOR_SIZE; if (divisor_top_digits == 0) { @@ -172,12 +172,12 @@ static inline void bc_standard_div( } size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = POW_10_LUT[divisor_top_digits - 1]; - BC_VECTOR divider = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; + BC_VECTOR tmp_divisor = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; for (i = 0; i < quot_arr_size; i++) { - BC_VECTOR divided = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; + BC_VECTOR tmp_numerator = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ - if (div_carry == 0 && divided < divider) { + if (div_carry == 0 && tmp_numerator < tmp_divisor) { quot_vector[quot_top_index - i] = 0; div_carry = numerator_vector[numerator_top_index - i]; numerator_vector[numerator_top_index - i] = 0; @@ -186,13 +186,13 @@ static inline void bc_standard_div( /* * Determine the temporary quotient. - * "divided" is numerator_high in the previous example. The maximum value of numerator_high is divisor_high * B, - * so here it is divider * BC_VECTOR_BOUNDARY_NUM. - * The number of digits in divider is BC_VECTOR_SIZE + 1, so even if divider is at its maximum value, + * "tmp_numerator" is numerator_high in the previous example. The maximum value of numerator_high is divisor_high * B, + * so here it is tmp_divisor * BC_VECTOR_BOUNDARY_NUM. + * The number of digits in tmp_divisor is BC_VECTOR_SIZE + 1, so even if tmp_divisor is at its maximum value, * it will never overflow here. */ - divided += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; - BC_VECTOR quot_guess = divided / divider; + tmp_numerator += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; + BC_VECTOR quot_guess = tmp_numerator / tmp_divisor; /* For sub, add the remainder to the high-order digit */ numerator_vector[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; From 0102defc24013d7e9ac35e23093072a74a2e32a8 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 09:00:56 +0900 Subject: [PATCH 07/34] Delete unnecessary i and j definitions --- ext/bcmath/libbcmath/src/div.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index f9b4d1f1924d9..a4a5daea31c89 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -74,7 +74,6 @@ static inline void bc_standard_div( BC_VECTOR *divisor_vector, size_t divisor_arr_size, size_t divisor_len, BC_VECTOR *quot_vector, size_t quot_arr_size ) { - size_t i, j; 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; @@ -173,7 +172,7 @@ static inline void bc_standard_div( size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = POW_10_LUT[divisor_top_digits - 1]; BC_VECTOR tmp_divisor = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; - for (i = 0; i < quot_arr_size; i++) { + for (size_t i = 0; i < quot_arr_size; i++) { BC_VECTOR tmp_numerator = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ @@ -212,6 +211,7 @@ static inline void bc_standard_div( BC_VECTOR sub; BC_VECTOR borrow = 0; BC_VECTOR *numerator_calc_bottom = numerator_vector + numerator_arr_size - divisor_arr_size - i; + size_t j; for (j = 0; j < divisor_arr_size - 1; j++) { sub = divisor_vector[j] * quot_guess + borrow; BC_VECTOR sub_low = sub % BC_VECTOR_BOUNDARY_NUM; @@ -250,8 +250,6 @@ static inline void bc_standard_div( static void bc_do_div(char *n1, size_t n1_readable_len, size_t n1_bottom_extension, char *n2, size_t n2_len, bc_num *quot, size_t quot_len) { - size_t i; - size_t n2_arr_size = (n2_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; size_t n1_arr_size = (n1_readable_len + n1_bottom_extension + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; size_t quot_arr_size = n1_arr_size - n2_arr_size + 1; @@ -277,7 +275,7 @@ static void bc_do_div(char *n1, size_t n1_readable_len, size_t n1_bottom_extensi n1_read = MIN(n1_bottom_read_len, n1_readable_len); base = POW_10_LUT[n1_bottom_extension]; n1_vector[n1_vector_count] = 0; - for (i = 0; i < n1_read; i++) { + for (size_t i = 0; i < n1_read; i++) { n1_vector[n1_vector_count] += *n1 * base; base *= BASE; n1--; @@ -302,6 +300,7 @@ static void bc_do_div(char *n1, size_t n1_readable_len, size_t n1_bottom_extensi char *qptr = (*quot)->n_value; char *qend = qptr + quot_len - 1; + size_t i; for (i = 0; i < quot_real_arr_size - 1; i++) { #if BC_VECTOR_SIZE == 4 bc_write_bcd_representation(quot_vector[i], qend - 3); From 523b3dc57907267571fb89d7c92a849625bca840 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 09:03:17 +0900 Subject: [PATCH 08/34] Fixed comments --- ext/bcmath/libbcmath/src/div.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index a4a5daea31c89..85dd0469b5faf 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -125,11 +125,11 @@ static inline void bc_standard_div( * The error E can be expressed by the following formula. * * E = (numerator_high * B^k) / (divisor_high * B^k) - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) - * E = numerator_high / divisor_high - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) - * E = (numerator_high * (divisor_high * B^k + divisor_low) - (numerator_high * B^k + numerator_low) * divisor_high) / (divisor_high * (divisor_high * B^k + divisor_low)) - * E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * = numerator_high / divisor_high - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) + * = (numerator_high * (divisor_high * B^k + divisor_low) - (numerator_high * B^k + numerator_low) * divisor_high) / (divisor_high * (divisor_high * B^k + divisor_low)) + * = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * - * Find the error MAX_E when the error E is maximum. + * Find the error MAX_E when the error E is maximum (0 <= E <= MAX_E). * First, numerator_high, which only exists in the numerator, uses its maximum value. * Considering carry-back, numerator_high can be expressed as follows. * numerator_high = divisor_high * B @@ -137,19 +137,19 @@ static inline void bc_standard_div( * numerator_low = 0 * * MAX_E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) - * MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * * divisor_low uses the maximum value. * divisor_low = B^k - 1 * MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) - * MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) + * = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) * * divisor_high uses the minimum value, but want to see how the number of digits affects the error, so represent * the minimum value as: * divisor_high = 10^x (any x ∈ ℕ) * Since B = 10^b, the formula becomes: * MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) - * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) + * = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) * * Now let's make an approximation. Remove -1 from the numerator. Approximate the numerator to be * large and the denominator to be small, such that MAX_E is less than this expression. @@ -157,8 +157,8 @@ static inline void bc_standard_div( * * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) * MAX_E < (10^x * 10^b * (10^b)^k) / (10^x * 10^x * (10^b)^k) - * MAX_E < 10^b / 10^x - * MAX_E < 10^(b - x) + * < 10^b / 10^x + * < 10^(b - x) * * Therefore, found that if the number of digits in base B and divisor_high are equal, the error will be less * than 1 regardless of the number of digits in the value of k. From 62a0ecf68b844c10589593f9e4c6dd97ad613d02 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 11:47:08 +0900 Subject: [PATCH 09/34] use <=> --- ext/bcmath/libbcmath/src/div.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 85dd0469b5faf..02b9f9007f0af 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -125,9 +125,9 @@ static inline void bc_standard_div( * The error E can be expressed by the following formula. * * E = (numerator_high * B^k) / (divisor_high * B^k) - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) - * = numerator_high / divisor_high - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) - * = (numerator_high * (divisor_high * B^k + divisor_low) - (numerator_high * B^k + numerator_low) * divisor_high) / (divisor_high * (divisor_high * B^k + divisor_low)) - * = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * <=> E = numerator_high / divisor_high - (numerator_high * B^k + numerator_low) / (divisor_high * B^k + divisor_low) + * <=> E = (numerator_high * (divisor_high * B^k + divisor_low) - (numerator_high * B^k + numerator_low) * divisor_high) / (divisor_high * (divisor_high * B^k + divisor_low)) + * <=> E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * * Find the error MAX_E when the error E is maximum (0 <= E <= MAX_E). * First, numerator_high, which only exists in the numerator, uses its maximum value. @@ -137,19 +137,19 @@ static inline void bc_standard_div( * numerator_low = 0 * * MAX_E = (numerator_high * divisor_low - divisor_high * numerator_low) / (divisor_high * (divisor_high * B^k + divisor_low)) - * = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) + * <=> MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) * * divisor_low uses the maximum value. * divisor_low = B^k - 1 * MAX_E = (divisor_high * B * divisor_low) / (divisor_high * (divisor_high * B^k + divisor_low)) - * = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) + * <=> MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) * * divisor_high uses the minimum value, but want to see how the number of digits affects the error, so represent * the minimum value as: * divisor_high = 10^x (any x ∈ ℕ) * Since B = 10^b, the formula becomes: * MAX_E = (divisor_high * B * (B^k - 1)) / (divisor_high * (divisor_high * B^k + B^k - 1)) - * = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) + * <=> MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) * * Now let's make an approximation. Remove -1 from the numerator. Approximate the numerator to be * large and the denominator to be small, such that MAX_E is less than this expression. @@ -157,8 +157,8 @@ static inline void bc_standard_div( * * MAX_E = (10^x * 10^b * ((10^b)^k - 1)) / (10^x * (10^x * (10^b)^k + (10^b)^k - 1)) * MAX_E < (10^x * 10^b * (10^b)^k) / (10^x * 10^x * (10^b)^k) - * < 10^b / 10^x - * < 10^(b - x) + * <=> MAX_E < 10^b / 10^x + * <=> MAX_E < 10^(b - x) * * Therefore, found that if the number of digits in base B and divisor_high are equal, the error will be less * than 1 regardless of the number of digits in the value of k. From 93e34b77211096a9d44539ff8aa658f58284e264 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 11:49:09 +0900 Subject: [PATCH 10/34] Char made const --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 02b9f9007f0af..3eea3c4b5b637 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -248,7 +248,7 @@ static inline void bc_standard_div( } } -static void bc_do_div(char *n1, size_t n1_readable_len, size_t n1_bottom_extension, char *n2, size_t n2_len, bc_num *quot, size_t quot_len) +static void bc_do_div(const char *n1, size_t n1_readable_len, size_t n1_bottom_extension, const char *n2, size_t n2_len, bc_num *quot, size_t quot_len) { size_t n2_arr_size = (n2_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; size_t n1_arr_size = (n1_readable_len + n1_bottom_extension + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; From 8972f5809869917f7eac8f288bd3ae6b15f3be91 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:04:51 +0900 Subject: [PATCH 11/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 3eea3c4b5b637..2853f1b72d1e9 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -42,7 +42,7 @@ static const BC_VECTOR POW_10_LUT[9] = { }; /* - * Use this function when the number of array elements in divisor is 1, that is, when divisor is not divided into multiple array elements. + * This function should be used when the divisor is not split into multiple chunks, i.e. when the size of the array is one. * The algorithm can be simplified if divisor is not divided, and this function is an optimization of bc_standard_div. */ static inline void bc_fast_div( From 960fc8ef99320876c1ebc3f50f227b87c95c4805 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:05:13 +0900 Subject: [PATCH 12/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 2853f1b72d1e9..b4fed564c549a 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -43,7 +43,8 @@ static const BC_VECTOR POW_10_LUT[9] = { /* * This function should be used when the divisor is not split into multiple chunks, i.e. when the size of the array is one. - * The algorithm can be simplified if divisor is not divided, and this function is an optimization of bc_standard_div. + * 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_vector, size_t numerator_arr_size, BC_VECTOR divisor_vector, BC_VECTOR *quot_vector, size_t quot_arr_size) From 5aeec239efd81ed249974cf913bdc5ec8279e259 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:08:42 +0900 Subject: [PATCH 13/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index b4fed564c549a..c813b9ef6c66e 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -68,7 +68,7 @@ static inline void bc_fast_div( } /* - * Used when the number of elements in the divisor vector is 2 or more. + * Used when the divisor is split into 2 or more chunks. */ static inline void bc_standard_div( BC_VECTOR *numerator_vector, size_t numerator_arr_size, From 919668527975986c9fc939ffe2d9a249c4475f82 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:12:10 +0900 Subject: [PATCH 14/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index c813b9ef6c66e..d3d15a3fed3ad 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -82,7 +82,7 @@ static inline void bc_standard_div( BC_VECTOR div_carry = 0; /* - * Consider the error E between the actual quotient and the temporary quotient. + * Errors might occur between the true quotient and the temporary quotient calculated using only the high order digits. * For example, the quotient of 240 / 121 is 1, but if calculate the quotient only using the high-order digits (24 / 12), * it will be 2. In this case, 2 is the temporary quotient. Here, the error E is 1. * Also note that for example 2400000 / 120, there will be 5 divisions. From 5b70842c05c317f98df601754200c5d70470d7b1 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 12:26:37 +0900 Subject: [PATCH 15/34] Fixed var names --- ext/bcmath/libbcmath/src/div.c | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index d3d15a3fed3ad..3d2b2263130e2 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -172,12 +172,12 @@ static inline void bc_standard_div( } size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = POW_10_LUT[divisor_top_digits - 1]; - BC_VECTOR tmp_divisor = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; + BC_VECTOR divisor_high_part = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; for (size_t i = 0; i < quot_arr_size; i++) { - BC_VECTOR tmp_numerator = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; + BC_VECTOR numerator_high_part = numerator_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ - if (div_carry == 0 && tmp_numerator < tmp_divisor) { + if (div_carry == 0 && numerator_high_part < divisor_high_part) { quot_vector[quot_top_index - i] = 0; div_carry = numerator_vector[numerator_top_index - i]; numerator_vector[numerator_top_index - i] = 0; @@ -186,13 +186,13 @@ static inline void bc_standard_div( /* * Determine the temporary quotient. - * "tmp_numerator" is numerator_high in the previous example. The maximum value of numerator_high is divisor_high * B, - * so here it is tmp_divisor * BC_VECTOR_BOUNDARY_NUM. - * The number of digits in tmp_divisor is BC_VECTOR_SIZE + 1, so even if tmp_divisor is at its maximum value, + * The maximum value of numerator_high is divisor_high * B in the previous example, so here numerator_high_part is + * divisor_high_part * BC_VECTOR_BOUNDARY_NUM. + * The number of digits in divisor_high_part is BC_VECTOR_SIZE + 1, so even if divisor_high_part is at its maximum value, * it will never overflow here. */ - tmp_numerator += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; - BC_VECTOR quot_guess = tmp_numerator / tmp_divisor; + numerator_high_part += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; + BC_VECTOR quot_guess = numerator_high_part / divisor_high_part; /* For sub, add the remainder to the high-order digit */ numerator_vector[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; From 28bf55df6c70b99d179e18bddeaaffd41d077526 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Thu, 18 Jul 2024 12:28:52 +0900 Subject: [PATCH 16/34] Fixed confusing comments --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 3d2b2263130e2..4b398699e7cff 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -54,7 +54,7 @@ static inline void bc_fast_div( for (size_t i = 0; i < quot_arr_size - 1; i++) { if (numerator_vector[numerator_top_index - i] < divisor_vector) { quot_vector[quot_top_index - i] = 0; - /* numerator_vector[numerator_top_index - i] is always a smaller number than divisor_vector. Therefore, there will be no overflow. */ + /* numerator_vector[numerator_top_index - i] < divisor_vector, so there will be no overflow. */ numerator_vector[numerator_top_index - i - 1] += numerator_vector[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM; continue; } From eeb72cac8cab3f87d41bd64f398ffd8183891087 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Thu, 18 Jul 2024 12:32:11 +0900 Subject: [PATCH 17/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 4b398699e7cff..f6cbf011bc3fa 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -83,8 +83,11 @@ static inline void bc_standard_div( /* * Errors might occur between the true quotient and the temporary quotient calculated using only the high order digits. - * For example, the quotient of 240 / 121 is 1, but if calculate the quotient only using the high-order digits (24 / 12), - * it will be 2. In this case, 2 is the temporary quotient. Here, the error E is 1. + * For example, the true quotient of 240 / 121 is 1. + * However, if it is calculated using only the high-order digits (24 / 12), we would get 2. + * We refer to this value of 2 as the temporary quotient. + * We also define E to be error between the true quotient and the temporary quotient, + * which in this case, is 1. * Also note that for example 2400000 / 120, there will be 5 divisions. * * Another example: Calculating 99900000000 / 10009999 with base 10000. From 401ad242b326b6293d05d0c6a8bfbf0d457129a7 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Fri, 19 Jul 2024 08:56:43 +0900 Subject: [PATCH 18/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index f6cbf011bc3fa..ecaad4837a5ac 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -90,10 +90,10 @@ static inline void bc_standard_div( * which in this case, is 1. * Also note that for example 2400000 / 120, there will be 5 divisions. * - * Another example: Calculating 99900000000 / 10009999 with base 10000. - * The actual quotient is 9980, but if calculate it using only the high-order digits: - * 999 / 1000 = 0, Numbers that are not divisible are carried forward to the next division. - * 9990000 / 1000 = 9990 (Here, the temporary quotient is 9990 and the error E is 10.) + * Another example: Calculating 999_0000_0000 / 1000_9999 with base 10000. + * The true quotient is 9980, but if it is calculated using only the high-order digits (999 / 1000), we would get 0 + * If the temporary quotient is 0, we need to carry the digits to the next division, which is 999_0000 / 1000. + * The new temporary quotient we get is 9990, with error E = 10. * * If calculate the temporary quotient using only one array of numerator and divisor, the error E can be larger than 1. * In other words, in the restoring division, the count of additions for restore increases significantly. From e0b83adc223ae609061e031e27d80d554b385466 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Fri, 19 Jul 2024 08:52:52 +0900 Subject: [PATCH 19/34] Fixed vector to vectors --- ext/bcmath/libbcmath/src/div.c | 60 +++++++++++++++++----------------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index ecaad4837a5ac..f140d2ad5ea40 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -47,33 +47,33 @@ static const BC_VECTOR POW_10_LUT[9] = { * This function is therefore an optimized version of bc_standard_div(). */ static inline void bc_fast_div( - BC_VECTOR *numerator_vector, size_t numerator_arr_size, BC_VECTOR divisor_vector, BC_VECTOR *quot_vector, size_t quot_arr_size) + BC_VECTOR *numerator_vectors, size_t numerator_arr_size, BC_VECTOR divisor_vector, BC_VECTOR *quot_vectors, size_t quot_arr_size) { size_t numerator_top_index = numerator_arr_size - 1; size_t quot_top_index = quot_arr_size - 1; for (size_t i = 0; i < quot_arr_size - 1; i++) { - if (numerator_vector[numerator_top_index - i] < divisor_vector) { - quot_vector[quot_top_index - i] = 0; - /* numerator_vector[numerator_top_index - i] < divisor_vector, so there will be no overflow. */ - numerator_vector[numerator_top_index - i - 1] += numerator_vector[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM; + 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; continue; } - quot_vector[quot_top_index - i] = numerator_vector[numerator_top_index - i] / divisor_vector; - numerator_vector[numerator_top_index - i] -= divisor_vector * quot_vector[quot_top_index - i]; - numerator_vector[numerator_top_index - i - 1] += numerator_vector[numerator_top_index - i] * BC_VECTOR_BOUNDARY_NUM; - numerator_vector[numerator_top_index - i] = 0; + 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; } /* last */ - quot_vector[0] = numerator_vector[0] / divisor_vector; + quot_vectors[0] = numerator_vectors[0] / divisor_vector; } /* * Used when the divisor is split into 2 or more chunks. */ static inline void bc_standard_div( - BC_VECTOR *numerator_vector, size_t numerator_arr_size, - BC_VECTOR *divisor_vector, size_t divisor_arr_size, size_t divisor_len, - BC_VECTOR *quot_vector, size_t quot_arr_size + BC_VECTOR *numerator_vectors, size_t numerator_arr_size, + 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; @@ -175,15 +175,15 @@ static inline void bc_standard_div( } size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = POW_10_LUT[divisor_top_digits - 1]; - BC_VECTOR divisor_high_part = divisor_vector[divisor_top_index] * high_part_shift + divisor_vector[divisor_top_index - 1] / low_part_shift; + 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_vector[numerator_top_index - i] * high_part_shift + numerator_vector[numerator_top_index - i - 1] / low_part_shift; + BC_VECTOR numerator_high_part = numerator_vectors[numerator_top_index - i] * high_part_shift + numerator_vectors[numerator_top_index - i - 1] / low_part_shift; /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ if (div_carry == 0 && numerator_high_part < divisor_high_part) { - quot_vector[quot_top_index - i] = 0; - div_carry = numerator_vector[numerator_top_index - i]; - numerator_vector[numerator_top_index - i] = 0; + quot_vectors[quot_top_index - i] = 0; + div_carry = numerator_vectors[numerator_top_index - i]; + numerator_vectors[numerator_top_index - i] = 0; continue; } @@ -198,26 +198,26 @@ static inline void bc_standard_div( BC_VECTOR quot_guess = numerator_high_part / divisor_high_part; /* For sub, add the remainder to the high-order digit */ - numerator_vector[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; + numerator_vectors[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; /* * The temporary quotient can only have errors in the "too large" direction. * So if the temporary quotient is 0 here, the quotient is 0. */ if (quot_guess == 0) { - quot_vector[quot_top_index - i] = 0; - div_carry = numerator_vector[numerator_top_index - i]; - numerator_vector[numerator_top_index - i] = 0; + quot_vectors[quot_top_index - i] = 0; + div_carry = numerator_vectors[numerator_top_index - i]; + numerator_vectors[numerator_top_index - i] = 0; continue; } /* sub */ BC_VECTOR sub; BC_VECTOR borrow = 0; - BC_VECTOR *numerator_calc_bottom = numerator_vector + numerator_arr_size - divisor_arr_size - i; + BC_VECTOR *numerator_calc_bottom = numerator_vectors + numerator_arr_size - divisor_arr_size - i; size_t j; for (j = 0; j < divisor_arr_size - 1; j++) { - sub = divisor_vector[j] * quot_guess + borrow; + sub = divisor_vectors[j] * quot_guess + borrow; BC_VECTOR sub_low = sub % BC_VECTOR_BOUNDARY_NUM; borrow = sub / BC_VECTOR_BOUNDARY_NUM; @@ -229,7 +229,7 @@ static inline void bc_standard_div( } } /* last digit sub */ - sub = divisor_vector[j] * quot_guess + borrow; + sub = divisor_vectors[j] * quot_guess + borrow; bool neg_remainder = numerator_calc_bottom[j] < sub; numerator_calc_bottom[j] -= sub; @@ -238,17 +238,17 @@ static inline void bc_standard_div( if (neg_remainder) { quot_guess--; for (j = 0; j < divisor_arr_size - 1; j++) { - numerator_calc_bottom[j] += divisor_vector[j] + carry; + numerator_calc_bottom[j] += divisor_vectors[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_vector[j] + carry; + numerator_calc_bottom[j] += divisor_vectors[j] + carry; } - quot_vector[quot_top_index - i] = quot_guess; - div_carry = numerator_vector[numerator_top_index - i]; - numerator_vector[numerator_top_index - i] = 0; + quot_vectors[quot_top_index - i] = quot_guess; + div_carry = numerator_vectors[numerator_top_index - i]; + numerator_vectors[numerator_top_index - i] = 0; } } From 06b6f5624797f9ac64629081f5b49ba23ccbaf81 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Fri, 19 Jul 2024 08:53:18 +0900 Subject: [PATCH 20/34] Removed unnecessary comments --- ext/bcmath/libbcmath/src/div.c | 1 - 1 file changed, 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index f140d2ad5ea40..d4e6d81fd60e6 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -88,7 +88,6 @@ static inline void bc_standard_div( * We refer to this value of 2 as the temporary quotient. * We also define E to be error between the true quotient and the temporary quotient, * which in this case, is 1. - * Also note that for example 2400000 / 120, there will be 5 divisions. * * Another example: Calculating 999_0000_0000 / 1000_9999 with base 10000. * The true quotient is 9980, but if it is calculated using only the high-order digits (999 / 1000), we would get 0 From dfb1e4fb0b36a5d05160245b7ee37f73e7639a85 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Fri, 19 Jul 2024 08:54:58 +0900 Subject: [PATCH 21/34] Calculate numerator_high_part first --- ext/bcmath/libbcmath/src/div.c | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index d4e6d81fd60e6..8d70e6fe66096 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -178,14 +178,6 @@ static inline void bc_standard_div( 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; - /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ - if (div_carry == 0 && 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; - continue; - } - /* * Determine the temporary quotient. * The maximum value of numerator_high is divisor_high * B in the previous example, so here numerator_high_part is @@ -194,6 +186,15 @@ static inline void bc_standard_div( * it will never overflow here. */ numerator_high_part += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; + + /* If it is clear that divisor is greater in this loop, then the quotient is 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; + continue; + } + BC_VECTOR quot_guess = numerator_high_part / divisor_high_part; /* For sub, add the remainder to the high-order digit */ From 859a9dbf5cf63e282a7b22fda02dc96039d086c7 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Fri, 19 Jul 2024 08:55:35 +0900 Subject: [PATCH 22/34] Fixed a comment --- ext/bcmath/libbcmath/src/div.c | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 8d70e6fe66096..cd3a0a5ee58f5 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -67,9 +67,7 @@ static inline void bc_fast_div( quot_vectors[0] = numerator_vectors[0] / divisor_vector; } -/* - * Used when the divisor is split into 2 or more chunks. - */ +/* Used when the divisor is split into 2 or more chunks. */ static inline void bc_standard_div( BC_VECTOR *numerator_vectors, size_t numerator_arr_size, BC_VECTOR *divisor_vectors, size_t divisor_arr_size, size_t divisor_len, From f24a72bb68cf8af7d2add328643749a58b9f201c Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Fri, 19 Jul 2024 09:11:24 +0900 Subject: [PATCH 23/34] Fixed a comment --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index cd3a0a5ee58f5..5a72e99a06612 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -185,7 +185,7 @@ static inline void bc_standard_div( */ numerator_high_part += div_carry * BC_VECTOR_BOUNDARY_NUM * high_part_shift; - /* If it is clear that divisor is greater in this loop, then the quotient is 0. */ + /* 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]; From 3cca73df462db70809dbe611ba977fddf5c9fabe Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Sat, 20 Jul 2024 18:15:26 +0900 Subject: [PATCH 24/34] Added a comment that using restoring division --- ext/bcmath/libbcmath/src/div.c | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 5a72e99a06612..35ac6554847af 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -67,7 +67,11 @@ static inline void bc_fast_div( quot_vectors[0] = numerator_vectors[0] / divisor_vector; } -/* Used when the divisor is split into 2 or more chunks. */ +/* + * Used when the divisor is split into 2 or more chunks. + * This use the restoring division algorithm. + * https://en.wikipedia.org/wiki/Division_algorithm#Restoring_division + */ static inline void bc_standard_div( BC_VECTOR *numerator_vectors, size_t numerator_arr_size, BC_VECTOR *divisor_vectors, size_t divisor_arr_size, size_t divisor_len, From bd9b1bb59ba5f7f66b14b893248157f33e42c2b9 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Sat, 20 Jul 2024 18:17:20 +0900 Subject: [PATCH 25/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 35ac6554847af..65e08511f5424 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -203,8 +203,8 @@ static inline void bc_standard_div( numerator_vectors[numerator_top_index - i] += div_carry * BC_VECTOR_BOUNDARY_NUM; /* - * The temporary quotient can only have errors in the "too large" direction. - * So if the temporary quotient is 0 here, the quotient is 0. + * It is impossible for the temporary quotient to underestimate the true quotient. + * Therefore a temporary quotient of 0 implies the true quotient it also 0. */ if (quot_guess == 0) { quot_vectors[quot_top_index - i] = 0; From f07abf459c61523cd1d3d57933d7bbbb331fbfe5 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Sat, 20 Jul 2024 18:18:14 +0900 Subject: [PATCH 26/34] Added a blank line Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 1 + 1 file changed, 1 insertion(+) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 65e08511f5424..daf12df8f47d0 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -174,6 +174,7 @@ static inline void bc_standard_div( if (divisor_top_digits == 0) { divisor_top_digits = BC_VECTOR_SIZE; } + size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = 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; From 8ca6f1ed55b442330924cf1e5a03a09284ba3fd6 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Sat, 20 Jul 2024 18:19:27 +0900 Subject: [PATCH 27/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 1 + 1 file changed, 1 insertion(+) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index daf12df8f47d0..c3a3207a93bb0 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -170,6 +170,7 @@ static inline void bc_standard_div( * * Here, B is BC_VECTOR_BOUNDARY_NUM, so adjust the number of digits in high part of divisor to be BC_VECTOR_SIZE + 1. */ + /* the number of usable digits, thus divisor_len % BC_VECTOR_SIZE == 0 implies we have BC_VECTOR_SIZE usable digits */ size_t divisor_top_digits = divisor_len % BC_VECTOR_SIZE; if (divisor_top_digits == 0) { divisor_top_digits = BC_VECTOR_SIZE; From ffb6e8ea3d7142e4082c69b4f241f273f7347048 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Sat, 20 Jul 2024 18:38:09 +0900 Subject: [PATCH 28/34] Fixed var names --- ext/bcmath/libbcmath/src/div.c | 76 +++++++++++++++++----------------- 1 file changed, 39 insertions(+), 37 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index c3a3207a93bb0..80835d660dcc3 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -175,7 +175,7 @@ static inline void bc_standard_div( if (divisor_top_digits == 0) { divisor_top_digits = BC_VECTOR_SIZE; } - + size_t high_part_shift = POW_10_LUT[BC_VECTOR_SIZE - divisor_top_digits + 1]; size_t low_part_shift = 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; @@ -256,52 +256,54 @@ static inline void bc_standard_div( } } -static void bc_do_div(const char *n1, size_t n1_readable_len, size_t n1_bottom_extension, const char *n2, size_t n2_len, bc_num *quot, size_t quot_len) -{ - size_t n2_arr_size = (n2_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; - size_t n1_arr_size = (n1_readable_len + n1_bottom_extension + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; - size_t quot_arr_size = n1_arr_size - n2_arr_size + 1; +static void bc_do_div( + const char *numerator, size_t numerator_readable_len, size_t numerator_bottom_extension, + const char *divisor, size_t divisor_len, bc_num *quot, size_t quot_len +) { + size_t divisor_arr_size = (divisor_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; + size_t numerator_arr_size = (numerator_readable_len + numerator_bottom_extension + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE; + size_t quot_arr_size = numerator_arr_size - divisor_arr_size + 1; size_t quot_real_arr_size = MIN(quot_arr_size, (quot_len + BC_VECTOR_SIZE - 1) / BC_VECTOR_SIZE); - BC_VECTOR *n1_vector = safe_emalloc(n1_arr_size + n2_arr_size + quot_arr_size, sizeof(BC_VECTOR), 0); - BC_VECTOR *n2_vector = n1_vector + n1_arr_size; - BC_VECTOR *quot_vector = n2_vector + n2_arr_size; + BC_VECTOR *numerator_vectors = safe_emalloc(numerator_arr_size + divisor_arr_size + quot_arr_size, sizeof(BC_VECTOR), 0); + BC_VECTOR *divisor_vectors = numerator_vectors + numerator_arr_size; + BC_VECTOR *quot_vectors = divisor_vectors + divisor_arr_size; /* Fill with zeros and convert as many vector elements as needed */ - size_t n1_vector_count = 0; - while (n1_bottom_extension >= BC_VECTOR_SIZE) { - n1_vector[n1_vector_count] = 0; - n1_bottom_extension -= BC_VECTOR_SIZE; - n1_vector_count++; + size_t numerator_vector_count = 0; + while (numerator_bottom_extension >= BC_VECTOR_SIZE) { + numerator_vectors[numerator_vector_count] = 0; + numerator_bottom_extension -= BC_VECTOR_SIZE; + numerator_vector_count++; } - size_t n1_bottom_read_len = BC_VECTOR_SIZE - n1_bottom_extension; + size_t numerator_bottom_read_len = BC_VECTOR_SIZE - numerator_bottom_extension; size_t base; - size_t n1_read = 0; - if (n1_bottom_read_len < BC_VECTOR_SIZE) { - n1_read = MIN(n1_bottom_read_len, n1_readable_len); - base = POW_10_LUT[n1_bottom_extension]; - n1_vector[n1_vector_count] = 0; - for (size_t i = 0; i < n1_read; i++) { - n1_vector[n1_vector_count] += *n1 * base; + size_t numerator_read = 0; + if (numerator_bottom_read_len < BC_VECTOR_SIZE) { + numerator_read = MIN(numerator_bottom_read_len, numerator_readable_len); + base = POW_10_LUT[numerator_bottom_extension]; + numerator_vectors[numerator_vector_count] = 0; + for (size_t i = 0; i < numerator_read; i++) { + numerator_vectors[numerator_vector_count] += *numerator * base; base *= BASE; - n1--; + numerator--; } - n1_vector_count++; + numerator_vector_count++; } - /* Bulk convert n1 and n2 to vectors */ - if (n1_readable_len > n1_read) { - bc_convert_to_vector(n1_vector + n1_vector_count, n1, n1_readable_len - n1_read); + /* Bulk convert numerator and divisor to vectors */ + if (numerator_readable_len > numerator_read) { + bc_convert_to_vector(numerator_vectors + numerator_vector_count, numerator, numerator_readable_len - numerator_read); } - bc_convert_to_vector(n2_vector, n2, n2_len); + bc_convert_to_vector(divisor_vectors, divisor, divisor_len); /* Do the divide */ - if (n2_arr_size == 1) { - bc_fast_div(n1_vector, n1_arr_size, n2_vector[0], quot_vector, quot_arr_size); + if (divisor_arr_size == 1) { + bc_fast_div(numerator_vectors, numerator_arr_size, divisor_vectors[0], quot_vectors, quot_arr_size); } else { - bc_standard_div(n1_vector, n1_arr_size, n2_vector, n2_arr_size, n2_len, quot_vector, quot_arr_size); + bc_standard_div(numerator_vectors, numerator_arr_size, divisor_vectors, divisor_arr_size, divisor_len, quot_vectors, quot_arr_size); } /* Convert to bc_num */ @@ -311,21 +313,21 @@ static void bc_do_div(const char *n1, size_t n1_readable_len, size_t n1_bottom_e size_t i; for (i = 0; i < quot_real_arr_size - 1; i++) { #if BC_VECTOR_SIZE == 4 - bc_write_bcd_representation(quot_vector[i], qend - 3); + bc_write_bcd_representation(quot_vectors[i], qend - 3); qend -= 4; #else - bc_write_bcd_representation(quot_vector[i] / 10000, qend - 7); - bc_write_bcd_representation(quot_vector[i] % 10000, qend - 3); + bc_write_bcd_representation(quot_vectors[i] / 10000, qend - 7); + bc_write_bcd_representation(quot_vectors[i] % 10000, qend - 3); qend -= 8; #endif } while (qend >= qptr) { - *qend-- = quot_vector[i] % BASE; - quot_vector[i] /= BASE; + *qend-- = quot_vectors[i] % BASE; + quot_vectors[i] /= BASE; } - efree(n1_vector); + efree(numerator_vectors); } bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, size_t scale) From 72f7f3a1c7f810fa6495ecf09f6e6f6b3d94d6c9 Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Mon, 22 Jul 2024 11:06:19 +0900 Subject: [PATCH 29/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 80835d660dcc3..0ca906cb794c3 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -96,8 +96,7 @@ static inline void bc_standard_div( * If the temporary quotient is 0, we need to carry the digits to the next division, which is 999_0000 / 1000. * The new temporary quotient we get is 9990, with error E = 10. * - * If calculate the temporary quotient using only one array of numerator and divisor, the error E can be larger than 1. - * In other words, in the restoring division, the count of additions for restore increases significantly. + * Because we use the restoring division we need to perform E restorations, which can be significant if E is large. * * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to * at most one, adjust the number of high-order digits used to calculate the temporary quotient as follows. From b88956d80b5afbbe024658eecaab9ef0a2f9e9ab Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Mon, 22 Jul 2024 11:07:47 +0900 Subject: [PATCH 30/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 0ca906cb794c3..a0382c9a8e7e9 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -205,7 +205,7 @@ static inline void bc_standard_div( /* * It is impossible for the temporary quotient to underestimate the true quotient. - * Therefore a temporary quotient of 0 implies the true quotient it also 0. + * Therefore a temporary quotient of 0 implies the true quotient is also 0. */ if (quot_guess == 0) { quot_vectors[quot_top_index - i] = 0; From 3f3679aa7fd9bb52b8e8310c9eb8f9599569a27f Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Mon, 22 Jul 2024 11:08:21 +0900 Subject: [PATCH 31/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index a0382c9a8e7e9..8d64f1311adb3 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -298,7 +298,7 @@ static void bc_do_div( } bc_convert_to_vector(divisor_vectors, divisor, divisor_len); - /* Do the divide */ + /* Do the division */ if (divisor_arr_size == 1) { bc_fast_div(numerator_vectors, numerator_arr_size, divisor_vectors[0], quot_vectors, quot_arr_size); } else { From 71a3ba3de58592ee9567ae5ee68994561208ed2b Mon Sep 17 00:00:00 2001 From: Saki Takamachi <34942839+SakiTakamachi@users.noreply.github.com> Date: Mon, 22 Jul 2024 11:09:29 +0900 Subject: [PATCH 32/34] Fixed a comment Co-authored-by: Gina Peter Banyard --- ext/bcmath/libbcmath/src/div.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 8d64f1311adb3..084763feba816 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -98,11 +98,10 @@ static inline void bc_standard_div( * * Because we use the restoring division we need to perform E restorations, which can be significant if E is large. * - * Therefore, in order to keep the error within 1 and to limit the number of additions required for restoration to - * at most one, adjust the number of high-order digits used to calculate the temporary quotient as follows. - * - Adjust the number of digits of divisor used in the calculation to BC_VECTOR_SIZE + 1 digit. The missing digits are - * filled in from the next array element. - * - Add digits to numerator in the same way as the number of digits adjusted by divisor. + * Therefore, for the error E to be at most 1 we adjust the number of high-order digits used to calculate the temporary quotient as follows: + * - Include BC_VECTOR_SIZE + 1 digits of the divisor used in the calculation of the temporary quotient. + The missing digits are filled in from the next array element. + * - Adjust the number of digits in the numerator similarly to what was done for the divisor. * * e.g. * numerator = 123456780000 From 68bcff9eb855af39c3efa19f6cc0e41c8263f5d4 Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Mon, 22 Jul 2024 11:19:17 +0900 Subject: [PATCH 33/34] Renamed var names --- ext/bcmath/libbcmath/src/div.c | 170 ++++++++++++++++----------------- 1 file changed, 85 insertions(+), 85 deletions(-) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index 084763feba816..f2465b639cac6 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -328,146 +328,146 @@ static void bc_do_div( efree(numerator_vectors); } -bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, size_t scale) +bool bc_divide(bc_num numerator, bc_num divisor, bc_num *quot, size_t scale) { /* divide by zero */ - if (bc_is_zero(n2)) { + if (bc_is_zero(divisor)) { return false; } bc_free_num(quot); - /* If n1 is zero, the quotient is always zero. */ - if (bc_is_zero(n1)) { + /* If numerator is zero, the quotient is always zero. */ + if (bc_is_zero(numerator)) { *quot = bc_copy_num(BCG(_zero_)); return true; } - char *n1ptr = n1->n_value; - char *n1end = n1ptr + n1->n_len + n1->n_scale - 1; - size_t n1_len = n1->n_len; - size_t n1_scale = n1->n_scale; - - char *n2ptr = n2->n_value; - char *n2end = n2ptr + n2->n_len + n2->n_scale - 1; - size_t n2_len = n2->n_len; - size_t n2_scale = n2->n_scale; - size_t n2_int_right_zeros = 0; - - /* remove n2 trailing zeros */ - while (*n2end == 0 && n2_scale > 0) { - n2end--; - n2_scale--; + char *numeratorptr = numerator->n_value; + char *numeratorend = numeratorptr + numerator->n_len + numerator->n_scale - 1; + size_t numerator_len = numerator->n_len; + size_t numerator_scale = numerator->n_scale; + + char *divisorptr = divisor->n_value; + char *divisorend = divisorptr + divisor->n_len + divisor->n_scale - 1; + size_t divisor_len = divisor->n_len; + size_t divisor_scale = divisor->n_scale; + size_t divisor_int_right_zeros = 0; + + /* remove divisor trailing zeros */ + while (*divisorend == 0 && divisor_scale > 0) { + divisorend--; + divisor_scale--; } - while (*n2end == 0) { - n2end--; - n2_int_right_zeros++; + while (*divisorend == 0) { + divisorend--; + divisor_int_right_zeros++; } - if (*n1ptr == 0 && n1_len == 1) { - n1ptr++; - n1_len = 0; + if (*numeratorptr == 0 && numerator_len == 1) { + numeratorptr++; + numerator_len = 0; } - size_t n1_top_extension = 0; - size_t n1_bottom_extension = 0; - if (n2_scale > 0) { + size_t numerator_top_extension = 0; + size_t numerator_bottom_extension = 0; + if (divisor_scale > 0) { /* - * e.g. n2_scale = 4 - * n2 = .0002, to be 2 or n2 = 200.001, to be 200001 - * n1 = .03, to be 300 or n1 = .000003, to be .03 - * n1 may become longer than the original data length due to the addition of + * e.g. divisor_scale = 4 + * divisor = .0002, to be 2 or divisor = 200.001, to be 200001 + * numerator = .03, to be 300 or numerator = .000003, to be .03 + * numerator may become longer than the original data length due to the addition of * trailing zeros in the integer part. */ - n1_len += n2_scale; - n1_bottom_extension = n1_scale < n2_scale ? n2_scale - n1_scale : 0; - n1_scale = n1_scale > n2_scale ? n1_scale - n2_scale : 0; - n2_len += n2_scale; - n2_scale = 0; - } else if (n2_int_right_zeros > 0) { + numerator_len += divisor_scale; + numerator_bottom_extension = numerator_scale < divisor_scale ? divisor_scale - numerator_scale : 0; + numerator_scale = numerator_scale > divisor_scale ? numerator_scale - divisor_scale : 0; + divisor_len += divisor_scale; + divisor_scale = 0; + } else if (divisor_int_right_zeros > 0) { /* - * e.g. n2_int_right_zeros = 4 - * n2 = 2000, to be 2 - * n1 = 30, to be .03 or n1 = 30000, to be 30 - * Also, n1 may become longer than the original data length due to the addition of + * e.g. divisor_int_right_zeros = 4 + * divisor = 2000, to be 2 + * numerator = 30, to be .03 or numerator = 30000, to be 30 + * Also, numerator may become longer than the original data length due to the addition of * leading zeros in the fractional part. */ - n1_top_extension = n1_len < n2_int_right_zeros ? n2_int_right_zeros - n1_len : 0; - n1_len = n1_len > n2_int_right_zeros ? n1_len - n2_int_right_zeros : 0; - n1_scale += n2_int_right_zeros; - n2_len -= n2_int_right_zeros; - n2_scale = 0; + numerator_top_extension = numerator_len < divisor_int_right_zeros ? divisor_int_right_zeros - numerator_len : 0; + numerator_len = numerator_len > divisor_int_right_zeros ? numerator_len - divisor_int_right_zeros : 0; + numerator_scale += divisor_int_right_zeros; + divisor_len -= divisor_int_right_zeros; + divisor_scale = 0; } - /* remove n1 leading zeros */ - while (*n1ptr == 0 && n1_len > 0) { - n1ptr++; - n1_len--; + /* remove numerator leading zeros */ + while (*numeratorptr == 0 && numerator_len > 0) { + numeratorptr++; + numerator_len--; } - /* remove n2 leading zeros */ - while (*n2ptr == 0) { - n2ptr++; - n2_len--; + /* remove divisor leading zeros */ + while (*divisorptr == 0) { + divisorptr++; + divisor_len--; } /* Considering the scale specification, the quotient is always 0 if this condition is met */ - if (n2_len > n1_len + scale) { + if (divisor_len > numerator_len + scale) { *quot = bc_copy_num(BCG(_zero_)); return true; } - /* set scale to n1 */ - if (n1_scale > scale) { - size_t scale_diff = n1_scale - scale; - if (n1_bottom_extension > scale_diff) { - n1_bottom_extension -= scale_diff; + /* set scale to numerator */ + if (numerator_scale > scale) { + size_t scale_diff = numerator_scale - scale; + if (numerator_bottom_extension > scale_diff) { + numerator_bottom_extension -= scale_diff; } else { - n1_bottom_extension = 0; - n1end -= scale_diff - n1_bottom_extension; + numerator_bottom_extension = 0; + numeratorend -= scale_diff - numerator_bottom_extension; } } else { - n1_bottom_extension += scale - n1_scale; + numerator_bottom_extension += scale - numerator_scale; } - n1_scale = scale; + numerator_scale = scale; - /* Length of n1 data that can be read */ - size_t n1_readable_len = n1end - n1ptr + 1; + /* Length of numerator data that can be read */ + size_t numerator_readable_len = numeratorend - numeratorptr + 1; - /* If n2 is 1 here, return the result of adjusting the decimal point position of n1. */ - if (n2_len == 1 && *n2ptr == 1) { - if (n1_len == 0) { - n1_len = 1; - n1_top_extension++; + /* If divisor is 1 here, return the result of adjusting the decimal point position of numerator. */ + if (divisor_len == 1 && *divisorptr == 1) { + if (numerator_len == 0) { + numerator_len = 1; + numerator_top_extension++; } - size_t quot_scale = n1_scale > n1_bottom_extension ? n1_scale - n1_bottom_extension : 0; - n1_bottom_extension = n1_scale < n1_bottom_extension ? n1_bottom_extension - n1_scale : 0; + size_t quot_scale = numerator_scale > numerator_bottom_extension ? numerator_scale - numerator_bottom_extension : 0; + numerator_bottom_extension = numerator_scale < numerator_bottom_extension ? numerator_bottom_extension - numerator_scale : 0; - *quot = bc_new_num_nonzeroed(n1_len, quot_scale); + *quot = bc_new_num_nonzeroed(numerator_len, quot_scale); char *qptr = (*quot)->n_value; - for (size_t i = 0; i < n1_top_extension; i++) { + for (size_t i = 0; i < numerator_top_extension; i++) { *qptr++ = 0; } - memcpy(qptr, n1ptr, n1_readable_len); - qptr += n1_readable_len; - for (size_t i = 0; i < n1_bottom_extension; i++) { + memcpy(qptr, numeratorptr, numerator_readable_len); + qptr += numerator_readable_len; + for (size_t i = 0; i < numerator_bottom_extension; i++) { *qptr++ = 0; } - (*quot)->n_sign = n1->n_sign == n2->n_sign ? PLUS : MINUS; + (*quot)->n_sign = numerator->n_sign == divisor->n_sign ? PLUS : MINUS; return true; } size_t quot_full_len; - if (n2_len > n1_len) { + if (divisor_len > numerator_len) { *quot = bc_new_num_nonzeroed(1, scale); quot_full_len = 1 + scale; } else { - *quot = bc_new_num_nonzeroed(n1_len - n2_len + 1, scale); - quot_full_len = n1_len - n2_len + 1 + scale; + *quot = bc_new_num_nonzeroed(numerator_len - divisor_len + 1, scale); + quot_full_len = numerator_len - divisor_len + 1 + scale; } /* do divide */ - bc_do_div(n1end, n1_readable_len, n1_bottom_extension, n2end, n2_len, quot, quot_full_len); - (*quot)->n_sign = n1->n_sign == n2->n_sign ? PLUS : MINUS; + bc_do_div(numeratorend, numerator_readable_len, numerator_bottom_extension, divisorend, divisor_len, quot, quot_full_len); + (*quot)->n_sign = numerator->n_sign == divisor->n_sign ? PLUS : MINUS; _bc_rm_leading_zeros(*quot); return true; From d80c452d8cf84c9e116b5eac7003307b0596b36b Mon Sep 17 00:00:00 2001 From: Saki Takamachi Date: Mon, 22 Jul 2024 11:31:28 +0900 Subject: [PATCH 34/34] Added a early return case --- ext/bcmath/libbcmath/src/div.c | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/ext/bcmath/libbcmath/src/div.c b/ext/bcmath/libbcmath/src/div.c index f2465b639cac6..14477f35fb391 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -343,6 +343,17 @@ bool bc_divide(bc_num numerator, bc_num divisor, bc_num *quot, size_t scale) return true; } + /* If divisor is 1 / -1, the quotient's n_value is equal to numerator's n_value. */ + if (_bc_do_compare(divisor, BCG(_one_), scale, false) == BCMATH_EQUAL) { + size_t quot_scale = MIN(numerator->n_scale, scale); + *quot = bc_new_num_nonzeroed(numerator->n_len, quot_scale); + char *qptr = (*quot)->n_value; + memcpy(qptr, numerator->n_value, numerator->n_len + quot_scale); + (*quot)->n_sign = numerator->n_sign == divisor->n_sign ? PLUS : MINUS; + _bc_rm_leading_zeros(*quot); + return true; + } + char *numeratorptr = numerator->n_value; char *numeratorend = numeratorptr + numerator->n_len + numerator->n_scale - 1; size_t numerator_len = numerator->n_len;