Line data Source code
1 : /**
2 : * The MIT License (MIT)
3 : *
4 : * Copyright (c) 2021 RSK Labs Ltd
5 : *
6 : * Permission is hereby granted, free of charge, to any person obtaining a copy
7 : * of this software and associated documentation files (the "Software"), to
8 : * deal in the Software without restriction, including without limitation the
9 : * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
10 : * sell copies of the Software, and to permit persons to whom the Software is
11 : * furnished to do so, subject to the following conditions:
12 : *
13 : * The above copyright notice and this permission notice shall be included in
14 : * all copies or substantial portions of the Software.
15 : *
16 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19 : * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
21 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
22 : * IN THE SOFTWARE.
23 : */
24 :
25 : /* $Id: bigdigits.c $ */
26 :
27 : /***** BEGIN LICENSE BLOCK *****
28 : *
29 : * This Source Code Form is subject to the terms of the Mozilla Public
30 : * License, v. 2.0. If a copy of the MPL was not distributed with this
31 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
32 : *
33 : * Copyright (c) 2001-16 David Ireland, D.I. Management Services Pty Limited
34 : * <http://www.di-mgt.com.au/bigdigits.html>. All rights reserved.
35 : *
36 : ***** END LICENSE BLOCK *****/
37 : /*
38 : * Last updated:
39 : * $Date: 2016-03-31 09:51:00 $
40 : * $Revision: 2.6.1 $
41 : * $Author: dai $
42 : */
43 :
44 : /* Core code for BigDigits library "mp" functions */
45 :
46 : /* No asserts in stripped down version of this library */
47 : #define assert(x)
48 :
49 : #include "bigdigits.h"
50 :
51 : #define BITS_PER_HALF_DIGIT (BITS_PER_DIGIT / 2)
52 : #define LOHALF(x) ((DIGIT_T)((x) & MAX_HALF_DIGIT))
53 : #define HIHALF(x) ((DIGIT_T)((x) >> BITS_PER_HALF_DIGIT & MAX_HALF_DIGIT))
54 : #define TOHIGH(x) ((DIGIT_T)((x) << BITS_PER_HALF_DIGIT))
55 :
56 1491 : DIGIT_T mpSetZero(volatile DIGIT_T a[], size_t ndigits)
57 : { /* Sets a = 0 */
58 :
59 : /* Prevent optimiser ignoring this */
60 : volatile DIGIT_T optdummy;
61 1491 : volatile DIGIT_T *p = a;
62 :
63 14914 : while (ndigits--)
64 13423 : a[ndigits] = 0;
65 :
66 1491 : optdummy = *p;
67 1491 : return optdummy;
68 : }
69 :
70 833 : void mpSetDigit(DIGIT_T a[], DIGIT_T d, size_t ndigits)
71 : { /* Sets a = d where d is a single digit */
72 : size_t i;
73 :
74 7497 : for (i = 1; i < ndigits; i++)
75 : {
76 6664 : a[i] = 0;
77 : }
78 833 : a[0] = d;
79 833 : }
80 :
81 : /** Returns number of significant digits in a */
82 710 : size_t mpSizeof(const DIGIT_T a[], size_t ndigits)
83 : {
84 3168 : while(ndigits--)
85 : {
86 3168 : if (a[ndigits] != 0)
87 710 : return (++ndigits);
88 : }
89 0 : return 0;
90 : }
91 :
92 38 : void mpSetEqual(DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
93 : { /* Sets a = b */
94 : size_t i;
95 :
96 380 : for (i = 0; i < ndigits; i++)
97 : {
98 342 : a[i] = b[i];
99 : }
100 38 : }
101 :
102 228 : DIGIT_T mpAdd(DIGIT_T w[], const DIGIT_T u[], const DIGIT_T v[], size_t ndigits)
103 : {
104 : /* Calculates w = u + v
105 : where w, u, v are multiprecision integers of ndigits each
106 : Returns carry if overflow. Carry = 0 or 1.
107 :
108 : Ref: Knuth Vol 2 Ch 4.3.1 p 266 Algorithm A.
109 : */
110 :
111 : DIGIT_T k;
112 : size_t j;
113 :
114 : /* w can't be the same as v
115 : Stop if assert is working, else return error (overflow = -1)
116 : */
117 228 : if (w == v) {
118 : assert(w != v);
119 0 : return MAX_DIGIT;
120 : }
121 :
122 : /* Step A1. Initialise */
123 228 : k = 0;
124 :
125 2280 : for (j = 0; j < ndigits; j++)
126 : {
127 : /* Step A2. Add digits w_j = (u_j + v_j + k)
128 : Set k = 1 if carry (overflow) occurs
129 : */
130 2052 : w[j] = u[j] + k;
131 2052 : if (w[j] < k)
132 0 : k = 1;
133 : else
134 2052 : k = 0;
135 :
136 2052 : w[j] += v[j];
137 2052 : if (w[j] < v[j])
138 238 : k++;
139 :
140 : } /* Step A3. Loop on j */
141 :
142 228 : return k; /* w_n = k */
143 : }
144 :
145 22 : static int QhatTooBig(DIGIT_T qhat, DIGIT_T rhat,
146 : DIGIT_T vn2, DIGIT_T ujn2)
147 : { /* Returns true if Qhat is too big
148 : i.e. if (Qhat * Vn-2) > (b.Rhat + Uj+n-2)
149 : */
150 : DIGIT_T t[2];
151 :
152 22 : spMultiply(t, qhat, vn2);
153 22 : if (t[1] < rhat)
154 22 : return 0;
155 0 : else if (t[1] > rhat)
156 0 : return 1;
157 0 : else if (t[0] > ujn2)
158 0 : return 1;
159 :
160 0 : return 0;
161 : }
162 :
163 44 : static DIGIT_T mpMultSub(DIGIT_T wn, DIGIT_T w[], const DIGIT_T v[],
164 : DIGIT_T q, size_t n)
165 : { /* Compute w = w - qv
166 : where w = (WnW[n-1]...W[0])
167 : return modified Wn.
168 : */
169 : DIGIT_T k, t[2];
170 : size_t i;
171 :
172 44 : if (q == 0) /* No change */
173 22 : return wn;
174 :
175 22 : k = 0;
176 :
177 198 : for (i = 0; i < n; i++)
178 : {
179 176 : spMultiply(t, q, v[i]);
180 176 : w[i] -= k;
181 176 : if (w[i] > MAX_DIGIT - k)
182 0 : k = 1;
183 : else
184 176 : k = 0;
185 176 : w[i] -= t[0];
186 176 : if (w[i] > MAX_DIGIT - t[0])
187 149 : k++;
188 176 : k += t[1];
189 : }
190 :
191 : /* Cope with Wn not stored in array w[0..n-1] */
192 22 : wn -= k;
193 :
194 22 : return wn;
195 : }
196 :
197 355 : int mpDivide(DIGIT_T q[], DIGIT_T r[], const DIGIT_T u[],
198 : size_t udigits, DIGIT_T v[], size_t vdigits)
199 : { /* Computes quotient q = u / v and remainder r = u mod v
200 : where q, r, u are multiple precision digits
201 : all of udigits and the divisor v is vdigits.
202 :
203 : Ref: Knuth Vol 2 Ch 4.3.1 p 272 Algorithm D.
204 :
205 : Do without extra storage space, i.e. use r[] for
206 : normalised u[], unnormalise v[] at end, and cope with
207 : extra digit Uj+n added to u after normalisation.
208 :
209 : WARNING: this trashes q and r first, so cannot do
210 : u = u / v or v = u mod v.
211 : It also changes v temporarily so cannot make it const.
212 : */
213 : size_t shift;
214 : int n, m, j;
215 : DIGIT_T bitmask, overflow;
216 : DIGIT_T qhat, rhat, t[2];
217 : DIGIT_T *uu, *ww;
218 : int qhatOK, cmp;
219 :
220 : /* Clear q and r */
221 355 : mpSetZero(q, udigits);
222 355 : mpSetZero(r, udigits);
223 :
224 : /* Work out exact sizes of u and v */
225 355 : n = (int)mpSizeof(v, vdigits);
226 355 : m = (int)mpSizeof(u, udigits);
227 355 : m -= n;
228 :
229 : /* Catch special cases */
230 355 : if (n == 0)
231 0 : return -1; /* Error: divide by zero */
232 :
233 355 : if (n == 1)
234 : { /* Use short division instead */
235 295 : r[0] = mpShortDiv(q, u, v[0], udigits);
236 295 : return 0;
237 : }
238 :
239 60 : if (m < 0)
240 : { /* v > u, so just set q = 0 and r = u */
241 0 : mpSetEqual(r, u, udigits);
242 0 : return 0;
243 : }
244 :
245 60 : if (m == 0)
246 : { /* u and v are the same length */
247 38 : cmp = mpCompare(u, v, (size_t)n);
248 38 : if (cmp < 0)
249 : { /* v > u, as above */
250 38 : mpSetEqual(r, u, udigits);
251 38 : return 0;
252 : }
253 0 : else if (cmp == 0)
254 : { /* v == u, so set q = 1 and r = 0 */
255 0 : mpSetDigit(q, 1, udigits);
256 0 : return 0;
257 : }
258 : }
259 :
260 : /* In Knuth notation, we have:
261 : Given
262 : u = (Um+n-1 ... U1U0)
263 : v = (Vn-1 ... V1V0)
264 : Compute
265 : q = u/v = (QmQm-1 ... Q0)
266 : r = u mod v = (Rn-1 ... R1R0)
267 : */
268 :
269 : /* Step D1. Normalise */
270 : /* Requires high bit of Vn-1
271 : to be set, so find most signif. bit then shift left,
272 : i.e. d = 2^shift, u' = u * d, v' = v * d.
273 : */
274 22 : bitmask = HIBITMASK;
275 22 : for (shift = 0; shift < BITS_PER_DIGIT; shift++)
276 : {
277 22 : if (v[n-1] & bitmask)
278 22 : break;
279 0 : bitmask >>= 1;
280 : }
281 :
282 : /* Normalise v in situ - NB only shift non-zero digits */
283 22 : overflow = mpShiftLeft(v, v, shift, n);
284 :
285 : /* Copy normalised dividend u*d into r */
286 22 : overflow = mpShiftLeft(r, u, shift, n + m);
287 22 : uu = r; /* Use ptr to keep notation constant */
288 :
289 22 : t[0] = overflow; /* Extra digit Um+n */
290 :
291 : /* Step D2. Initialise j. Set j = m */
292 66 : for (j = m; j >= 0; j--)
293 : {
294 : /* Step D3. Set Qhat = [(b.Uj+n + Uj+n-1)/Vn-1]
295 : and Rhat = remainder */
296 44 : qhatOK = 0;
297 44 : t[1] = t[0]; /* This is Uj+n */
298 44 : t[0] = uu[j+n-1];
299 44 : overflow = spDivide(&qhat, &rhat, t, v[n-1]);
300 :
301 : /* Test Qhat */
302 44 : if (overflow)
303 : { /* Qhat == b so set Qhat = b - 1 */
304 0 : qhat = MAX_DIGIT;
305 0 : rhat = uu[j+n-1];
306 0 : rhat += v[n-1];
307 0 : if (rhat < v[n-1]) /* Rhat >= b, so no re-test */
308 0 : qhatOK = 1;
309 : }
310 : /* [VERSION 2: Added extra test "qhat && "] */
311 44 : if (qhat && !qhatOK && QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
312 : { /* If Qhat.Vn-2 > b.Rhat + Uj+n-2
313 : decrease Qhat by one, increase Rhat by Vn-1
314 : */
315 0 : qhat--;
316 0 : rhat += v[n-1];
317 : /* Repeat this test if Rhat < b */
318 0 : if (!(rhat < v[n-1]))
319 0 : if (QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
320 0 : qhat--;
321 : }
322 :
323 :
324 : /* Step D4. Multiply and subtract */
325 44 : ww = &uu[j];
326 44 : overflow = mpMultSub(t[1], ww, v, qhat, (size_t)n);
327 :
328 : /* Step D5. Test remainder. Set Qj = Qhat */
329 44 : q[j] = qhat;
330 44 : if (overflow)
331 : { /* Step D6. Add back if D4 was negative */
332 0 : q[j]--;
333 0 : overflow = mpAdd(ww, ww, v, (size_t)n);
334 : }
335 :
336 44 : t[0] = uu[j+n-1]; /* Uj+n on next round */
337 :
338 : } /* Step D7. Loop on j */
339 :
340 : /* Clear high digits in uu */
341 44 : for (j = n; j < m+n; j++)
342 22 : uu[j] = 0;
343 :
344 : /* Step D8. Unnormalise. */
345 :
346 22 : mpShiftRight(r, r, shift, n);
347 22 : mpShiftRight(v, v, shift, n);
348 :
349 22 : return 0;
350 : }
351 :
352 295 : DIGIT_T mpShortDiv(DIGIT_T q[], const DIGIT_T u[], DIGIT_T v,
353 : size_t ndigits)
354 : {
355 : /* Calculates quotient q = u div v
356 : Returns remainder r = u mod v
357 : where q, u are multiprecision integers of ndigits each
358 : and r, v are single precision digits.
359 :
360 : Makes no assumptions about normalisation.
361 :
362 : Ref: Knuth Vol 2 Ch 4.3.1 Exercise 16 p625
363 : */
364 : size_t j;
365 : DIGIT_T t[2], r;
366 : size_t shift;
367 : DIGIT_T bitmask, overflow, *uu;
368 :
369 295 : if (ndigits == 0) return 0;
370 295 : if (v == 0) return 0; /* Divide by zero error */
371 :
372 : /* Normalise first */
373 : /* Requires high bit of V
374 : to be set, so find most signif. bit then shift left,
375 : i.e. d = 2^shift, u' = u * d, v' = v * d.
376 : */
377 295 : bitmask = HIBITMASK;
378 9017 : for (shift = 0; shift < BITS_PER_DIGIT; shift++)
379 : {
380 9017 : if (v & bitmask)
381 295 : break;
382 8722 : bitmask >>= 1;
383 : }
384 :
385 295 : v <<= shift;
386 295 : overflow = mpShiftLeft(q, u, shift, ndigits);
387 295 : uu = q;
388 :
389 : /* Step S1 - modified for extra digit. */
390 295 : r = overflow; /* New digit Un */
391 295 : j = ndigits;
392 2950 : while (j--)
393 : {
394 : /* Step S2. */
395 2655 : t[1] = r;
396 2655 : t[0] = uu[j];
397 2655 : overflow = spDivide(&q[j], &r, t, v);
398 : }
399 :
400 : /* Unnormalise */
401 295 : r >>= shift;
402 :
403 295 : return r;
404 : }
405 :
406 339 : DIGIT_T mpShiftLeft(DIGIT_T a[], const DIGIT_T *b,
407 : size_t shift, size_t ndigits)
408 : { /* Computes a = b << shift */
409 : /* [v2.1] Modified to cope with shift > BITS_PERDIGIT */
410 : size_t i, y, nw, bits;
411 : DIGIT_T mask, carry, nextcarry;
412 :
413 : /* Do we shift whole digits? */
414 339 : if (shift >= BITS_PER_DIGIT)
415 : {
416 0 : nw = shift / BITS_PER_DIGIT;
417 0 : i = ndigits;
418 0 : while (i--)
419 : {
420 0 : if (i >= nw)
421 0 : a[i] = b[i-nw];
422 : else
423 0 : a[i] = 0;
424 : }
425 : /* Call again to shift bits inside digits */
426 0 : bits = shift % BITS_PER_DIGIT;
427 0 : carry = b[ndigits-nw] << bits;
428 0 : if (bits)
429 0 : carry |= mpShiftLeft(a, a, bits, ndigits);
430 0 : return carry;
431 : }
432 : else
433 : {
434 339 : bits = shift;
435 : }
436 :
437 : /* Construct mask = high bits set */
438 339 : mask = ~(~(DIGIT_T)0 >> bits);
439 :
440 339 : y = BITS_PER_DIGIT - bits;
441 339 : carry = 0;
442 3368 : for (i = 0; i < ndigits; i++)
443 : {
444 3029 : nextcarry = (b[i] & mask) >> y;
445 3029 : a[i] = b[i] << bits | carry;
446 3029 : carry = nextcarry;
447 : }
448 :
449 339 : return carry;
450 : }
451 :
452 44 : DIGIT_T mpShiftRight(DIGIT_T a[], const DIGIT_T b[], size_t shift, size_t ndigits)
453 : { /* Computes a = b >> shift */
454 : /* [v2.1] Modified to cope with shift > BITS_PERDIGIT */
455 : size_t i, y, nw, bits;
456 : DIGIT_T mask, carry, nextcarry;
457 :
458 : /* Do we shift whole digits? */
459 44 : if (shift >= BITS_PER_DIGIT)
460 : {
461 0 : nw = shift / BITS_PER_DIGIT;
462 0 : for (i = 0; i < ndigits; i++)
463 : {
464 0 : if ((i+nw) < ndigits)
465 0 : a[i] = b[i+nw];
466 : else
467 0 : a[i] = 0;
468 : }
469 : /* Call again to shift bits inside digits */
470 0 : bits = shift % BITS_PER_DIGIT;
471 0 : carry = b[nw-1] >> bits;
472 0 : if (bits)
473 0 : carry |= mpShiftRight(a, a, bits, ndigits);
474 0 : return carry;
475 : }
476 : else
477 : {
478 44 : bits = shift;
479 : }
480 :
481 : /* Construct mask to set low bits */
482 : /* (thanks to Jesse Chisholm for suggesting this improved technique) */
483 44 : mask = ~(~(DIGIT_T)0 << bits);
484 :
485 44 : y = BITS_PER_DIGIT - bits;
486 44 : carry = 0;
487 44 : i = ndigits;
488 396 : while (i--)
489 : {
490 352 : nextcarry = (b[i] & mask) << y;
491 352 : a[i] = b[i] >> bits | carry;
492 352 : carry = nextcarry;
493 : }
494 :
495 44 : return carry;
496 : }
497 :
498 : /** Returns sign of (a - b) as 0, +1 or -1 (constant-time) */
499 1094 : int mpCompare_ct(const DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
500 : {
501 : /* All these vars are either 0 or 1 */
502 1094 : unsigned int gt = 0;
503 1094 : unsigned int lt = 0;
504 1094 : unsigned int mask = 1; /* Set to zero once first inequality found */
505 : unsigned int c;
506 :
507 10940 : while (ndigits--) {
508 9846 : gt |= (a[ndigits] > b[ndigits]) & mask;
509 9846 : lt |= (a[ndigits] < b[ndigits]) & mask;
510 9846 : c = (gt | lt);
511 9846 : mask &= (c-1); /* Unchanged if c==0 or mask==0, else mask=0 */
512 : }
513 :
514 1094 : return (int)gt - (int)lt; /* EQ=0 GT=+1 LT=-1 */
515 : }
516 :
517 : /** Returns sign of (a - b) as 0, +1 or -1. Not constant-time. */
518 38 : int mpCompare(const DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
519 : {
520 : /* if (ndigits == 0) return 0; // deleted [v2.5] */
521 :
522 38 : while (ndigits--)
523 : {
524 38 : if (a[ndigits] > b[ndigits])
525 0 : return 1; /* GT */
526 38 : if (a[ndigits] < b[ndigits])
527 38 : return -1; /* LT */
528 : }
529 :
530 0 : return 0; /* EQ */
531 : }
532 :
533 198 : int spMultiply(DIGIT_T p[2], DIGIT_T x, DIGIT_T y)
534 : { /* Computes p = x * y */
535 : /* Ref: Arbitrary Precision Computation
536 : http://numbers.computation.free.fr/Constants/constants.html
537 :
538 : high p1 p0 low
539 : +--------+--------+--------+--------+
540 : | x1*y1 | x0*y0 |
541 : +--------+--------+--------+--------+
542 : +-+--------+--------+
543 : |1| (x0*y1 + x1*y1) |
544 : +-+--------+--------+
545 : ^carry from adding (x0*y1+x1*y1) together
546 : +-+
547 : |1|< carry from adding LOHALF t
548 : +-+ to high half of p0
549 : */
550 : DIGIT_T x0, y0, x1, y1;
551 : DIGIT_T t, u, carry;
552 :
553 : /* Split each x,y into two halves
554 : x = x0 + B*x1
555 : y = y0 + B*y1
556 : where B = 2^16, half the digit size
557 : Product is
558 : xy = x0y0 + B(x0y1 + x1y0) + B^2(x1y1)
559 : */
560 :
561 198 : x0 = LOHALF(x);
562 198 : x1 = HIHALF(x);
563 198 : y0 = LOHALF(y);
564 198 : y1 = HIHALF(y);
565 :
566 : /* Calc low part - no carry */
567 198 : p[0] = x0 * y0;
568 :
569 : /* Calc middle part */
570 198 : t = x0 * y1;
571 198 : u = x1 * y0;
572 198 : t += u;
573 198 : if (t < u)
574 0 : carry = 1;
575 : else
576 198 : carry = 0;
577 :
578 : /* This carry will go to high half of p[1]
579 : + high half of t into low half of p[1] */
580 198 : carry = TOHIGH(carry) + HIHALF(t);
581 :
582 : /* Add low half of t to high half of p[0] */
583 198 : t = TOHIGH(t);
584 198 : p[0] += t;
585 198 : if (p[0] < t)
586 0 : carry++;
587 :
588 198 : p[1] = x1 * y1;
589 198 : p[1] += carry;
590 :
591 :
592 198 : return 0;
593 : }
594 :
595 : /* spDivide */
596 :
597 : #define B (MAX_HALF_DIGIT + 1)
598 :
599 5398 : static void spMultSub(DIGIT_T uu[2], DIGIT_T qhat, DIGIT_T v1, DIGIT_T v0)
600 : {
601 : /* Compute uu = uu - q(v1v0)
602 : where uu = u3u2u1u0, u3 = 0
603 : and u_n, v_n are all half-digits
604 : even though v1, v2 are passed as full digits.
605 : */
606 : DIGIT_T p0, p1, t;
607 :
608 5398 : p0 = qhat * v0;
609 5398 : p1 = qhat * v1;
610 5398 : t = p0 + TOHIGH(LOHALF(p1));
611 5398 : uu[0] -= t;
612 5398 : if (uu[0] > MAX_DIGIT - t)
613 4646 : uu[1]--; /* Borrow */
614 5398 : uu[1] -= HIHALF(p1);
615 5398 : }
616 :
617 2699 : DIGIT_T spDivide(DIGIT_T *q, DIGIT_T *r, const DIGIT_T u[2], DIGIT_T v)
618 : { /* Computes quotient q = u / v, remainder r = u mod v
619 : where u is a double digit
620 : and q, v, r are single precision digits.
621 : Returns high digit of quotient (max value is 1)
622 : CAUTION: Assumes normalised such that v1 >= b/2
623 : where b is size of HALF_DIGIT
624 : i.e. the most significant bit of v should be one
625 :
626 : In terms of half-digits in Knuth notation:
627 : (q2q1q0) = (u4u3u2u1u0) / (v1v0)
628 : (r1r0) = (u4u3u2u1u0) mod (v1v0)
629 : for m = 2, n = 2 where u4 = 0
630 : q2 is either 0 or 1.
631 : We set q = (q1q0) and return q2 as "overflow"
632 : */
633 : DIGIT_T qhat, rhat, t, v0, v1, u0, u1, u2, u3;
634 : DIGIT_T uu[2], q2;
635 :
636 : /* Check for normalisation */
637 2699 : if (!(v & HIBITMASK))
638 : { /* Stop if assert is working, else return error */
639 : assert(v & HIBITMASK);
640 0 : *q = *r = 0;
641 0 : return MAX_DIGIT;
642 : }
643 :
644 : /* Split up into half-digits */
645 2699 : v0 = LOHALF(v);
646 2699 : v1 = HIHALF(v);
647 2699 : u0 = LOHALF(u[0]);
648 2699 : u1 = HIHALF(u[0]);
649 2699 : u2 = LOHALF(u[1]);
650 2699 : u3 = HIHALF(u[1]);
651 :
652 : /* Do three rounds of Knuth Algorithm D Vol 2 p272 */
653 :
654 : /* ROUND 1. Set j = 2 and calculate q2 */
655 : /* Estimate qhat = (u4u3)/v1 = 0 or 1
656 : then set (u4u3u2) -= qhat(v1v0)
657 : where u4 = 0.
658 : */
659 : /* [Replaced in Version 2] -->
660 : qhat = u3 / v1;
661 : if (qhat > 0)
662 : {
663 : rhat = u3 - qhat * v1;
664 : t = TOHIGH(rhat) | u2;
665 : if (qhat * v0 > t)
666 : qhat--;
667 : }
668 : <-- */
669 2699 : qhat = (u3 < v1 ? 0 : 1);
670 2699 : if (qhat > 0)
671 : { /* qhat is one, so no need to mult */
672 0 : rhat = u3 - v1;
673 : /* t = r.b + u2 */
674 0 : t = TOHIGH(rhat) | u2;
675 0 : if (v0 > t)
676 0 : qhat--;
677 : }
678 :
679 2699 : uu[1] = 0; /* (u4) */
680 2699 : uu[0] = u[1]; /* (u3u2) */
681 2699 : if (qhat > 0)
682 : {
683 : /* (u4u3u2) -= qhat(v1v0) where u4 = 0 */
684 0 : spMultSub(uu, qhat, v1, v0);
685 0 : if (HIHALF(uu[1]) != 0)
686 : { /* Add back */
687 0 : qhat--;
688 0 : uu[0] += v;
689 0 : uu[1] = 0;
690 : }
691 : }
692 2699 : q2 = qhat;
693 :
694 : /* ROUND 2. Set j = 1 and calculate q1 */
695 : /* Estimate qhat = (u3u2) / v1
696 : then set (u3u2u1) -= qhat(v1v0)
697 : */
698 2699 : t = uu[0];
699 2699 : qhat = t / v1;
700 2699 : rhat = t - qhat * v1;
701 : /* Test on v0 */
702 2699 : t = TOHIGH(rhat) | u1;
703 2699 : if ((qhat == B) || (qhat * v0 > t))
704 : {
705 0 : qhat--;
706 0 : rhat += v1;
707 0 : t = TOHIGH(rhat) | u1;
708 0 : if ((rhat < B) && (qhat * v0 > t))
709 0 : qhat--;
710 : }
711 :
712 : /* Multiply and subtract
713 : (u3u2u1)' = (u3u2u1) - qhat(v1v0)
714 : */
715 2699 : uu[1] = HIHALF(uu[0]); /* (0u3) */
716 2699 : uu[0] = TOHIGH(LOHALF(uu[0])) | u1; /* (u2u1) */
717 2699 : spMultSub(uu, qhat, v1, v0);
718 2699 : if (HIHALF(uu[1]) != 0)
719 : { /* Add back */
720 0 : qhat--;
721 0 : uu[0] += v;
722 0 : uu[1] = 0;
723 : }
724 :
725 : /* q1 = qhat */
726 2699 : *q = TOHIGH(qhat);
727 :
728 : /* ROUND 3. Set j = 0 and calculate q0 */
729 : /* Estimate qhat = (u2u1) / v1
730 : then set (u2u1u0) -= qhat(v1v0)
731 : */
732 2699 : t = uu[0];
733 2699 : qhat = t / v1;
734 2699 : rhat = t - qhat * v1;
735 : /* Test on v0 */
736 2699 : t = TOHIGH(rhat) | u0;
737 2699 : if ((qhat == B) || (qhat * v0 > t))
738 : {
739 0 : qhat--;
740 0 : rhat += v1;
741 0 : t = TOHIGH(rhat) | u0;
742 0 : if ((rhat < B) && (qhat * v0 > t))
743 0 : qhat--;
744 : }
745 :
746 : /* Multiply and subtract
747 : (u2u1u0)" = (u2u1u0)' - qhat(v1v0)
748 : */
749 2699 : uu[1] = HIHALF(uu[0]); /* (0u2) */
750 2699 : uu[0] = TOHIGH(LOHALF(uu[0])) | u0; /* (u1u0) */
751 2699 : spMultSub(uu, qhat, v1, v0);
752 2699 : if (HIHALF(uu[1]) != 0)
753 : { /* Add back */
754 0 : qhat--;
755 0 : uu[0] += v;
756 0 : uu[1] = 0;
757 : }
758 :
759 : /* q0 = qhat */
760 2699 : *q |= LOHALF(qhat);
761 :
762 : /* Remainder is in (u1u0) i.e. uu[0] */
763 2699 : *r = uu[0];
764 2699 : return q2;
765 : }
766 :
767 : // Platform-dependent code
768 : #ifndef HSM_PLATFORM_LEDGER
769 :
770 : #include "hal/log.h"
771 :
772 2022 : void LOG_BIGD_HEX(const char *prefix,
773 : const DIGIT_T *a,
774 : size_t len,
775 : const char *suffix) {
776 2022 : if (prefix) {
777 2022 : LOG("%s", prefix);
778 : }
779 : /* Trim leading digits which are zero */
780 11728 : while (len--) {
781 11708 : if (a[len] != 0)
782 2002 : break;
783 : }
784 2022 : len++;
785 2022 : if (0 == len)
786 20 : len = 1;
787 : /* print first digit without leading zeros */
788 2022 : LOG("0x%" PRIxBIGD, a[--len]);
789 8512 : while (len--) {
790 6490 : LOG("%08" PRIxBIGD, a[len]);
791 : }
792 2022 : if (suffix) {
793 2022 : LOG("%s", suffix);
794 : }
795 2022 : }
796 :
797 : #endif
|
