@@ -9,7 +9,7 @@ hp_vec2 initialize_hp_vec2(number x, number y) {
99
1010number null_number () {
1111 number res;
12- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
12+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
1313 res.limb [i] = 0u ;
1414 }
1515 res.sign = 1 ;
@@ -18,7 +18,7 @@ number null_number() {
1818
1919number infinite_number () {
2020 number res;
21- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
21+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
2222 res.limb [i] = (1u << 31u ) - 1u ;
2323 }
2424 res.sign = 1 ;
@@ -52,7 +52,7 @@ uint lo(uint a) {
5252}
5353
5454int compare_abs (number a, number b) {
55- for (int i = LIMB_SIZE - 1 ; i >= 0 ; --i) {
55+ for (int i = NUMBER_OF_LIMBS - 1 ; i >= 0 ; --i) {
5656 if (a.limb [i] > b.limb [i]) return 1 ;
5757 if (a.limb [i] < b.limb [i]) return -1 ;
5858 }
@@ -61,7 +61,7 @@ int compare_abs(number a, number b) {
6161
6262bool is_zero (number a) {
6363 bool notzero = false ;
64- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
64+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
6565 notzero = a.limb [i] != 0u || notzero;
6666 }
6767 return !notzero;
@@ -77,7 +77,7 @@ uvec2 sum_with_carry(uint a, uint b) {
7777number abs_sum (number a, number b) {
7878 number c = null_number ();
7979 uint carry = 0u ;
80- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
80+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
8181 uvec2 res = sum_with_carry (sum_with_carry (a.limb [i], b.limb [i]).x , carry);
8282 uint sum = res.x ;
8383 carry = res.y ;
@@ -89,7 +89,7 @@ number abs_sum(number a, number b) {
8989number abs_hp_sub (number a, number b) {
9090 number c = null_number ();
9191 uint borrow = 0u ;
92- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
92+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
9393 uint sub = a.limb [i] - b.limb [i] - borrow;
9494 borrow = uint ((a.limb [i] < b.limb [i] || (a.limb [i] == b.limb [i] && borrow > 0u )));
9595 c.limb [i] = sub;
@@ -149,14 +149,14 @@ number hp_mult(number a, number b) {
149149 c.limb = infinite_number ().limb ;
150150 return c;
151151 }
152- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
152+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
153153 if (a.limb [i] == 0u ) continue ;
154154 uint carry = 0u ;
155155
156- for (int j = 0 ; j < LIMB_SIZE ; ++j) {
156+ for (int j = 0 ; j < NUMBER_OF_LIMBS ; ++j) {
157157 int target = i + j - FRACTIONAL_SIZE;
158158
159- if (target >= LIMB_SIZE ) break ;
159+ if (target >= NUMBER_OF_LIMBS ) break ;
160160
161161 uvec2 prod = product_with_remainder (a.limb [i], b.limb [j]);
162162 if (target < 0 ) {
@@ -184,9 +184,9 @@ number shift_left(number a, int shift) {
184184 number c = null_number ();
185185 c.sign = a.sign ;
186186
187- if (limb_shift >= LIMB_SIZE ) return null_number ();
187+ if (limb_shift >= NUMBER_OF_LIMBS ) return null_number ();
188188
189- for (int i = LIMB_SIZE - 1 ; i >= limb_shift; --i) {
189+ for (int i = NUMBER_OF_LIMBS - 1 ; i >= limb_shift; --i) {
190190 int target = i;
191191 int source = i - limb_shift;
192192
@@ -211,14 +211,14 @@ number shift_right(number a, int shift) {
211211 number c = null_number ();
212212 c.sign = a.sign ;
213213
214- if (limb_shift >= LIMB_SIZE ) return null_number ();
215- for (int i = 0 ; i < LIMB_SIZE - limb_shift; ++i) {
214+ if (limb_shift >= NUMBER_OF_LIMBS ) return null_number ();
215+ for (int i = 0 ; i < NUMBER_OF_LIMBS - limb_shift; ++i) {
216216 int target = i;
217217 int source = i + limb_shift;
218218
219219 uint val = a.limb [source] >> bit_shift;
220220
221- if (source < LIMB_SIZE - 1 && bit_shift > 0 ) {
221+ if (source < NUMBER_OF_LIMBS - 1 && bit_shift > 0 ) {
222222 val |= a.limb [source + 1 ] << (32 - bit_shift);
223223 }
224224 c.limb [target] = val;
@@ -227,7 +227,7 @@ number shift_right(number a, int shift) {
227227}
228228
229229int find_msb (number a) {
230- for (int i = LIMB_SIZE - 1 ; i >= 0 ; --i) {
230+ for (int i = NUMBER_OF_LIMBS - 1 ; i >= 0 ; --i) {
231231 uint x = a.limb [i];
232232 if (x != 0u ) {
233233 int bit_pos = 0 ;
@@ -243,12 +243,15 @@ int find_msb(number a) {
243243}
244244
245245uint get_half (number a, int index) {
246+ if (index >= NUMBER_OF_LIMBS * 2 || index < 0 ) {
247+ return 0u ;
248+ }
246249 uint l = a.limb [index / 2 ];
247250 return uint (index % 2 == 0 ) * hi (l) + uint (index % 2 == 1 ) * lo (l);
248251}
249252
250253void set_half (number& a, int index, uint val) {
251- if (index >= LIMB_SIZE * 2 || index < 0 ) return ;
254+ if (index >= NUMBER_OF_LIMBS * 2 || index < 0 ) return ;
252255 int limb_idx = index / 2 ;
253256 if (index % 2 == 1 ) {
254257 a.limb [limb_idx] = (a.limb [limb_idx] & 0x0000FFFFu ) | (val << 16 );
@@ -263,7 +266,7 @@ number mult_scalar_16(number a, uint b_16) {
263266 c.sign = a.sign ;
264267 uint carry = 0u ;
265268
266- for (int i = 0 ; i < LIMB_SIZE ; ++i) {
269+ for (int i = 0 ; i < NUMBER_OF_LIMBS ; ++i) {
267270 if (a.limb [i] == 0u && carry == 0u ) continue ;
268271
269272 uvec2 prod = product_with_remainder (a.limb [i], b_16);
@@ -279,12 +282,12 @@ number hp_div(number n, number d) {
279282 number q = null_number ();
280283 q.sign = n.sign * d.sign ;
281284
282- n = shift_left (n, FRACTIONAL_SIZE * 32 ) ;
285+ uint u_halves[NUMBER_OF_LIMBS * 2 + FRACTIONAL_SIZE * 2 + 1 ] ;
283286
284287 if (is_zero (d)) {
285288 q = infinite_number ();
286289 return q;
287- };
290+ }
288291 if (d.is_infinite ) return q;
289292
290293 int msb_d = find_msb (d);
@@ -293,50 +296,76 @@ number hp_div(number n, number d) {
293296 int msb_n = find_msb (n);
294297 if (msb_n == -1 ) return q;
295298
299+ for (int i = 0 ; i < NUMBER_OF_LIMBS*2 + FRACTIONAL_SIZE * 2 + 1 ; ++i) {
300+ u_halves[i] = 0u ;
301+ }
302+
303+ for (int i = 0 ; i < NUMBER_OF_LIMBS * 2 ; ++i) {
304+ u_halves[i + FRACTIONAL_SIZE * 2 ] = get_half (n, i);
305+ }
306+
296307 int len_v = (msb_d >> 4 ) + 1 ;
297- int len_u = (msb_n >> 4 ) + 1 ;
308+ int len_u_unshifted = (( msb_n + FRACTIONAL_SIZE * 32 ) >> 4 ) + 1 ;
298309
299- if (len_u < len_v) return q;
310+ if (len_u_unshifted < len_v) return q;
300311
301312 if (len_v == 1 ) {
302313 uint v0 = get_half (d, 0 );
303314 uint rem = 0u ;
304- for (int i = len_u - 1 ; i >= 0 ; --i) {
305- uint dividend = (rem << 16 ) | get_half (n, i) ;
315+ for (int i = len_u_unshifted - 1 ; i >= 0 ; --i) {
316+ uint dividend = (rem << 16 ) | u_halves[i] ;
306317 uint q_i = dividend / v0;
307318 rem = dividend % v0;
308- set_half (q, i, q_i);
319+
320+ if (i < NUMBER_OF_LIMBS * 2 ) {
321+ set_half (q, i, q_i);
322+ }
309323 }
310324 return q;
311325 }
326+
312327 int shift = 15 - (msb_d % 16 );
313- number u = shift_left (n, shift);
314328 number v = shift_left (d, shift);
315329
330+ if (shift > 0 ) {
331+ uint carry_shift = 0u ;
332+ for (int i = 0 ; i < 25 ; ++i) {
333+ uint val = (u_halves[i] << shift) | carry_shift;
334+ carry_shift = u_halves[i] >> (16 - shift);
335+ u_halves[i] = val & 0xFFFFu ;
336+ }
337+ }
338+
339+ int len_u = ((msb_n + FRACTIONAL_SIZE * 32 + shift) >> 4 ) + 1 ;
316340 int m = len_u - len_v;
317341 int n_len = len_v;
318342
319343 uint v_n1 = get_half (v, n_len - 1 );
320344 uint v_n2 = get_half (v, n_len - 2 );
321345
322346 for (int j = m; j >= 0 ; --j) {
323- uint u_jn = get_half (u, j + n_len) ;
324- uint u_jn1 = get_half (u, j + n_len - 1 ) ;
325- uint u_jn2 = get_half (u, j + n_len - 2 ) ;
347+ uint u_jn = u_halves[ j + n_len] ;
348+ uint u_jn1 = u_halves[ j + n_len - 1 ] ;
349+ uint u_jn2 = (j + n_len >= 2 ) ? u_halves[ j + n_len - 2 ] : 0u ;
326350
327351 uint dividend = (u_jn << 16 ) | u_jn1;
328352 uint q_hat, r_hat;
329353
330354 if (u_jn == v_n1) {
331355 q_hat = 0xFFFFu ;
332356 r_hat = u_jn1 + v_n1;
333- }
334- else {
335- q_hat = dividend / v_n1;
336- r_hat = dividend % v_n1;
357+ } else {
358+ if (v_n1 == 0 ) {
359+ q_hat = (1 << 30 + 1 );
360+ r_hat = 1 ;
361+ }
362+ else {
363+ q_hat = dividend / (v_n1);
364+ r_hat = dividend % v_n1;
365+ }
337366 }
338367
339- while (r_hat < 0x10000u && (q_hat * v_n2) >((r_hat << 16 ) | u_jn2)) {
368+ while (r_hat < 0x10000u && (q_hat * v_n2) > ((r_hat << 16 ) | u_jn2)) {
340369 q_hat--;
341370 r_hat += v_n1;
342371 }
@@ -348,26 +377,31 @@ number hp_div(number n, number d) {
348377 k = p >> 16 ;
349378 uint p_lo = p & 0xFFFFu ;
350379
351- uint u_ji = get_half (u, j + i) ;
380+ uint u_ji = u_halves[ j + i] ;
352381 int diff = int (u_ji) - int (p_lo) - int (borrow);
353382
354- set_half (u, j + i, uint (diff) & 0xFFFFu ) ;
383+ u_halves[ j + i] = uint (diff) & 0xFFFFu ;
355384 borrow = (diff < 0 ) ? 1u : 0u ;
356385 }
357- int final_diff = int (get_half (u, j + n_len)) - int (k) - int (borrow);
358- set_half (u, j + n_len, uint (final_diff) & 0xFFFFu );
386+
387+ int final_diff = int (u_halves[j + n_len]) - int (k) - int (borrow);
388+ u_halves[j + n_len] = uint (final_diff) & 0xFFFFu ;
389+
359390 if (final_diff < 0 ) {
360391 q_hat--;
361- uint carry_hp_add = 0u ;
392+ uint carry_add = 0u ;
362393 for (int i = 0 ; i < n_len; ++i) {
363- uint sum = get_half (u, j + i) + get_half (v, i) + carry_hp_add ;
364- set_half (u, j + i, sum & 0xFFFFu ) ;
365- carry_hp_add = sum >> 16 ;
394+ uint sum = u_halves[ j + i] + get_half (v, i) + carry_add ;
395+ u_halves[ j + i] = sum & 0xFFFFu ;
396+ carry_add = sum >> 16 ;
366397 }
367- uint sum_last = get_half (u, j + n_len) + carry_hp_add;
368- set_half (u, j + n_len, sum_last & 0xFFFFu );
398+ uint sum_last = u_halves[j + n_len] + carry_add;
399+ u_halves[j + n_len] = sum_last & 0xFFFFu ;
400+ }
401+
402+ if (j < NUMBER_OF_LIMBS * 2 ) {
403+ set_half (q, j, q_hat);
369404 }
370- set_half (q, j, q_hat);
371405 }
372406 return q;
373407}
@@ -376,7 +410,7 @@ number div_uint(number n, uint d) {
376410 number q = null_number ();
377411 q.sign = n.sign ;
378412 uint rem = 0u ;
379- for (int i = LIMB_SIZE * 2 - 1 ; i >= 0 ; --i) {
413+ for (int i = NUMBER_OF_LIMBS * 2 - 1 ; i >= 0 ; --i) {
380414 uint dividend = (rem << 16 ) | get_half (n, i);
381415 uint q_i = dividend / d;
382416 rem = dividend % d;
@@ -476,7 +510,7 @@ number hp_sqrt(number x) {
476510
477511 number n_k_next = null_number ();
478512
479- for (int i = 0 ; i < (LIMB_SIZE / 23 + 2 ); ++i) {
513+ for (int i = 0 ; i < (NUMBER_OF_LIMBS / 23 * 32 + 2 ); ++i) {
480514 number div_term = hp_div (x, n_k);
481515 n_k_next = hp_add (n_k, div_term);
482516 n_k_next = shift_right (n_k_next, 1 );
@@ -488,4 +522,3 @@ number hp_sqrt(number x) {
488522 return n_k;
489523}
490524
491- inline number::number () : limb(LIMB_SIZE, 0 ), sign(1 ), is_infinite(false ) {}
0 commit comments