Skip to content

Commit 11ac1ea

Browse files
committed
wip
1 parent 683c823 commit 11ac1ea

File tree

3 files changed

+121
-39
lines changed

3 files changed

+121
-39
lines changed

ext/bcmath/libbcmath/src/div.c

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -429,3 +429,27 @@ bool bc_divide(bc_num numerator, bc_num divisor, bc_num *quot, size_t scale)
429429
*quot = bc_copy_num(BCG(_zero_));
430430
return true;
431431
}
432+
433+
bool bc_divide_vector(
434+
BC_VECTOR *numerator_vectors, size_t numerator_arr_size,
435+
const BC_VECTOR *divisor_vectors, size_t divisor_arr_size,
436+
BC_VECTOR *quot_vectors, size_t quot_arr_size
437+
) {
438+
ZEND_ASSERT(divisor_vectors[divisor_arr_size - 1] != 0);
439+
ZEND_ASSERT(quot_arr_size == numerator_arr_size - divisor_arr_size + 1);
440+
441+
size_t divisor_size = (divisor_arr_size - 1) * BC_VECTOR_SIZE;
442+
BC_VECTOR tmp_divisor_top = divisor_vectors[divisor_arr_size - 1];
443+
while (tmp_divisor_top > 0) {
444+
divisor_size++;
445+
tmp_divisor_top /= BASE;
446+
}
447+
448+
/* Do the division */
449+
if (divisor_arr_size == 1) {
450+
bc_fast_div(numerator_vectors, numerator_arr_size, divisor_vectors[0], quot_vectors, quot_arr_size);
451+
} else {
452+
bc_standard_div(numerator_vectors, numerator_arr_size, divisor_vectors, divisor_arr_size, divisor_size, quot_vectors, quot_arr_size);
453+
}
454+
return true;
455+
}

ext/bcmath/libbcmath/src/private.h

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,10 @@ bc_num _bc_do_sub (bc_num n1, bc_num n2);
8787
void bc_multiply_vector(
8888
const BC_VECTOR *n1_vector, size_t n1_arr_size, const BC_VECTOR *n2_vector, size_t n2_arr_size,
8989
BC_VECTOR *prod_vector, size_t prod_arr_size);
90+
bool bc_divide_vector(
91+
BC_VECTOR *numerator_vectors, size_t numerator_arr_size,
92+
const BC_VECTOR *divisor_vectors, size_t divisor_arr_size,
93+
BC_VECTOR *quot_vectors, size_t quot_arr_size);
9094
void _bc_rm_leading_zeros (bc_num num);
9195

9296
#endif

ext/bcmath/libbcmath/src/sqrt.c

Lines changed: 93 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030
*************************************************************************/
3131

3232
#include "bcmath.h"
33+
#include "convert.h"
3334
#include <stdbool.h>
3435
#include <stddef.h>
3536
#include "private.h"
@@ -100,7 +101,7 @@ bool bc_sqrt(bc_num *num, size_t scale)
100101
size_t num_calc_full_len = (*num)->n_len + (rscale + 1) * 2;
101102

