|
30 | 30 | *************************************************************************/
|
31 | 31 |
|
32 | 32 | #include "bcmath.h"
|
| 33 | +#include "convert.h" |
33 | 34 | #include <stdbool.h>
|
34 | 35 | #include <stddef.h>
|
35 | 36 | #include "private.h"
|
@@ -109,57 +110,126 @@ static inline void bc_fast_sqrt(bc_num *num, size_t rscale)
|
109 | 110 | *num = ret;
|
110 | 111 | }
|
111 | 112 |
|
112 |
| -static inline void bc_standard_sqrt(bc_num *num, size_t rscale, bcmath_compare_result num_cmp_one) |
| 113 | +static inline void bc_standard_sqrt(bc_num *num, size_t rscale, size_t num_calc_full_len) |
113 | 114 | {
|
114 |
| - bc_num guess; |
115 |
| - size_t cscale; |
116 |
| - /* Calculate the initial guess. */ |
117 |
| - if (num_cmp_one == BCMATH_RIGHT_GREATER) { |
118 |
| - /* The number is between 0 and 1. Guess should start at 1. */ |
119 |
| - guess = bc_copy_num(BCG(_one_)); |
120 |
| - cscale = (*num)->n_scale; |
| 115 | + /* allocate memory */ |
| 116 | + size_t n_arr_size = BC_ARR_SIZE_FROM_LEN(num_calc_full_len); |
| 117 | + |
| 118 | + size_t guess_len = ((*num)->n_len + 1) / 2; |
| 119 | + size_t guess_scale = rscale + 1; |
| 120 | + size_t guess_full_len = guess_len + guess_scale; |
| 121 | + /* Since add the old guess and the new guess together during the calculation, |
| 122 | + * there is a chance of overflow, so allocate an extra size. */ |
| 123 | + size_t guess_arr_size = BC_ARR_SIZE_FROM_LEN(guess_full_len) + 1; |
| 124 | + |
| 125 | + size_t allocate_size = n_arr_size * 2 + guess_arr_size * 3; |
| 126 | + BC_VECTOR *buf = safe_emalloc(allocate_size, sizeof(BC_VECTOR), 0); |
| 127 | + |
| 128 | + BC_VECTOR *n_vector = buf; |
| 129 | + /* In division by successive approximation, the numerator is modified during the computation, |
| 130 | + * so it must be copied each time. */ |
| 131 | + BC_VECTOR *n_vector_copy = n_vector + n_arr_size; |
| 132 | + BC_VECTOR *guess_vector = n_vector_copy + n_arr_size; |
| 133 | + BC_VECTOR *guess1_vector = guess_vector + guess_arr_size; |
| 134 | + BC_VECTOR *tmp_div_ret_vector = guess1_vector + guess_arr_size; |
| 135 | + |
| 136 | + /* convert num to n_vector */ |
| 137 | + size_t n_full_len = (*num)->n_len + (*num)->n_scale; |
| 138 | + const char *nend = (*num)->n_value + n_full_len - 1; |
| 139 | + size_t n_extend_zeros = num_calc_full_len - n_full_len; |
| 140 | + |
| 141 | + bc_convert_to_vector_with_zero_pad(n_vector, nend, n_full_len, n_extend_zeros); |
| 142 | + |
| 143 | + /* Prepare guess_vector (Temporary implementation) */ |
| 144 | + for (size_t i = 0; i < guess_arr_size - 2; i++) { |
| 145 | + guess_vector[i] = BC_VECTOR_BOUNDARY_NUM - 1; |
| 146 | + } |
| 147 | + if (guess_full_len % BC_VECTOR_SIZE == 0) { |
| 148 | + guess_vector[guess_arr_size - 2] = BC_VECTOR_BOUNDARY_NUM - 1; |
121 | 149 | } else {
|
122 |
| - /* The number is greater than 1. Guess should start at 10^(exp/2). */ |
123 |
| - bc_init_num(&guess); |
124 |
| - bc_int2num(&guess, 10); |
125 |
| - |
126 |
| - bc_int2num(&guess1, (*num)->n_len); |
127 |
| - bc_multiply_ex(guess1, point5, &guess1, 0); |
128 |
| - guess1->n_scale = 0; |
129 |
| - bc_raise_bc_exponent(guess, guess1, &guess, 0); |
130 |
| - bc_free_num (&guess1); |
131 |
| - cscale = 3; |
| 150 | + guess_vector[guess_arr_size - 2] = 0; |
| 151 | + for (size_t i = 0; i < guess_full_len % BC_VECTOR_SIZE; i++) { |
| 152 | + guess_vector[guess_arr_size - 2] *= BASE; |
| 153 | + guess_vector[guess_arr_size - 2] += 9; |
| 154 | + } |
132 | 155 | }
|
| 156 | + guess_vector[guess_arr_size - 1] = 0; |
133 | 157 |
|
134 |
| - bc_num guess1 = NULL; |
135 |
| - bc_num point5 = bc_new_num (1, 1); |
136 |
| - point5->n_value[1] = 5; |
137 |
| - bc_num diff = NULL; |
| 158 | + size_t quot_size = n_arr_size - (guess_arr_size - 1) + 1; |
138 | 159 |
|
| 160 | + BC_VECTOR two[1] = { 2 }; |
| 161 | + |
| 162 | + /** |
| 163 | + * Newton's algorithm. Iterative expression is `x_{n+1} = (x_n + a / x_n) / 2` |
| 164 | + * If break down the calculation into detailed steps, it looks like this: |
| 165 | + * 1. quot = a / x_n |
| 166 | + * 2. add = x_n + quot1 |
| 167 | + * 3. x_{n+1} = add / 2 |
| 168 | + * 4. repeat until the difference between the `x_n` and `x_{n+1}` is less than or equal to 1. |
| 169 | + */ |
139 | 170 | bool done = false;
|
140 |
| - while (!done) { |
141 |
| - bc_free_num (&guess1); |
142 |
| - guess1 = bc_copy_num(guess); |
143 |
| - bc_divide(*num, guess, &guess, cscale); |
144 |
| - bc_add_ex(guess, guess1, &guess, 0); |
145 |
| - bc_multiply_ex(guess, point5, &guess, cscale); |
146 |
| - bc_sub_ex(guess, guess1, &diff, cscale + 1); |
147 |
| - if (bc_is_near_zero(diff, cscale)) { |
148 |
| - if (cscale < rscale + 1) { |
149 |
| - cscale = MIN (cscale * 3, rscale + 1); |
| 171 | + do { |
| 172 | + /* Since the value changes during division by successive approximation, use a copied version of it. */ |
| 173 | + for (size_t i = 0; i < n_arr_size; i++) { |
| 174 | + n_vector_copy[i] = n_vector[i]; |
| 175 | + } |
| 176 | + |
| 177 | + /* 1. quot = a / x_n */ |
| 178 | + bool div_ret = bc_divide_vector(n_vector_copy, n_arr_size, guess_vector, guess_arr_size - 1, tmp_div_ret_vector, quot_size); |
| 179 | + ZEND_ASSERT(div_ret); |
| 180 | + |
| 181 | + BC_VECTOR *tmp_vptr = guess1_vector; |
| 182 | + guess1_vector = guess_vector; |
| 183 | + guess_vector = tmp_vptr; |
| 184 | + |
| 185 | + /* 2. add = x_n + quot1 */ |
| 186 | + int carry = 0; |
| 187 | + for (size_t i = 0; i < guess_arr_size - 1; i++) { |
| 188 | + guess_vector[i] = guess1_vector[i] + tmp_div_ret_vector[i] + carry; |
| 189 | + if (guess_vector[i] >= BC_VECTOR_BOUNDARY_NUM) { |
| 190 | + guess_vector[i] -= BC_VECTOR_BOUNDARY_NUM; |
| 191 | + carry = 1; |
150 | 192 | } else {
|
151 |
| - done = true; |
| 193 | + carry = 0; |
152 | 194 | }
|
153 | 195 | }
|
154 |
| - } |
| 196 | + guess_vector[guess_arr_size - 1] = carry; |
| 197 | + |
| 198 | + /* 3. x_{n+1} = add / 2 */ |
| 199 | + div_ret = bc_divide_vector(guess_vector, guess_arr_size, two, 1, tmp_div_ret_vector, guess_arr_size); |
| 200 | + ZEND_ASSERT(div_ret); |
| 201 | + |
| 202 | + for (size_t i = 0; i < guess_arr_size; i++) { |
| 203 | + guess_vector[i] = tmp_div_ret_vector[i]; |
| 204 | + } |
155 | 205 |
|
156 |
| - /* Assign the number and clean up. */ |
157 |
| - bc_free_num (num); |
158 |
| - bc_divide(guess, BCG(_one_), num, rscale); |
159 |
| - bc_free_num (&guess); |
160 |
| - bc_free_num (&guess1); |
161 |
| - bc_free_num (&point5); |
162 |
| - bc_free_num (&diff); |
| 206 | + /* 4. repeat until the difference between the `x_n` and `x_{n+1}` is less than or equal to 1. */ |
| 207 | + size_t diff = guess_vector[0] > guess1_vector[0] ? guess_vector[0] - guess1_vector[0] : guess1_vector[0] - guess_vector[0]; |
| 208 | + if (diff <= 1) { |
| 209 | + bool is_same = true; |
| 210 | + for (size_t i = 1; i < guess_arr_size - 1; i++) { |
| 211 | + if (guess_vector[i] != guess1_vector[i]) { |
| 212 | + is_same = false; |
| 213 | + break; |
| 214 | + } |
| 215 | + } |
| 216 | + done = is_same; |
| 217 | + } |
| 218 | + } while (!done); |
| 219 | + |
| 220 | + bc_num ret = bc_new_num_nonzeroed(guess_len, guess_scale); |
| 221 | + char *rptr = ret->n_value; |
| 222 | + char *rend = rptr + guess_full_len - 1; |
| 223 | + |
| 224 | + bc_convert_vector_to_char(guess_vector, rptr, rend, guess_arr_size - 1); |
| 225 | + |
| 226 | + ret->n_scale = rscale; |
| 227 | + _bc_rm_leading_zeros(ret); |
| 228 | + |
| 229 | + bc_free_num(num); |
| 230 | + *num = ret; |
| 231 | + |
| 232 | + efree(buf); |
163 | 233 | }
|
164 | 234 |
|
165 | 235 | bool bc_sqrt(bc_num *num, size_t scale)
|
@@ -192,8 +262,7 @@ bool bc_sqrt(bc_num *num, size_t scale)
|
192 | 262 | if (num_calc_full_len < MAX_LENGTH_OF_LONG) {
|
193 | 263 | bc_fast_sqrt(num, rscale);
|
194 | 264 | } else {
|
195 |
| - bc_standard_sqrt(num, rscale, num_cmp_one); |
| 265 | + bc_standard_sqrt(num, rscale, num_calc_full_len); |
196 | 266 | }
|
197 |
| - |
198 | 267 | return true;
|
199 | 268 | }
|
0 commit comments