wip

Author: SakiTakamachi

ext/bcmath/libbcmath/src/div.c

@@ -429,3 +429,27 @@ bool bc_divide(bc_num numerator, bc_num divisor, bc_num *quot, size_t scale)
 	*quot = bc_copy_num(BCG(_zero_));
 	return true;
 }
+
+bool bc_divide_vector(
+	BC_VECTOR *numerator_vectors, size_t numerator_arr_size,
+	const BC_VECTOR *divisor_vectors, size_t divisor_arr_size,
+	BC_VECTOR *quot_vectors, size_t quot_arr_size
+) {
+	ZEND_ASSERT(divisor_vectors[divisor_arr_size - 1] != 0);
+	ZEND_ASSERT(quot_arr_size >= numerator_arr_size - divisor_arr_size + 1);
+
+	size_t divisor_size = (divisor_arr_size - 1) * BC_VECTOR_SIZE;
+	BC_VECTOR tmp_divisor_top = divisor_vectors[divisor_arr_size - 1];
+	while (tmp_divisor_top > 0) {
+		divisor_size++;
+		tmp_divisor_top /= BASE;
+	}
+
+	/* 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 {
+		bc_standard_div(numerator_vectors, numerator_arr_size, divisor_vectors, divisor_arr_size, divisor_size, quot_vectors, quot_arr_size);
+	}
+	return true;
+}

ext/bcmath/libbcmath/src/private.h

@@ -87,6 +87,10 @@ bc_num _bc_do_sub (bc_num n1, bc_num n2);
 void bc_multiply_vector(
 	const BC_VECTOR *n1_vector, size_t n1_arr_size, const BC_VECTOR *n2_vector, size_t n2_arr_size,
 	BC_VECTOR *prod_vector, size_t prod_arr_size);
+bool bc_divide_vector(
+	BC_VECTOR *numerator_vectors, size_t numerator_arr_size,
+	const BC_VECTOR *divisor_vectors, size_t divisor_arr_size,
+	BC_VECTOR *quot_vectors, size_t quot_arr_size);
 void _bc_rm_leading_zeros (bc_num num);
 
 #endif

ext/bcmath/libbcmath/src/sqrt.c

@@ -30,6 +30,7 @@
 *************************************************************************/
 
 #include "bcmath.h"
+#include "convert.h"
 #include <stdbool.h>
 #include <stddef.h>
 #include "private.h"
@@ -100,7 +101,7 @@ bool bc_sqrt(bc_num *num, size_t scale)
 	size_t 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) {
+	if (0 && num_calc_full_len < MAX_LENGTH_OF_LONG) {
 		/* Fast path */
 
 		BC_VECTOR n_vector = 0;
@@ -140,53 +141,106 @@ bool bc_sqrt(bc_num *num, size_t scale)
 	} else {
 		/* Standard path */
 
-		bc_num guess;
-		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. */
-			guess = bc_copy_num(BCG(_one_));
-			cscale = (*num)->n_scale;
+		/* convert num to n_vector */
+		size_t n_vector_size = BC_ARR_SIZE_FROM_LEN(num_calc_full_len);
+		BC_VECTOR *buf = safe_emalloc(n_vector_size, sizeof(BC_VECTOR) * 2, 0);
+		BC_VECTOR *n_vector = buf;
+		BC_VECTOR *n_vector_copy = n_vector + n_vector_size;
+
+		size_t n_full_len = (*num)->n_len + (*num)->n_scale;
+		const char *nend = (*num)->n_value + n_full_len - 1;
+		size_t n_extend_zeros = num_calc_full_len - n_full_len;
+
+		bc_convert_to_vector_with_zero_pad(n_vector, nend, n_full_len, n_extend_zeros);
+
+		/* init vectors */
+		size_t guess_len = ((*num)->n_len + 1) / 2;
+		size_t guess_scale = rscale + 1;
+		size_t guess_full_len = guess_len + guess_scale;
+		size_t guess_vector_size = BC_ARR_SIZE_FROM_LEN(guess_full_len);
+
+		BC_VECTOR *buf2 = safe_emalloc(guess_vector_size, sizeof(BC_VECTOR) * 2, 0);
+		BC_VECTOR *guess_vector = buf2;
+		BC_VECTOR *guess1_vector = guess_vector + guess_vector_size;
+
+		size_t guess_quot_size = guess_vector_size + 1;
+		BC_VECTOR *tmp_div_ret_vector = safe_emalloc(guess_quot_size, sizeof(BC_VECTOR), 0);
+
+		for (size_t i = 0; i < guess_vector_size - 1; i++) {
+			guess_vector[i] = BC_VECTOR_BOUNDARY_NUM - 1;
+		}
+		if (guess_full_len % BC_VECTOR_SIZE == 0) {
+			guess_vector[guess_vector_size - 1] = BC_VECTOR_BOUNDARY_NUM / 10;
 		} else {
-			/* The number is greater than 1.  Guess should start at 10^(exp/2). */
-			/* If just divide size_t by 2 it will not overflow. */
-			size_t exponent_for_initial_guess = (size_t) (*num)->n_len >> 1;
-
-			/* 10^n is a 1 followed by n zeros. */
-			guess = bc_new_num(exponent_for_initial_guess + 1, 0);
-			guess->n_value[0] = 1;
-			cscale = 3;
+			guess_vector[guess_vector_size - 1] = 0;
+			for (size_t i = 0; i < guess_full_len % BC_VECTOR_SIZE; i++) {
+				guess_vector[guess_vector_size - 1] *= BASE;
+				guess_vector[guess_vector_size - 1] += 9;
+			}
 		}
 
-		bc_num guess1 = NULL;
-		bc_num point5 = bc_new_num (1, 1);
-		point5->n_value[1] = 5;
-		bc_num diff = NULL;
-
-		bool done = false;
-		while (!done) {
-			bc_free_num (&guess1);
-			guess1 = bc_copy_num(guess);
-			bc_divide(*num, guess, &guess, cscale);
-			bc_add_ex(guess, guess1, &guess, 0);
-			bc_multiply_ex(guess, point5, &guess, cscale);
-			bc_sub_ex(guess, guess1, &diff, cscale + 1);
-			if (bc_is_near_zero(diff, cscale)) {
-				if (cscale < rscale + 1) {
-					cscale = MIN (cscale * 3, rscale + 1);
+		BC_VECTOR two[1] = { 2 };
+
+		/* Newton's algorithm. */
+		//bool done = false;
+		for (size_t i = 0; i < 100000; i++) {
+		//do {
+			for (size_t i = 0; i < n_vector_size; i++) {
+				n_vector_copy[i] = n_vector[i];
+			}
+			bool div_ret = bc_divide_vector(n_vector_copy, n_vector_size, guess_vector, guess_vector_size, tmp_div_ret_vector, guess_quot_size);
+			ZEND_ASSERT(div_ret);
+
+			BC_VECTOR *tmp_vptr = guess1_vector;
+			guess1_vector = guess_vector;
+			guess_vector = tmp_vptr;
+			int carry = 0;
+			for (size_t i = 0; i < guess_vector_size; i++) {
+				guess_vector[i] = guess1_vector[i] + tmp_div_ret_vector[i] + carry;
+				if (guess_vector[i] >= BC_VECTOR_BOUNDARY_NUM) {
+					guess_vector[i] -= BC_VECTOR_BOUNDARY_NUM;
+					carry = 1;
 				} else {
-					done = true;
+					carry = 0;
 				}
 			}
+			ZEND_ASSERT(carry == 0);
+			div_ret = bc_divide_vector(guess_vector, guess_vector_size, two, 1, tmp_div_ret_vector, guess_vector_size);
+			ZEND_ASSERT(div_ret);
+
+			for (size_t i = 0; i < guess_vector_size; i++) {
+				guess_vector[i] = tmp_div_ret_vector[i];
+			}
+
+			//size_t diff = guess_vector[0] > guess1_vector[0] ? guess_vector[0] - guess1_vector[0] : guess1_vector[0] - guess_vector[0];
+			// if (diff <= 1) {
+			// 	bool is_same = true;
+			// 	for (size_t i = 1; i < guess_vector_size; i++) {
+			// 		if (guess_vector[i] != guess1_vector[i]) {
+			// 			is_same = false;
+			// 			break;
+			// 		}
+			// 	}
+			// 	done = is_same;
+			// }
+		//} while (!done);
 		}
 
-		/* Assign the number and clean up. */
-		bc_free_num (num);
-		bc_divide(guess, BCG(_one_), num, rscale);
-		bc_free_num (&guess);
-		bc_free_num (&guess1);
-		bc_free_num (&point5);
-		bc_free_num (&diff);
+		bc_num ret = bc_new_num_nonzeroed(guess_len, guess_scale);
+		char *rptr = ret->n_value;
+		char *rend = rptr + guess_full_len - 1;
+
+		bc_convert_vector_to_char(guess_vector, rptr, rend, guess_vector_size);
+
+		ret->n_scale = rscale;
+		_bc_rm_leading_zeros(ret);
+
+		bc_free_num(num);
+		*num = ret;
+
+		efree(buf);
+		efree(buf2);
+		efree(tmp_div_ret_vector);
 	}
 
 	return true;