create fast path and move old logic to standard path

Author: SakiTakamachi

ext/bcmath/libbcmath/src/sqrt.c

@@ -32,37 +32,78 @@
 #include "bcmath.h"
 #include <stdbool.h>
 #include <stddef.h>
+#include "private.h"
 
 /* Take the square root NUM and return it in NUM with SCALE digits
    after the decimal place. */
 
-bool bc_sqrt(bc_num *num, size_t scale)
+static inline BC_VECTOR bc_sqrt_get_pow_10(size_t exponent)
 {
-	/* Initial checks. */
-	if (bc_is_neg(*num)) {
-		/* Cannot take the square root of a negative number */
-		return false;
+	BC_VECTOR value = 1;
+	while (exponent >= 8) {
+		value *= BC_POW_10_LUT[8];
+		exponent -= 8;
 	}
-	/* Square root of 0 is 0 */
-	if (bc_is_zero(*num)) {
-		bc_free_num (num);
-		*num = bc_copy_num(BCG(_zero_));
-		return true;
+	value *= BC_POW_10_LUT[exponent];
+	return value;
+}
+
+static BC_VECTOR bc_fast_sqrt_vector(BC_VECTOR n_vector)
+{
+	/* Use sqrt() from <math.h> to determine the initial value. */
+	BC_VECTOR guess_vector = (BC_VECTOR) round(sqrt((double) n_vector));
+
+	/* Newton's algorithm. Iterative expression is `x_{n+1} = (x_n + a / x_n) / 2` */
+	BC_VECTOR guess1_vector;
+	size_t diff;
+	do {
+		guess1_vector = guess_vector;
+		guess_vector = (guess1_vector + n_vector / guess1_vector) / 2; /* Iterative expression */
+		diff = guess1_vector > guess_vector ? guess1_vector - guess_vector : guess_vector - guess1_vector;
+	} while (diff > 1);
+	return guess_vector;
+}
+
+static inline void bc_fast_sqrt(bc_num *num, size_t rscale, size_t num_calc_full_len, size_t leading_zeros)
+{
+	BC_VECTOR n_vector = 0;
+	for (size_t i = leading_zeros; i < (*num)->n_len + (*num)->n_scale; i++) {
+		n_vector = n_vector * BASE + (*num)->n_value[i];
+	}
+	/* When calculating the square root of a number using only integer operations,
+	 * need to adjust the digit scale accordingly.
+	 * Considering that the original number is the square of the result,
+	 * if the desired scale of the result is 5, the input number should be scaled
+	 * by twice that, i.e., scale 10. */
+	n_vector *= bc_sqrt_get_pow_10((rscale + 1) * 2 - (*num)->n_scale);
+
+	/* Get sqrt */
+	BC_VECTOR guess_vector = bc_fast_sqrt_vector(n_vector);
+
+	size_t ret_len;
+	if (leading_zeros > 0) {
+		ret_len = 1;
+	} else {
+		ret_len = ((num_calc_full_len + 1) / 2) - (rscale + 1);
 	}
 
-	bcmath_compare_result num_cmp_one = bc_compare(*num, BCG(_one_), (*num)->n_scale);
-	/* Square root of 1 is 1 */
-	if (num_cmp_one == BCMATH_EQUAL) {
-		bc_free_num (num);
-		*num = bc_copy_num(BCG(_one_));
-		return true;
+	bc_num ret = bc_new_num_nonzeroed(ret_len, rscale);
+	char *rptr = ret->n_value;
+	char *rend = rptr + ret_len + rscale - 1;
+
+	guess_vector /= BASE; /* Since the scale of guess_vector is rscale + 1, reduce the scale by 1. */
+	while (rend >= rptr) {
+		*rend-- = guess_vector % BASE;
+		guess_vector /= BASE;
 	}
+	bc_free_num(num);
+	*num = ret;
+}
 
-	/* Initialize the variables. */
-	size_t cscale;
+static inline void bc_standard_sqrt(bc_num *num, size_t rscale, bcmath_compare_result num_cmp_one)
+{
 	bc_num guess;
-	size_t rscale = MAX(scale, (*num)->n_scale);
-
+	size_t cscale;
 	/* Calculate the initial guess. */
 	if (num_cmp_one == BCMATH_RIGHT_GREATER) {
 		/* The number is between 0 and 1.  Guess should start at 1. */
@@ -86,7 +127,6 @@ bool bc_sqrt(bc_num *num, size_t scale)
 	point5->n_value[1] = 5;
 	bc_num diff = NULL;
 
-	/* Find the square root using Newton's algorithm. */
 	bool done = false;
 	while (!done) {
 		bc_free_num (&guess1);
@@ -111,5 +151,53 @@ bool bc_sqrt(bc_num *num, size_t scale)
 	bc_free_num (&guess1);
 	bc_free_num (&point5);
 	bc_free_num (&diff);
+}
+
+bool bc_sqrt(bc_num *num, size_t scale)
+{
+	/* Initial checks. */
+	if (bc_is_neg(*num)) {
+		/* Cannot take the square root of a negative number */
+		return false;
+	}
+	/* Square root of 0 is 0 */
+	if (bc_is_zero(*num)) {
+		bc_free_num (num);
+		*num = bc_copy_num(BCG(_zero_));
+		return true;
+	}
+
+	bcmath_compare_result num_cmp_one = bc_compare(*num, BCG(_one_), (*num)->n_scale);
+	/* Square root of 1 is 1 */
+	if (num_cmp_one == BCMATH_EQUAL) {
+		bc_free_num (num);
+		*num = bc_copy_num(BCG(_one_));
+		return true;
+	}
+
+	/* Initialize the variables. */
+	size_t rscale = MAX(scale, (*num)->n_scale);
+	size_t leading_zeros = 0;
+	size_t num_calc_full_len;
+	if (num_cmp_one == BCMATH_RIGHT_GREATER) {
+		const char *nptr = (*num)->n_value;
+		/* We know that num is not zero, so there will be no overrun. */
+		while (*nptr == 0) {
+			leading_zeros++;
+			nptr++;
+		}
+		size_t leading_zeros_in_frac = leading_zeros - 1;
+		num_calc_full_len = (rscale + 1 - leading_zeros_in_frac) * 2;
+	} else {
+		num_calc_full_len = (*num)->n_len + (rscale + 1) * 2;
+	}
+
+	/* Find the square root using Newton's algorithm. */
+	if (num_calc_full_len < MAX_LENGTH_OF_LONG) {
+		bc_fast_sqrt(num, rscale, num_calc_full_len, leading_zeros);
+	} else {
+		bc_standard_sqrt(num, rscale, num_cmp_one);
+	}
+
 	return true;
 }