102103
/* Find the square root using Newton's algorithm. */
103-
if (num_calc_full_len < MAX_LENGTH_OF_LONG) {
104+
if (0 && num_calc_full_len < MAX_LENGTH_OF_LONG) {
104105
/* Fast path */
105106

106107
BC_VECTOR n_vector = 0;
@@ -140,53 +141,106 @@ bool bc_sqrt(bc_num *num, size_t scale)
140141
} else {
141142
/* Standard path */
142143

143-
bc_num guess;
144-
size_t cscale;
145-
/* Calculate the initial guess. */
146-
if (num_cmp_one == BCMATH_RIGHT_GREATER) {
147-
/* The number is between 0 and 1. Guess should start at 1. */
148-
guess = bc_copy_num(BCG(_one_));
149-
cscale = (*num)->n_scale;
144+
/* convert num to n_vector */
145+
size_t n_vector_size = BC_ARR_SIZE_FROM_LEN(num_calc_full_len);
146+
BC_VECTOR *buf = safe_emalloc(n_vector_size, sizeof(BC_VECTOR) * 2, 0);
147+
BC_VECTOR *n_vector = buf;
148+
BC_VECTOR *n_vector_copy = n_vector + n_vector_size;
149+
150+
size_t n_full_len = (*num)->n_len + (*num)->n_scale;
151+
const char *nend = (*num)->n_value + n_full_len - 1;
152+
size_t n_extend_zeros = num_calc_full_len - n_full_len;
153+
154+
bc_convert_to_vector_with_zero_pad(n_vector, nend, n_full_len, n_extend_zeros);
155+
156+
/* init vectors */
157+
size_t guess_len = ((*num)->n_len + 1) / 2;
158+
size_t guess_scale = rscale + 1;
159+
size_t guess_full_len = guess_len + guess_scale;
160+
size_t guess_vector_size = BC_ARR_SIZE_FROM_LEN(guess_full_len) + 1;
161+
162+
BC_VECTOR *buf2 = safe_emalloc(guess_vector_size, sizeof(BC_VECTOR) * 3, 0);
163+
BC_VECTOR *guess_vector = buf2;
164+
BC_VECTOR *guess1_vector = guess_vector + guess_vector_size;
165+
BC_VECTOR *tmp_div_ret_vector = guess1_vector + guess_vector_size;
166+
167+
for (size_t i = 0; i < guess_vector_size - 2; i++) {
168+
guess_vector[i] = BC_VECTOR_BOUNDARY_NUM - 1;
169+
}
170+
if (guess_full_len % BC_VECTOR_SIZE == 0) {
171+
guess_vector[guess_vector_size - 2] = BC_VECTOR_BOUNDARY_NUM - 1;
150172
} else {
151-
/* The number is greater than 1. Guess should start at 10^(exp/2). */
152-
/* If just divide size_t by 2 it will not overflow. */
153-
size_t exponent_for_initial_guess = (size_t) (*num)->n_len >> 1;
154-
155-
/* 10^n is a 1 followed by n zeros. */
156-
guess = bc_new_num(exponent_for_initial_guess + 1, 0);
157-
guess->n_value[0] = 1;
158-
cscale = 3;
173+
guess_vector[guess_vector_size - 2] = 0;
174+
for (size_t i = 0; i < guess_full_len % BC_VECTOR_SIZE; i++) {
175+
guess_vector[guess_vector_size - 2] *= BASE;
176+
guess_vector[guess_vector_size - 2] += 9;
177+
}
159178
}
179+
guess_vector[guess_vector_size - 1] = 0;
180+
guess1_vector[guess_vector_size - 1] = 0;
160181

161-
bc_num guess1 = NULL;
162-
bc_num point5 = bc_new_num (1, 1);
163-
point5->n_value[1] = 5;
164-
bc_num diff = NULL;
182+
size_t quot_size = n_vector_size - (guess_vector_size - 1) + 1;
165183

184+
BC_VECTOR two[1] = { 2 };
185+
186+
/* Newton's algorithm. */
166187
bool done = false;
167-
while (!done) {
168-
bc_free_num (&guess1);
169-
guess1 = bc_copy_num(guess);
170-
bc_divide(*num, guess, &guess, cscale);
171-
bc_add_ex(guess, guess1, &guess, 0);
172-
bc_multiply_ex(guess, point5, &guess, cscale);
173-
bc_sub_ex(guess, guess1, &diff, cscale + 1);
174-
if (bc_is_near_zero(diff, cscale)) {
175-
if (cscale < rscale + 1) {
176-
cscale = MIN (cscale * 3, rscale + 1);
188+
do {
189+
for (size_t i = 0; i < n_vector_size; i++) {
190+
n_vector_copy[i] = n_vector[i];
191+
}
192+
bool div_ret = bc_divide_vector(n_vector_copy, n_vector_size, guess_vector, guess_vector_size - 1, tmp_div_ret_vector, quot_size);
193+
ZEND_ASSERT(div_ret);
194+
195+
BC_VECTOR *tmp_vptr = guess1_vector;
196+
guess1_vector = guess_vector;
197+
guess_vector = tmp_vptr;
198+
int carry = 0;
199+
for (size_t i = 0; i < guess_vector_size - 1; i++) {
200+
guess_vector[i] = guess1_vector[i] + tmp_div_ret_vector[i] + carry;
201+
if (guess_vector[i] >= BC_VECTOR_BOUNDARY_NUM) {
202+
guess_vector[i] -= BC_VECTOR_BOUNDARY_NUM;
203+
carry = 1;
177204
} else {
178-
done = true;
205+
carry = 0;
179206
}
180207
}
181-
}
208+
guess_vector[guess_vector_size - 1] = carry;
182209

183-
/* Assign the number and clean up. */
184-
bc_free_num (num);
185-
bc_divide(guess, BCG(_one_), num, rscale);
186-
bc_free_num (&guess);
187-
bc_free_num (&guess1);
188-
bc_free_num (&point5);
189-
bc_free_num (&diff);
210+
div_ret = bc_divide_vector(guess_vector, guess_vector_size, two, 1, tmp_div_ret_vector, guess_vector_size);
211+
ZEND_ASSERT(div_ret);
212+
213+
for (size_t i = 0; i < guess_vector_size; i++) {
214+
guess_vector[i] = tmp_div_ret_vector[i];
215+
}
216+
217+
size_t diff = guess_vector[0] > guess1_vector[0] ? guess_vector[0] - guess1_vector[0] : guess1_vector[0] - guess_vector[0];
218+
if (diff <= 1) {
219+
bool is_same = true;
220+
for (size_t i = 1; i < guess_vector_size - 1; i++) {
221+
if (guess_vector[i] != guess1_vector[i]) {
222+
is_same = false;
223+
break;
224+
}
225+
}
226+
done = is_same;
227+
}
228+
} while (!done);
229+
230+
bc_num ret = bc_new_num_nonzeroed(guess_len, guess_scale);
231+
char *rptr = ret->n_value;
232+
char *rend = rptr + guess_full_len - 1;
233+
234+
bc_convert_vector_to_char(guess_vector, rptr, rend, guess_vector_size - 1);
235+
236+
ret->n_scale = rscale;
237+
_bc_rm_leading_zeros(ret);
238+
239+
bc_free_num(num);
240+
*num = ret;
241+
242+
efree(buf);
243+
efree(buf2);
190244
}
191245

192246
return true;

0 commit comments

Comments
 (0)