@@ -180,10 +180,12 @@ static inline void bc_standard_sqrt(bc_num *num, size_t rscale, size_t num_calc_
180
180
}
181
181
guess_vector [guess_arr_size - 1 ] = 0 ;
182
182
183
- size_t quot_size = n_arr_size - (guess_arr_size - 1 ) + 1 ;
184
-
185
183
BC_VECTOR two [1 ] = { 2 };
186
184
185
+ /* The precision (number of vectors) used for the calculation.
186
+ * Since the initial value uses two vectors, the initial precision is set to 2. */
187
+ size_t guess_precision = 2 ;
188
+
187
189
/**
188
190
* Newton's algorithm. Iterative expression is `x_{n+1} = (x_n + a / x_n) / 2`
189
191
* If break down the calculation into detailed steps, it looks like this:
@@ -194,13 +196,23 @@ static inline void bc_standard_sqrt(bc_num *num, size_t rscale, size_t num_calc_
194
196
*/
195
197
bool done = false;
196
198
do {
199
+ size_t guess_offset = guess_arr_size - 1 - guess_precision ;
200
+ size_t n_offset = guess_offset * 2 ;
201
+
197
202
/* Since the value changes during division by successive approximation, use a copied version of it. */
198
- for (size_t i = 0 ; i < n_arr_size ; i ++ ) {
203
+ for (size_t i = n_offset ; i < n_arr_size ; i ++ ) {
199
204
n_vector_copy [i ] = n_vector [i ];
200
205
}
201
206
207
+ size_t n_precision = n_arr_size - n_offset ;
208
+ size_t quot_size = n_precision - (guess_precision ) + 1 ;
209
+
202
210
/* 1. quot = a / x_n */
203
- bool div_ret = bc_divide_vector (n_vector_copy , n_arr_size , guess_vector , guess_arr_size - 1 , tmp_div_ret_vector , quot_size );
211
+ bool div_ret = bc_divide_vector (
212
+ n_vector_copy + n_offset , n_precision ,
213
+ guess_vector + guess_offset , guess_precision ,
214
+ tmp_div_ret_vector + guess_offset , quot_size
215
+ );
204
216
ZEND_ASSERT (div_ret );
205
217
206
218
BC_VECTOR * tmp_vptr = guess1_vector ;
@@ -209,7 +221,7 @@ static inline void bc_standard_sqrt(bc_num *num, size_t rscale, size_t num_calc_
209
221
210
222
/* 2. add = x_n + quot1 */
211
223
int carry = 0 ;
212
- for (size_t i = 0 ; i < guess_arr_size - 1 ; i ++ ) {
224
+ for (size_t i = guess_offset ; i < guess_arr_size ; i ++ ) {
213
225
guess_vector [i ] = guess1_vector [i ] + tmp_div_ret_vector [i ] + carry ;
214
226
if (guess_vector [i ] >= BC_VECTOR_BOUNDARY_NUM ) {
215
227
guess_vector [i ] -= BC_VECTOR_BOUNDARY_NUM ;
@@ -221,24 +233,37 @@ static inline void bc_standard_sqrt(bc_num *num, size_t rscale, size_t num_calc_
221
233
guess_vector [guess_arr_size - 1 ] = carry ;
222
234
223
235
/* 3. x_{n+1} = add / 2 */
224
- div_ret = bc_divide_vector (guess_vector , guess_arr_size , two , 1 , tmp_div_ret_vector , guess_arr_size );
236
+ div_ret = bc_divide_vector (
237
+ guess_vector + guess_offset , guess_precision + 1 ,
238
+ two , 1 ,
239
+ tmp_div_ret_vector + guess_offset , guess_precision + 1
240
+ );
225
241
ZEND_ASSERT (div_ret );
226
242
227
- for (size_t i = 0 ; i < guess_arr_size ; i ++ ) {
243
+ for (size_t i = guess_offset ; i < guess_arr_size ; i ++ ) {
228
244
guess_vector [i ] = tmp_div_ret_vector [i ];
229
245
}
230
246
231
247
/* 4. repeat until the difference between the `x_n` and `x_{n+1}` is less than or equal to 1. */
232
- size_t diff = guess_vector [0 ] > guess1_vector [0 ] ? guess_vector [0 ] - guess1_vector [0 ] : guess1_vector [0 ] - guess_vector [0 ];
248
+ size_t diff = guess_vector [guess_offset ] > guess1_vector [guess_offset ]
249
+ ? guess_vector [guess_offset ] - guess1_vector [guess_offset ]
250
+ : guess1_vector [guess_offset ] - guess_vector [guess_offset ];
233
251
if (diff <= 1 ) {
234
252
bool is_same = true;
235
- for (size_t i = 1 ; i < guess_arr_size - 1 ; i ++ ) {
253
+ for (size_t i = guess_offset + 1 ; i < guess_arr_size - 1 ; i ++ ) {
236
254
if (guess_vector [i ] != guess1_vector [i ]) {
237
255
is_same = false;
238
256
break ;
239
257
}
240
258
}
241
- done = is_same ;
259
+ if (is_same ) {
260
+ if (guess_precision < guess_arr_size - 1 ) {
261
+ /* If the precision has not yet reached the maximum number of digits, it will be increased. */
262
+ guess_precision = MIN (guess_precision * 3 , guess_arr_size - 1 );
263
+ } else {
264
+ done = is_same ;
265
+ }
266
+ }
242
267
}
243
268
} while (!done );
244
269
0 commit comments