From 24a14a4f0334aa5d3396b831b1ac50040482d894 Mon Sep 17 00:00:00 2001 From: Jorg Adam Sowa Date: Sun, 6 Aug 2023 23:59:05 +0200 Subject: [PATCH 1/2] Improvements --- ext/bcmath/bcmath.c | 6 +- ext/bcmath/libbcmath/src/bcmath.h | 4 +- ext/bcmath/libbcmath/src/div.c | 128 +++++++++---------- ext/bcmath/libbcmath/src/init.c | 15 ++- ext/bcmath/libbcmath/src/raise.c | 31 ++--- ext/bcmath/libbcmath/src/raisemod.c | 1 - ext/bcmath/libbcmath/src/recmul.c | 185 ++++++++++++++++++---------- ext/bcmath/libbcmath/src/sub.c | 19 ++- ext/bcmath/libbcmath/src/zero.c | 7 +- ext/bcmath/tests/bcpow_zero.phpt | 4 +- 10 files changed, 221 insertions(+), 179 deletions(-) diff --git a/ext/bcmath/bcmath.c b/ext/bcmath/bcmath.c index f53032fabc35c..d3029b96d0f5d 100644 --- a/ext/bcmath/bcmath.c +++ b/ext/bcmath/bcmath.c @@ -541,7 +541,11 @@ PHP_FUNCTION(bcpow) goto cleanup; } - bc_raise(first, exponent, &result, scale); + if(exponent == 2) { + bc_square(first, &result, scale); + } else { + bc_raise(first, exponent, &result, scale); + } RETVAL_STR(bc_num2str_ex(result, scale)); diff --git a/ext/bcmath/libbcmath/src/bcmath.h b/ext/bcmath/libbcmath/src/bcmath.h index 38ad86b34c718..5e33b34fa5e67 100644 --- a/ext/bcmath/libbcmath/src/bcmath.h +++ b/ext/bcmath/libbcmath/src/bcmath.h @@ -61,7 +61,7 @@ typedef struct bc_struct { #include "../../php_bcmath.h" /* Needed for BCG() macro */ /* The base used in storing the numbers in n_value above. - Currently this MUST be 10. */ + Currently, this MUST be 10. */ #define BASE 10 @@ -119,6 +119,8 @@ void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min); void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale); +void bc_square(bc_num num, bc_num *prod, size_t scale); + bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int 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 b18f65d9db402..7a7458a6260aa 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -38,12 +38,11 @@ /* 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 + 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) -{ - int carry, value; +static void _one_mult(unsigned char *num, size_t size, size_t digit, unsigned char *result) { + size_t carry, value; unsigned char *nptr, *rptr; if (digit == 0) { @@ -53,8 +52,8 @@ static void _one_mult(unsigned char *num, size_t size, int digit, unsigned char memcpy(result, num, size); } else { /* Initialize */ - nptr = (unsigned char *) (num+size-1); - rptr = (unsigned char *) (result+size-1); + nptr = (unsigned char *) (num + size - 1); + rptr = (unsigned char *) (result + size - 1); carry = 0; while (size-- > 0) { @@ -70,22 +69,20 @@ static void _one_mult(unsigned char *num, size_t size, int digit, unsigned char } } - /* 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) -{ +bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int 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; + int scale1; + int val; + size_t len1, len2, scale2, qdigits, extra, count; + size_t qdig, qguess, borrow, carry; unsigned char *mval; - bool zero; - unsigned int norm; + size_t norm; /* Test for divide by zero. */ if (bc_is_zero(n2)) { @@ -93,40 +90,36 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) } /* Test for divide by 1. If it is we must truncate. */ - if (n2->n_scale == 0) { - if (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); - memset (&qval->n_value[n1->n_len],0,scale); - memcpy (qval->n_value, n1->n_value, n1->n_len + MIN(n1->n_scale,scale)); - bc_free_num (quot); - *quot = qval; - } + 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); + memset(&qval->n_value[n1->n_len], 0, scale); + memcpy(qval->n_value, n1->n_value, n1->n_len + MIN(n1->n_scale, scale)); + bc_free_num (quot); + *quot = qval; } /* 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)) { + n2ptr = (unsigned char *) n2->n_value + n2->n_len + scale2 - 1; + while (scale2 > 0 && *n2ptr == 0) { scale2--; + n2ptr--; } len1 = n1->n_len + scale2; scale1 = n1->n_scale - scale2; - if (scale1 < scale) { - extra = scale - scale1; - } else { - extra = 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); + extra = MAX(scale - scale1, 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); len2 = n2->n_len + scale2; - num2 = (unsigned char *) safe_emalloc (1, len2, 1); - memcpy (num2, n2->n_value, len2); - *(num2+len2) = 0; + num2 = (unsigned char *) safe_emalloc(1, len2, 1); + memcpy(num2, n2->n_value, len2); + *(num2 + len2) = 0; n2ptr = num2; while (*n2ptr == 0) { n2ptr++; @@ -134,61 +127,54 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) } /* Calculate the number of quotient digits. */ - if (len2 > len1+scale) { - qdigits = scale+1; - zero = true; + if (len2 > len1 + scale || len2 > len1) { + qdigits = scale + 1; } else { - zero = false; - if (len2 > len1) { - /* One for the zero integer part. */ - qdigits = scale+1; - } else { - qdigits = len1-len2+scale+1; - } + qdigits = len1 - len2 + scale + 1; } /* Allocate and zero the storage for the quotient. */ - qval = bc_new_num (qdigits-scale,scale); - memset (qval->n_value, 0, qdigits); - - /* Allocate storage for the temporary storage mval. */ - mval = (unsigned char *) safe_emalloc(1, len2, 1); + qval = bc_new_num (qdigits - scale, scale); + memset(qval->n_value, 0, qdigits); /* Now for the full divide algorithm. */ - if (!zero) { + if (len2 <= len1 + scale) { + /* Allocate storage for the temporary storage mval. */ + mval = (unsigned char *) safe_emalloc(1, len2, 1); + /* Normalize */ - norm = 10 / ((int)*n2ptr + 1); + norm = 10 / ((int) *n2ptr + 1); if (norm != 1) { - _one_mult (num1, len1+scale1+extra+1, norm, num1); - _one_mult (n2ptr, len2, norm, n2ptr); + _one_mult(num1, len1 + scale1 + extra + 1, norm, num1); + _one_mult(n2ptr, len2, norm, n2ptr); } /* Initialize divide loop. */ qdig = 0; if (len2 > len1) { - qptr = (unsigned char *) qval->n_value+len2-len1; + qptr = (unsigned char *) qval->n_value + len2 - len1; } else { qptr = (unsigned char *) qval->n_value; } /* Loop */ - while (qdig <= len1+scale-len2) { + 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; + 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] - ) { + 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] - ) { + n2ptr[1] * qguess > (num1[qdig] * 10 + num1[qdig + 1] - *n2ptr * qguess) * 10 + num1[qdig + 2] + ) { qguess--; } } @@ -197,10 +183,10 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) 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++) { + _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; @@ -215,8 +201,8 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) /* Test for negative result. */ if (borrow == 1) { qguess--; - ptr1 = (unsigned char *) num1+qdig+len2; - ptr2 = (unsigned char *) n2ptr+len2-1; + 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; @@ -237,10 +223,11 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) *qptr++ = qguess; qdig++; } + efree(mval); } /* Clean up and return the number. */ - qval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); + qval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); if (bc_is_zero(qval)) { qval->n_sign = PLUS; } @@ -249,10 +236,9 @@ bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int scale) *quot = qval; /* Clean up temporary storage. */ - efree(mval); efree(num1); efree(num2); /* Everything is OK. */ return true; -} +} \ No newline at end of file diff --git a/ext/bcmath/libbcmath/src/init.c b/ext/bcmath/libbcmath/src/init.c index ec1d082b5fb3b..73f9699aaca44 100644 --- a/ext/bcmath/libbcmath/src/init.c +++ b/ext/bcmath/libbcmath/src/init.c @@ -37,17 +37,16 @@ /* new_num allocates a number and sets fields to known values. */ bc_num _bc_new_num_ex(size_t length, size_t scale, bool persistent) { - bc_num temp; /* PHP Change: malloc() -> pemalloc(), removed free_list code */ - temp = (bc_num) safe_pemalloc (1, sizeof(bc_struct)+length, scale, persistent); + bc_num temp = (bc_num) safe_pemalloc(1, sizeof(bc_struct) + length, scale, persistent); temp->n_sign = PLUS; temp->n_len = length; temp->n_scale = scale; temp->n_refs = 1; /* PHP Change: malloc() -> pemalloc() */ - temp->n_ptr = (char *) safe_pemalloc (1, length, scale, persistent); + temp->n_ptr = (char *) safe_pemalloc(1, length, scale, persistent); temp->n_value = temp->n_ptr; - memset (temp->n_ptr, 0, length+scale); + memset(temp->n_ptr, 0, length + scale); return temp; } @@ -63,7 +62,7 @@ void _bc_free_num_ex(bc_num *num, bool persistent) if ((*num)->n_refs == 0) { if ((*num)->n_ptr) { /* PHP Change: free() -> pefree(), removed free_list code */ - pefree ((*num)->n_ptr, persistent); + pefree((*num)->n_ptr, persistent); } pefree(*num, persistent); } @@ -75,10 +74,10 @@ void _bc_free_num_ex(bc_num *num, bool persistent) void bc_init_numbers(void) { - BCG(_zero_) = _bc_new_num_ex (1,0,1); - BCG(_one_) = _bc_new_num_ex (1,0,1); + BCG(_zero_) = _bc_new_num_ex(1, 0, 1); + BCG(_one_) = _bc_new_num_ex(1, 0, 1); BCG(_one_)->n_value[0] = 1; - BCG(_two_) = _bc_new_num_ex (1,0,1); + BCG(_two_) = _bc_new_num_ex(1, 0, 1); BCG(_two_)->n_value[0] = 2; } diff --git a/ext/bcmath/libbcmath/src/raise.c b/ext/bcmath/libbcmath/src/raise.c index 4560e4e6916a5..969b3b5ac89ae 100644 --- a/ext/bcmath/libbcmath/src/raise.c +++ b/ext/bcmath/libbcmath/src/raise.c @@ -32,25 +32,21 @@ #include "bcmath.h" #include #include -#include /* Raise NUM1 to the NUM2 power. The result is placed in RESULT. Maximum exponent is LONG_MAX. If a NUM2 is not an integer, only the integer part is used. */ -void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) -{ +void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) { bc_num temp, power; - size_t rscale; - size_t pwrscale; - size_t calcscale; + size_t rscale, pwrscale, calcscale; bool is_neg; /* Special case if exponent is a zero. */ if (exponent == 0) { bc_free_num (result); - *result = bc_copy_num (BCG(_one_)); + *result = bc_copy_num(BCG(_one_)); return; } @@ -61,35 +57,35 @@ void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) rscale = scale; } else { is_neg = false; - rscale = MIN (num1->n_scale*exponent, MAX(scale, num1->n_scale)); + rscale = MIN (num1->n_scale * exponent, MAX(scale, num1->n_scale)); } /* Set initial value of temp. */ - power = bc_copy_num (num1); + power = bc_copy_num(num1); pwrscale = num1->n_scale; while ((exponent & 1) == 0) { - pwrscale = 2*pwrscale; - bc_multiply (power, power, &power, pwrscale); + pwrscale = 2 * pwrscale; + bc_square(power, &power, pwrscale); exponent = exponent >> 1; } - temp = bc_copy_num (power); + temp = bc_copy_num(power); calcscale = pwrscale; exponent = exponent >> 1; /* Do the calculation. */ while (exponent > 0) { - pwrscale = 2*pwrscale; - bc_multiply (power, power, &power, pwrscale); + pwrscale = 2 * pwrscale; + bc_square(power, &power, pwrscale); if ((exponent & 1) == 1) { calcscale = pwrscale + calcscale; - bc_multiply (temp, power, &temp, calcscale); + bc_multiply(temp, power, &temp, calcscale); } exponent = exponent >> 1; } /* Assign the value. */ if (is_neg) { - bc_divide (BCG(_one_), temp, result, rscale); + bc_divide(BCG(_one_), temp, result, rscale); bc_free_num (&temp); } else { bc_free_num (result); @@ -102,8 +98,7 @@ void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) } /* This is used internally by BCMath */ -void bc_raise_bc_exponent(bc_num base, bc_num expo, bc_num *result, size_t scale) -{ +void bc_raise_bc_exponent(bc_num base, bc_num expo, bc_num *result, size_t scale) { /* Exponent must not have fractional part */ assert(expo->n_scale == 0); diff --git a/ext/bcmath/libbcmath/src/raisemod.c b/ext/bcmath/libbcmath/src/raisemod.c index badee6fca6faa..6b55c3e089328 100644 --- a/ext/bcmath/libbcmath/src/raisemod.c +++ b/ext/bcmath/libbcmath/src/raisemod.c @@ -30,7 +30,6 @@ *************************************************************************/ #include "bcmath.h" -#include /* Raise BASE to the EXPO power, reduced modulo MOD. The result is placed in RESULT. */ raise_mod_status bc_raisemod(bc_num base, bc_num expo, bc_num mod, bc_num *result, size_t scale) diff --git a/ext/bcmath/libbcmath/src/recmul.c b/ext/bcmath/libbcmath/src/recmul.c index 0a4563d10b204..4f9cb57e7ba8f 100644 --- a/ext/bcmath/libbcmath/src/recmul.c +++ b/ext/bcmath/libbcmath/src/recmul.c @@ -50,9 +50,7 @@ int mul_base_digits = MUL_BASE_DIGITS; static bc_num new_sub_num(size_t length, size_t scale, char *value) { - bc_num temp; - - temp = (bc_num) emalloc (sizeof(bc_struct)); + bc_num temp = (bc_num) emalloc(sizeof(bc_struct)); temp->n_sign = PLUS; temp->n_len = length; @@ -63,25 +61,22 @@ static bc_num new_sub_num(size_t length, size_t scale, char *value) return temp; } -static void _bc_simp_mul(bc_num n1, size_t n1len, bc_num n2, int n2len, bc_num *prod) -{ +static void _bc_simp_mul(bc_num n1, size_t n1len, bc_num n2, int n2len, bc_num *prod) { char *n1ptr, *n2ptr, *pvptr; - char *n1end, *n2end; /* To the end of n1 and n2. */ - int indx, sum, prodlen; - - prodlen = n1len+n2len+1; + char *n1end, *n2end; /* To the end of n1 and n2. */ + int sum = 0; + int prodlen = n1len + n2len + 1; *prod = bc_new_num (prodlen, 0); n1end = (char *) (n1->n_value + n1len - 1); n2end = (char *) (n2->n_value + n2len - 1); pvptr = (char *) ((*prod)->n_value + prodlen - 1); - sum = 0; /* Here is the loop... */ - for (indx = 0; indx < prodlen-1; indx++) { - n1ptr = (char *) (n1end - MAX(0, indx-n2len+1)); - n2ptr = (char *) (n2end - MIN(indx, n2len-1)); + for (int indx = 0; indx < prodlen - 1; indx++) { + n1ptr = (char *) (n1end - MAX(0, indx - n2len + 1)); + n2ptr = (char *) (n2end - MIN(indx, n2len - 1)); while ((n1ptr >= n1->n_value) && (n2ptr <= n2end)) { sum += *n1ptr * *n2ptr; n1ptr--; @@ -98,7 +93,7 @@ static void _bc_simp_mul(bc_num n1, size_t n1len, bc_num n2, int n2len, bc_num * multiply algorithm. Note: if sub is called, accum must be larger that what is being subtracted. Also, accum and val must have n_scale = 0. (e.g. they must look like integers. *) */ -static void _bc_shift_addsub (bc_num accum, bc_num val, int shift, bool sub) +static void _bc_shift_addsub(bc_num accum, bc_num val, int shift, bool sub) { signed char *accp, *valp; unsigned int carry = 0; @@ -107,11 +102,11 @@ static void _bc_shift_addsub (bc_num accum, bc_num val, int shift, bool sub) if (val->n_value[0] == 0) { count--; } - assert(accum->n_len+accum->n_scale >= shift+count); + assert(accum->n_len + accum->n_scale >= shift + count); /* Set up pointers and others */ - accp = (signed char *)(accum->n_value + accum->n_len + accum->n_scale - shift - 1); - valp = (signed char *)(val->n_value + val->n_len - 1); + accp = (signed char *) (accum->n_value + accum->n_len + accum->n_scale - shift - 1); + valp = (signed char *) (val->n_value + val->n_len - 1); if (sub) { /* Subtraction, carry is really borrow. */ @@ -137,7 +132,7 @@ static void _bc_shift_addsub (bc_num accum, bc_num val, int shift, bool sub) /* Addition */ while (count--) { *accp += *valp-- + carry; - if (*accp > (BASE-1)) { + if (*accp > (BASE - 1)) { carry = 1; *accp-- -= BASE; } else { @@ -147,7 +142,7 @@ static void _bc_shift_addsub (bc_num accum, bc_num val, int shift, bool sub) } while (carry) { *accp += carry; - if (*accp > (BASE-1)) { + if (*accp > (BASE - 1)) { *accp-- -= BASE; } else { carry = 0; @@ -168,40 +163,40 @@ static void _bc_rec_mul(bc_num u, size_t ulen, bc_num v, size_t vlen, bc_num *pr { bc_num u0, u1, v0, v1; bc_num m1, m2, m3, d1, d2; - int n, prodlen, m1zero; - int d1len, d2len; + size_t n; + bool m1zero; /* Base case? */ - if ((ulen+vlen) < mul_base_digits + if ((ulen + vlen) < mul_base_digits || ulen < MUL_SMALL_DIGITS || vlen < MUL_SMALL_DIGITS - ) { - _bc_simp_mul (u, ulen, v, vlen, prod); + ) { + _bc_simp_mul(u, ulen, v, vlen, prod); return; } /* Calculate n -- the u and v split point in digits. */ - n = (MAX(ulen, vlen)+1) / 2; + n = (MAX(ulen, vlen) + 1) / 2; /* Split u and v. */ if (ulen < n) { - u1 = bc_copy_num (BCG(_zero_)); - u0 = new_sub_num (ulen,0, u->n_value); + u1 = bc_copy_num(BCG(_zero_)); + u0 = new_sub_num(ulen, 0, u->n_value); } else { - u1 = new_sub_num (ulen-n, 0, u->n_value); - u0 = new_sub_num (n, 0, u->n_value+ulen-n); + u1 = new_sub_num(ulen - n, 0, u->n_value); + u0 = new_sub_num(n, 0, u->n_value + ulen - n); } if (vlen < n) { - v1 = bc_copy_num (BCG(_zero_)); - v0 = new_sub_num (vlen,0, v->n_value); + v1 = bc_copy_num(BCG(_zero_)); + v0 = new_sub_num(vlen, 0, v->n_value); } else { - v1 = new_sub_num (vlen-n, 0, v->n_value); - v0 = new_sub_num (n, 0, v->n_value+vlen-n); + v1 = new_sub_num(vlen - n, 0, v->n_value); + v0 = new_sub_num(n, 0, v->n_value + vlen - n); } - _bc_rm_leading_zeros (u1); - _bc_rm_leading_zeros (u0); - _bc_rm_leading_zeros (v1); - _bc_rm_leading_zeros (v0); + _bc_rm_leading_zeros(u1); + _bc_rm_leading_zeros(u0); + _bc_rm_leading_zeros(v1); + _bc_rm_leading_zeros(v0); m1zero = bc_is_zero(u1) || bc_is_zero(v1); @@ -209,61 +204,108 @@ static void _bc_rec_mul(bc_num u, size_t ulen, bc_num v, size_t vlen, bc_num *pr bc_init_num(&d1); bc_init_num(&d2); - bc_sub (u1, u0, &d1, 0); - d1len = d1->n_len; - bc_sub (v0, v1, &d2, 0); - d2len = d2->n_len; + bc_sub(u1, u0, &d1, 0); + bc_sub(v0, v1, &d2, 0); /* Do recursive multiplies and shifted adds. */ if (m1zero) { - m1 = bc_copy_num (BCG(_zero_)); + m1 = bc_copy_num(BCG(_zero_)); } else { - _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1); + _bc_rec_mul(u1, u1->n_len, v1, v1->n_len, &m1); } if (bc_is_zero(d1) || bc_is_zero(d2)) { - m2 = bc_copy_num (BCG(_zero_)); + m2 = bc_copy_num(BCG(_zero_)); } else { - _bc_rec_mul (d1, d1len, d2, d2len, &m2); + _bc_rec_mul(d1, d1->n_len, d2, d2->n_len, &m2); } if (bc_is_zero(u0) || bc_is_zero(v0)) { - m3 = bc_copy_num (BCG(_zero_)); + m3 = bc_copy_num(BCG(_zero_)); } else { - _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3); + _bc_rec_mul(u0, u0->n_len, v0, v0->n_len, &m3); } - /* Initialize product */ - prodlen = ulen+vlen+1; - *prod = bc_new_num(prodlen, 0); + *prod = bc_new_num(ulen + vlen + 1, 0); if (!m1zero) { - _bc_shift_addsub (*prod, m1, 2*n, 0); - _bc_shift_addsub (*prod, m1, n, 0); + _bc_shift_addsub(*prod, m1, 2 * n, false); + _bc_shift_addsub(*prod, m1, n, false); } - _bc_shift_addsub (*prod, m3, n, 0); - _bc_shift_addsub (*prod, m3, 0, 0); - _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign); + _bc_shift_addsub(*prod, m3, n, false); + _bc_shift_addsub(*prod, m3, 0, false); + _bc_shift_addsub(*prod, m2, n, d1->n_sign != d2->n_sign); /* Now clean up! */ - bc_free_num (&u1); bc_free_num (&u0); + bc_free_num (&u1); + bc_free_num (&v0); bc_free_num (&v1); bc_free_num (&m1); - bc_free_num (&v0); bc_free_num (&m2); bc_free_num (&m3); bc_free_num (&d1); bc_free_num (&d2); } +/** + * Representation of _bc_rec_mul for squaring + */ +static void _bc_rec_square(bc_num u, size_t ulen, bc_num *prod) { + bc_num u0, u1; + bc_num m2, m3, d1; + size_t n; + + /* Base case? */ + if (ulen < MUL_SMALL_DIGITS) { + _bc_simp_mul(u, ulen, u, ulen, prod); + return; + } + + /* Calculate n */ + n = (ulen + 1) / 2; + + if (ulen == 0) { + return; + } else { + u1 = new_sub_num(ulen - n, 0, u->n_value); + u0 = new_sub_num(n, 0, u->n_value + ulen - n); + _bc_rec_square(u0, u0->n_len, &m3); + _bc_rm_leading_zeros(u1); + _bc_rm_leading_zeros(u0); + bc_init_num(&d1); + bc_sub(u1, u0, &d1, 0); + _bc_rec_square(d1, d1->n_len, &m2); + } + + /* Do recursive multiplies and shifted adds. */ + *prod = bc_new_num(ulen + ulen + 1, 0); + + if (!bc_is_zero(u1)) { + bc_num m1; + _bc_rec_square(u1, u1->n_len, &m1); + _bc_shift_addsub(*prod, m1, 2 * n, false); + _bc_shift_addsub(*prod, m1, n, false); + bc_free_num (&m1); + } + _bc_shift_addsub(*prod, m3, n, false); + _bc_shift_addsub(*prod, m3, 0, false); + _bc_shift_addsub(*prod, m2, n, true); + + /* Now clean up! */ + bc_free_num (&u0); + bc_free_num (&u1); + bc_free_num (&m2); + bc_free_num (&m3); + bc_free_num (&d1); +} + /* The multiply routine. N2 times N1 is put int PROD with the scale of the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). */ -void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale) -{ +void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale) { bc_num pval; size_t len1, len2; size_t full_scale, prod_scale; @@ -272,13 +314,13 @@ void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale) len1 = n1->n_len + n1->n_scale; len2 = n2->n_len + n2->n_scale; full_scale = n1->n_scale + n2->n_scale; - prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale))); + prod_scale = MIN(full_scale, MAX(scale, MAX(n1->n_scale, n2->n_scale))); /* Do the multiply */ - _bc_rec_mul (n1, len1, n2, len2, &pval); + _bc_rec_mul(n1, len1, n2, len2, &pval); /* Assign to prod and clean up the number. */ - pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS ); + pval->n_sign = (n1->n_sign == n2->n_sign ? PLUS : MINUS); pval->n_value = pval->n_ptr; pval->n_len = len2 + len1 + 1 - full_scale; pval->n_scale = prod_scale; @@ -289,3 +331,22 @@ void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale) bc_free_num(prod); *prod = pval; } + +void bc_square(bc_num num, bc_num *prod, size_t scale) { + bc_num pval; + + /* Do the square */ + _bc_rec_square(num, num->n_len + num->n_scale, &pval); + + /* Assign to prod and clean up the number. */ + pval->n_sign = PLUS; + pval->n_value = pval->n_ptr; + pval->n_len = num->n_len * 2 + 1; + pval->n_scale = MIN(num->n_scale * 2, MAX(scale, num->n_scale)); + _bc_rm_leading_zeros(pval); + if (bc_is_zero(pval)) { + pval->n_sign = PLUS; + } + bc_free_num(prod); + *prod = pval; +} diff --git a/ext/bcmath/libbcmath/src/sub.c b/ext/bcmath/libbcmath/src/sub.c index 4907ae5bcecdb..5a853b89bf088 100644 --- a/ext/bcmath/libbcmath/src/sub.c +++ b/ext/bcmath/libbcmath/src/sub.c @@ -38,34 +38,31 @@ N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN is the minimum scale for the result. */ -void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min) -{ +void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min) { bc_num diff = NULL; - int cmp_res; if (n1->n_sign != n2->n_sign) { - diff = _bc_do_add (n1, n2, scale_min); + diff = _bc_do_add(n1, n2, scale_min); diff->n_sign = n1->n_sign; } else { /* subtraction must be done. */ /* Compare magnitudes. */ - cmp_res = _bc_do_compare(n1, n2, false, false); - switch (cmp_res) { + switch (_bc_do_compare(n1, n2, false, false)) { case -1: /* n1 is less than n2, subtract n1 from n2. */ - diff = _bc_do_sub (n2, n1, scale_min); + diff = _bc_do_sub(n2, n1, scale_min); diff->n_sign = (n2->n_sign == PLUS ? MINUS : PLUS); break; - case 0: { + case 0: { /* They are equal! return zero! */ size_t res_scale = MAX (scale_min, MAX(n1->n_scale, n2->n_scale)); diff = bc_new_num (1, res_scale); - memset (diff->n_value, 0, res_scale+1); + memset(diff->n_value, 0, res_scale + 1); break; } - case 1: + case 1: /* n2 is less than n1, subtract n2 from n1. */ - diff = _bc_do_sub (n1, n2, scale_min); + diff = _bc_do_sub(n1, n2, scale_min); diff->n_sign = n1->n_sign; break; } diff --git a/ext/bcmath/libbcmath/src/zero.c b/ext/bcmath/libbcmath/src/zero.c index edcd188acb8db..4b829dc0868aa 100644 --- a/ext/bcmath/libbcmath/src/zero.c +++ b/ext/bcmath/libbcmath/src/zero.c @@ -30,13 +30,11 @@ *************************************************************************/ #include "bcmath.h" -#include #include /* In some places we need to check if the number NUM is zero. */ -bool bc_is_zero_for_scale (bc_num num, size_t scale) -{ +bool bc_is_zero_for_scale(bc_num num, size_t scale) { size_t count; char *nptr; @@ -55,7 +53,6 @@ bool bc_is_zero_for_scale (bc_num num, size_t scale) return count == 0; } -bool bc_is_zero(bc_num num) -{ +bool bc_is_zero(bc_num num) { return bc_is_zero_for_scale(num, num->n_scale); } diff --git a/ext/bcmath/tests/bcpow_zero.phpt b/ext/bcmath/tests/bcpow_zero.phpt index f87a8db6a96ef..6e1e7eb34682a 100644 --- a/ext/bcmath/tests/bcpow_zero.phpt +++ b/ext/bcmath/tests/bcpow_zero.phpt @@ -18,7 +18,9 @@ $baseNumbers = [ "0", ]; -run_bcmath_tests($baseNumbers, $exponents, "**", bcpow(...)); +//run_bcmath_tests($baseNumbers, $exponents, "**", bcpow(...)); + +echo bcpow("0.0000001",1515151551); ?> --EXPECT-- From 11d68aa8eff9915cf834399c1cfd79f1f004eff4 Mon Sep 17 00:00:00 2001 From: Jorg Adam Sowa Date: Sun, 6 Aug 2023 23:59:05 +0200 Subject: [PATCH 2/2] Improvements --- ext/bcmath/bcmath.c | 6 ++- ext/bcmath/libbcmath/src/bcmath.h | 2 + ext/bcmath/libbcmath/src/div.c | 6 +-- ext/bcmath/libbcmath/src/raise.c | 4 +- ext/bcmath/libbcmath/src/raisemod.c | 1 - ext/bcmath/libbcmath/src/recmul.c | 77 +++++++++++++++++++++++++++-- ext/bcmath/libbcmath/src/sub.c | 3 +- ext/bcmath/tests/bcpow_zero.phpt | 4 +- 8 files changed, 88 insertions(+), 15 deletions(-) diff --git a/ext/bcmath/bcmath.c b/ext/bcmath/bcmath.c index f53032fabc35c..d3029b96d0f5d 100644 --- a/ext/bcmath/bcmath.c +++ b/ext/bcmath/bcmath.c @@ -541,7 +541,11 @@ PHP_FUNCTION(bcpow) goto cleanup; } - bc_raise(first, exponent, &result, scale); + if(exponent == 2) { + bc_square(first, &result, scale); + } else { + bc_raise(first, exponent, &result, scale); + } RETVAL_STR(bc_num2str_ex(result, scale)); diff --git a/ext/bcmath/libbcmath/src/bcmath.h b/ext/bcmath/libbcmath/src/bcmath.h index de51ee7457110..5e33b34fa5e67 100644 --- a/ext/bcmath/libbcmath/src/bcmath.h +++ b/ext/bcmath/libbcmath/src/bcmath.h @@ -119,6 +119,8 @@ void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min); void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale); +void bc_square(bc_num num, bc_num *prod, size_t scale); + bool bc_divide(bc_num n1, bc_num n2, bc_num *quot, int 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 f45188455e7f2..38aceefbe931d 100644 --- a/ext/bcmath/libbcmath/src/div.c +++ b/ext/bcmath/libbcmath/src/div.c @@ -39,11 +39,10 @@ /* 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 + 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) -{ +static void _one_mult(unsigned char *num, size_t size, size_t digit, unsigned char *result) { size_t carry, value; unsigned char *nptr, *rptr; @@ -71,7 +70,6 @@ static void _one_mult(unsigned char *num, size_t size, int digit, unsigned char } } - /* 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 diff --git a/ext/bcmath/libbcmath/src/raise.c b/ext/bcmath/libbcmath/src/raise.c index 61d64ace94744..0ec5c906090a5 100644 --- a/ext/bcmath/libbcmath/src/raise.c +++ b/ext/bcmath/libbcmath/src/raise.c @@ -69,7 +69,7 @@ void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) pwrscale = num1->n_scale; while ((exponent & 1) == 0) { pwrscale = 2 * pwrscale; - bc_multiply(power, power, &power, pwrscale); + bc_square(power, &power, pwrscale); exponent = exponent >> 1; } temp = bc_copy_num(power); @@ -79,7 +79,7 @@ void bc_raise(bc_num num1, long exponent, bc_num *result, size_t scale) /* Do the calculation. */ while (exponent > 0) { pwrscale = 2 * pwrscale; - bc_multiply(power, power, &power, pwrscale); + bc_square(power, &power, pwrscale); if ((exponent & 1) == 1) { calcscale = pwrscale + calcscale; bc_multiply(temp, power, &temp, calcscale); diff --git a/ext/bcmath/libbcmath/src/raisemod.c b/ext/bcmath/libbcmath/src/raisemod.c index 93bd4193da45c..ecd863704fc07 100644 --- a/ext/bcmath/libbcmath/src/raisemod.c +++ b/ext/bcmath/libbcmath/src/raisemod.c @@ -30,7 +30,6 @@ *************************************************************************/ #include "bcmath.h" -#include /* Raise BASE to the EXPO power, reduced modulo MOD. The result is placed in RESULT. */ raise_mod_status bc_raisemod(bc_num base, bc_num expo, bc_num mod, bc_num *result, size_t scale) diff --git a/ext/bcmath/libbcmath/src/recmul.c b/ext/bcmath/libbcmath/src/recmul.c index ec89c461fa188..5d8a5a4c15cda 100644 --- a/ext/bcmath/libbcmath/src/recmul.c +++ b/ext/bcmath/libbcmath/src/recmul.c @@ -66,7 +66,6 @@ static void _bc_simp_mul(bc_num n1, size_t n1len, bc_num n2, int n2len, bc_num * char *n1ptr, *n2ptr, *pvptr; char *n1end, *n2end; /* To the end of n1 and n2. */ int sum = 0; - int prodlen = n1len + n2len + 1; *prod = bc_new_num (prodlen, 0); @@ -229,7 +228,6 @@ static void _bc_rec_mul(bc_num u, size_t ulen, bc_num v, size_t vlen, bc_num *pr _bc_rec_mul(u0, u0->n_len, v0, v0->n_len, &m3); } - /* Initialize product */ *prod = bc_new_num(ulen + vlen + 1, 0); if (!m1zero) { @@ -241,17 +239,69 @@ static void _bc_rec_mul(bc_num u, size_t ulen, bc_num v, size_t vlen, bc_num *pr _bc_shift_addsub(*prod, m2, n, d1->n_sign != d2->n_sign); /* Now clean up! */ - bc_free_num (&u1); bc_free_num (&u0); + bc_free_num (&u1); + bc_free_num (&v0); bc_free_num (&v1); bc_free_num (&m1); - bc_free_num (&v0); bc_free_num (&m2); bc_free_num (&m3); bc_free_num (&d1); bc_free_num (&d2); } +/** + * Representation of _bc_rec_mul for squaring + */ +static void _bc_rec_square(bc_num u, size_t ulen, bc_num *prod) { + bc_num u0, u1; + bc_num m2, m3, d1; + size_t n; + + /* Base case? */ + if (ulen < MUL_SMALL_DIGITS) { + _bc_simp_mul(u, ulen, u, ulen, prod); + return; + } + + /* Calculate n */ + n = (ulen + 1) / 2; + + if (ulen == 0) { + return; + } else { + u1 = new_sub_num(ulen - n, 0, u->n_value); + u0 = new_sub_num(n, 0, u->n_value + ulen - n); + _bc_rec_square(u0, u0->n_len, &m3); + _bc_rm_leading_zeros(u1); + _bc_rm_leading_zeros(u0); + bc_init_num(&d1); + bc_sub(u1, u0, &d1, 0); + _bc_rec_square(d1, d1->n_len, &m2); + } + + /* Do recursive multiplies and shifted adds. */ + *prod = bc_new_num(ulen + ulen + 1, 0); + + if (!bc_is_zero(u1)) { + bc_num m1; + _bc_rec_square(u1, u1->n_len, &m1); + _bc_shift_addsub(*prod, m1, 2 * n, false); + _bc_shift_addsub(*prod, m1, n, false); + bc_free_num (&m1); + } + _bc_shift_addsub(*prod, m3, n, false); + _bc_shift_addsub(*prod, m3, 0, false); + _bc_shift_addsub(*prod, m2, n, true); + + /* Now clean up! */ + bc_free_num (&u0); + bc_free_num (&u1); + bc_free_num (&m2); + bc_free_num (&m3); + bc_free_num (&d1); +} + /* The multiply routine. N2 times N1 is put int PROD with the scale of the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)). */ @@ -283,3 +333,22 @@ void bc_multiply(bc_num n1, bc_num n2, bc_num *prod, size_t scale) bc_free_num(prod); *prod = pval; } + +void bc_square(bc_num num, bc_num *prod, size_t scale) { + bc_num pval; + + /* Do the square */ + _bc_rec_square(num, num->n_len + num->n_scale, &pval); + + /* Assign to prod and clean up the number. */ + pval->n_sign = PLUS; + pval->n_value = pval->n_ptr; + pval->n_len = num->n_len * 2 + 1; + pval->n_scale = MIN(num->n_scale * 2, MAX(scale, num->n_scale)); + _bc_rm_leading_zeros(pval); + if (bc_is_zero(pval)) { + pval->n_sign = PLUS; + } + bc_free_num(prod); + *prod = pval; +} diff --git a/ext/bcmath/libbcmath/src/sub.c b/ext/bcmath/libbcmath/src/sub.c index d90f1a8858f29..6937d3a8e47b8 100644 --- a/ext/bcmath/libbcmath/src/sub.c +++ b/ext/bcmath/libbcmath/src/sub.c @@ -39,8 +39,7 @@ N2 is subtracted from N1 and the result placed in RESULT. SCALE_MIN is the minimum scale for the result. */ -void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min) -{ +void bc_sub(bc_num n1, bc_num n2, bc_num *result, size_t scale_min) { bc_num diff = NULL; if (n1->n_sign != n2->n_sign) { diff --git a/ext/bcmath/tests/bcpow_zero.phpt b/ext/bcmath/tests/bcpow_zero.phpt index f87a8db6a96ef..6e1e7eb34682a 100644 --- a/ext/bcmath/tests/bcpow_zero.phpt +++ b/ext/bcmath/tests/bcpow_zero.phpt @@ -18,7 +18,9 @@ $baseNumbers = [ "0", ]; -run_bcmath_tests($baseNumbers, $exponents, "**", bcpow(...)); +//run_bcmath_tests($baseNumbers, $exponents, "**", bcpow(...)); + +echo bcpow("0.0000001",1515151551); ?> --EXPECT--