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 : // Disable static analyzer for this third-party library file
52 : #ifndef __clang_analyzer__
53 :
54 : #define BITS_PER_HALF_DIGIT (BITS_PER_DIGIT / 2)
55 : #define LOHALF(x) ((DIGIT_T)((x) & MAX_HALF_DIGIT))
56 : #define HIHALF(x) ((DIGIT_T)((x) >> BITS_PER_HALF_DIGIT & MAX_HALF_DIGIT))
57 : #define TOHIGH(x) ((DIGIT_T)((x) << BITS_PER_HALF_DIGIT))
58 :
59 : /** Computes p = x * y, where x and y are single digits */
60 : static int spMultiply(DIGIT_T p[2], DIGIT_T x, DIGIT_T y);
61 :
62 : /** Computes quotient q = u div v, remainder r = u mod v, where q, r and v are single digits */
63 : static DIGIT_T spDivide(DIGIT_T *q, DIGIT_T *r, const DIGIT_T u[2], DIGIT_T v);
64 :
65 : /** Computes quotient q = u div d, returns remainder */
66 : static DIGIT_T mpShortDiv(DIGIT_T q[], const DIGIT_T u[], DIGIT_T d, size_t ndigits);
67 :
68 323637 : DIGIT_T mpSetZero(volatile DIGIT_T a[], size_t ndigits)
69 : { /* Sets a = 0 */
70 :
71 : /* Prevent optimiser ignoring this */
72 : volatile DIGIT_T optdummy;
73 323637 : volatile DIGIT_T *p = a;
74 :
75 2946524 : while (ndigits--)
76 2622887 : a[ndigits] = 0;
77 :
78 323637 : optdummy = *p;
79 323637 : return optdummy;
80 : }
81 :
82 881 : void mpSetDigit(DIGIT_T a[], DIGIT_T d, size_t ndigits)
83 : { /* Sets a = d where d is a single digit */
84 : size_t i;
85 :
86 7809 : for (i = 1; i < ndigits; i++)
87 : {
88 6928 : a[i] = 0;
89 : }
90 881 : a[0] = d;
91 881 : }
92 :
93 : /** Returns number of significant digits in a */
94 32156 : size_t mpSizeof(const DIGIT_T a[], size_t ndigits)
95 : {
96 35420 : while(ndigits--)
97 : {
98 35318 : if (a[ndigits] != 0)
99 32054 : return (++ndigits);
100 : }
101 102 : return 0;
102 : }
103 :
104 1184 : void mpSetEqual(DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
105 : { /* Sets a = b */
106 : size_t i;
107 :
108 9628 : for (i = 0; i < ndigits; i++)
109 : {
110 8444 : a[i] = b[i];
111 : }
112 1184 : }
113 :
114 2628 : DIGIT_T mpAdd(DIGIT_T w[], const DIGIT_T u[], const DIGIT_T v[], size_t ndigits)
115 : {
116 : /* Calculates w = u + v
117 : where w, u, v are multiprecision integers of ndigits each
118 : Returns carry if overflow. Carry = 0 or 1.
119 :
120 : Ref: Knuth Vol 2 Ch 4.3.1 p 266 Algorithm A.
121 : */
122 :
123 : DIGIT_T k;
124 : size_t j;
125 :
126 : /* w can't be the same as v
127 : Stop if assert is working, else return error (overflow = -1)
128 : */
129 2628 : if (w == v) {
130 : assert(w != v);
131 0 : return MAX_DIGIT;
132 : }
133 :
134 : /* Step A1. Initialise */
135 2628 : k = 0;
136 :
137 20280 : for (j = 0; j < ndigits; j++)
138 : {
139 : /* Step A2. Add digits w_j = (u_j + v_j + k)
140 : Set k = 1 if carry (overflow) occurs
141 : */
142 17652 : w[j] = u[j] + k;
143 17652 : if (w[j] < k)
144 0 : k = 1;
145 : else
146 17652 : k = 0;
147 :
148 17652 : w[j] += v[j];
149 17652 : if (w[j] < v[j])
150 7972 : k++;
151 :
152 : } /* Step A3. Loop on j */
153 :
154 2628 : return k; /* w_n = k */
155 : }
156 :
157 59604 : static int QhatTooBig(DIGIT_T qhat, DIGIT_T rhat,
158 : DIGIT_T vn2, DIGIT_T ujn2)
159 : { /* Returns true if Qhat is too big
160 : i.e. if (Qhat * Vn-2) > (b.Rhat + Uj+n-2)
161 : */
162 : DIGIT_T t[2];
163 :
164 59604 : spMultiply(t, qhat, vn2);
165 59604 : if (t[1] < rhat)
166 44991 : return 0;
167 14613 : else if (t[1] > rhat)
168 14613 : return 1;
169 0 : else if (t[0] > ujn2)
170 0 : return 1;
171 :
172 0 : return 0;
173 : }
174 :
175 56122 : static DIGIT_T mpMultSub(DIGIT_T wn, DIGIT_T w[], const DIGIT_T v[],
176 : DIGIT_T q, size_t n)
177 : { /* Compute w = w - qv
178 : where w = (WnW[n-1]...W[0])
179 : return modified Wn.
180 : */
181 : DIGIT_T k, t[2];
182 : size_t i;
183 :
184 56122 : if (q == 0) /* No change */
185 5533 : return wn;
186 :
187 50589 : k = 0;
188 :
189 272646 : for (i = 0; i < n; i++)
190 : {
191 222057 : spMultiply(t, q, v[i]);
192 222057 : w[i] -= k;
193 222057 : if (w[i] > MAX_DIGIT - k)
194 34852 : k = 1;
195 : else
196 187205 : k = 0;
197 222057 : w[i] -= t[0];
198 222057 : if (w[i] > MAX_DIGIT - t[0])
199 104425 : k++;
200 222057 : k += t[1];
201 : }
202 :
203 : /* Cope with Wn not stored in array w[0..n-1] */
204 50589 : wn -= k;
205 :
206 50589 : return wn;
207 : }
208 :
209 16033 : int mpDivide(DIGIT_T q[], DIGIT_T r[], const DIGIT_T u[],
210 : size_t udigits, DIGIT_T v[], size_t vdigits)
211 : { /* Computes quotient q = u / v and remainder r = u mod v
212 : where q, r, u are multiple precision digits
213 : all of udigits and the divisor v is vdigits.
214 :
215 : Ref: Knuth Vol 2 Ch 4.3.1 p 272 Algorithm D.
216 :
217 : Do without extra storage space, i.e. use r[] for
218 : normalised u[], unnormalise v[] at end, and cope with
219 : extra digit Uj+n added to u after normalisation.
220 :
221 : WARNING: this trashes q and r first, so cannot do
222 : u = u / v or v = u mod v.
223 : It also changes v temporarily so cannot make it const.
224 : */
225 : size_t shift;
226 : int n, m, j;
227 : DIGIT_T bitmask, overflow;
228 : DIGIT_T qhat, rhat, t[2];
229 : DIGIT_T *uu, *ww;
230 : int qhatOK, cmp;
231 :
232 : /* Clear q and r */
233 16033 : mpSetZero(q, udigits);
234 16033 : mpSetZero(r, udigits);
235 :
236 : /* Work out exact sizes of u and v */
237 16033 : n = (int)mpSizeof(v, vdigits);
238 16033 : m = (int)mpSizeof(u, udigits);
239 16033 : m -= n;
240 :
241 : /* Catch special cases */
242 16033 : if (n == 0)
243 78 : return -1; /* Error: divide by zero */
244 :
245 15955 : if (n == 1)
246 : { /* Use short division instead */
247 2695 : r[0] = mpShortDiv(q, u, v[0], udigits);
248 2695 : return 0;
249 : }
250 :
251 13260 : if (m < 0)
252 : { /* v > u, so just set q = 0 and r = u */
253 0 : mpSetEqual(r, u, udigits);
254 0 : return 0;
255 : }
256 :
257 13260 : if (m == 0)
258 : { /* u and v are the same length */
259 2238 : cmp = mpCompare(u, v, (size_t)n);
260 2238 : if (cmp < 0)
261 : { /* v > u, as above */
262 1160 : mpSetEqual(r, u, udigits);
263 1160 : return 0;
264 : }
265 1078 : else if (cmp == 0)
266 : { /* v == u, so set q = 1 and r = 0 */
267 0 : mpSetDigit(q, 1, udigits);
268 0 : return 0;
269 : }
270 : }
271 :
272 : /* In Knuth notation, we have:
273 : Given
274 : u = (Um+n-1 ... U1U0)
275 : v = (Vn-1 ... V1V0)
276 : Compute
277 : q = u/v = (QmQm-1 ... Q0)
278 : r = u mod v = (Rn-1 ... R1R0)
279 : */
280 :
281 : /* Step D1. Normalise */
282 : /* Requires high bit of Vn-1
283 : to be set, so find most signif. bit then shift left,
284 : i.e. d = 2^shift, u' = u * d, v' = v * d.
285 : */
286 12100 : bitmask = HIBITMASK;
287 24769 : for (shift = 0; shift < BITS_PER_DIGIT; shift++)
288 : {
289 24769 : if (v[n-1] & bitmask)
290 12100 : break;
291 12669 : bitmask >>= 1;
292 : }
293 :
294 : /* Normalise v in situ - NB only shift non-zero digits */
295 12100 : overflow = mpShiftLeft(v, v, shift, n);
296 :
297 : /* Copy normalised dividend u*d into r */
298 12100 : overflow = mpShiftLeft(r, u, shift, n + m);
299 12100 : uu = r; /* Use ptr to keep notation constant */
300 :
301 12100 : t[0] = overflow; /* Extra digit Um+n */
302 :
303 : /* Step D2. Initialise j. Set j = m */
304 68222 : for (j = m; j >= 0; j--)
305 : {
306 : /* Step D3. Set Qhat = [(b.Uj+n + Uj+n-1)/Vn-1]
307 : and Rhat = remainder */
308 56122 : qhatOK = 0;
309 56122 : t[1] = t[0]; /* This is Uj+n */
310 56122 : t[0] = uu[j+n-1];
311 56122 : overflow = spDivide(&qhat, &rhat, t, v[n-1]);
312 :
313 : /* Test Qhat */
314 56122 : if (overflow)
315 : { /* Qhat == b so set Qhat = b - 1 */
316 0 : qhat = MAX_DIGIT;
317 0 : rhat = uu[j+n-1];
318 0 : rhat += v[n-1];
319 0 : if (rhat < v[n-1]) /* Rhat >= b, so no re-test */
320 0 : qhatOK = 1;
321 : }
322 : /* [VERSION 2: Added extra test "qhat && "] */
323 56122 : if (qhat && !qhatOK && QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
324 : { /* If Qhat.Vn-2 > b.Rhat + Uj+n-2
325 : decrease Qhat by one, increase Rhat by Vn-1
326 : */
327 14249 : qhat--;
328 14249 : rhat += v[n-1];
329 : /* Repeat this test if Rhat < b */
330 14249 : if (!(rhat < v[n-1]))
331 9015 : if (QhatTooBig(qhat, rhat, v[n-2], uu[j+n-2]))
332 364 : qhat--;
333 : }
334 :
335 :
336 : /* Step D4. Multiply and subtract */
337 56122 : ww = &uu[j];
338 56122 : overflow = mpMultSub(t[1], ww, v, qhat, (size_t)n);
339 :
340 : /* Step D5. Test remainder. Set Qj = Qhat */
341 56122 : q[j] = qhat;
342 56122 : if (overflow)
343 : { /* Step D6. Add back if D4 was negative */
344 0 : q[j]--;
345 0 : overflow = mpAdd(ww, ww, v, (size_t)n);
346 : }
347 :
348 56122 : t[0] = uu[j+n-1]; /* Uj+n on next round */
349 :
350 : } /* Step D7. Loop on j */
351 :
352 : /* Clear high digits in uu */
353 56122 : for (j = n; j < m+n; j++)
354 44022 : uu[j] = 0;
355 :
356 : /* Step D8. Unnormalise. */
357 :
358 12100 : mpShiftRight(r, r, shift, n);
359 12100 : mpShiftRight(v, v, shift, n);
360 :
361 12100 : return 0;
362 : }
363 :
364 2695 : static DIGIT_T mpShortDiv(DIGIT_T q[], const DIGIT_T u[], DIGIT_T v,
365 : size_t ndigits)
366 : {
367 : /* Calculates quotient q = u div v
368 : Returns remainder r = u mod v
369 : where q, u are multiprecision integers of ndigits each
370 : and r, v are single precision digits.
371 :
372 : Makes no assumptions about normalisation.
373 :
374 : Ref: Knuth Vol 2 Ch 4.3.1 Exercise 16 p625
375 : */
376 : size_t j;
377 : DIGIT_T t[2], r;
378 : size_t shift;
379 : DIGIT_T bitmask, overflow, *uu;
380 :
381 2695 : if (ndigits == 0) return 0;
382 2695 : if (v == 0) return 0; /* Divide by zero error */
383 :
384 : /* Normalise first */
385 : /* Requires high bit of V
386 : to be set, so find most signif. bit then shift left,
387 : i.e. d = 2^shift, u' = u * d, v' = v * d.
388 : */
389 2695 : bitmask = HIBITMASK;
390 13756 : for (shift = 0; shift < BITS_PER_DIGIT; shift++)
391 : {
392 13756 : if (v & bitmask)
393 2695 : break;
394 11061 : bitmask >>= 1;
395 : }
396 :
397 2695 : v <<= shift;
398 2695 : overflow = mpShiftLeft(q, u, shift, ndigits);
399 2695 : uu = q;
400 :
401 : /* Step S1 - modified for extra digit. */
402 2695 : r = overflow; /* New digit Un */
403 2695 : j = ndigits;
404 20950 : while (j--)
405 : {
406 : /* Step S2. */
407 18255 : t[1] = r;
408 18255 : t[0] = uu[j];
409 18255 : overflow = spDivide(&q[j], &r, t, v);
410 : }
411 :
412 : /* Unnormalise */
413 2695 : r >>= shift;
414 :
415 2695 : return r;
416 : }
417 :
418 210495 : DIGIT_T mpShiftLeft(DIGIT_T a[], const DIGIT_T *b,
419 : size_t shift, size_t ndigits)
420 : { /* Computes a = b << shift */
421 : /* [v2.1] Modified to cope with shift > BITS_PERDIGIT */
422 : size_t i, y, nw, bits;
423 : DIGIT_T mask, carry, nextcarry;
424 :
425 : /* Do we shift whole digits? */
426 210495 : if (shift >= BITS_PER_DIGIT)
427 : {
428 59400 : nw = shift / BITS_PER_DIGIT;
429 59400 : i = ndigits;
430 630600 : while (i--)
431 : {
432 571200 : if (i >= nw)
433 494600 : a[i] = b[i-nw];
434 : else
435 76600 : a[i] = 0;
436 : }
437 : /* Call again to shift bits inside digits */
438 59400 : bits = shift % BITS_PER_DIGIT;
439 59400 : carry = b[ndigits-nw] << bits;
440 59400 : if (bits)
441 56400 : carry |= mpShiftLeft(a, a, bits, ndigits);
442 59400 : return carry;
443 : }
444 : else
445 : {
446 151095 : bits = shift;
447 : }
448 :
449 : /* Construct mask = high bits set */
450 151095 : mask = ~(~(DIGIT_T)0 >> bits);
451 :
452 151095 : y = BITS_PER_DIGIT - bits;
453 151095 : carry = 0;
454 1367432 : for (i = 0; i < ndigits; i++)
455 : {
456 1216337 : nextcarry = (b[i] & mask) >> y;
457 1216337 : a[i] = b[i] << bits | carry;
458 1216337 : carry = nextcarry;
459 : }
460 :
461 151095 : return carry;
462 : }
463 :
464 207800 : DIGIT_T mpShiftRight(DIGIT_T a[], const DIGIT_T b[], size_t shift, size_t ndigits)
465 : { /* Computes a = b >> shift */
466 : /* [v2.1] Modified to cope with shift > BITS_PERDIGIT */
467 : size_t i, y, nw, bits;
468 : DIGIT_T mask, carry, nextcarry;
469 :
470 : /* Do we shift whole digits? */
471 207800 : if (shift >= BITS_PER_DIGIT)
472 : {
473 59400 : nw = shift / BITS_PER_DIGIT;
474 630600 : for (i = 0; i < ndigits; i++)
475 : {
476 571200 : if ((i+nw) < ndigits)
477 494600 : a[i] = b[i+nw];
478 : else
479 76600 : a[i] = 0;
480 : }
481 : /* Call again to shift bits inside digits */
482 59400 : bits = shift % BITS_PER_DIGIT;
483 59400 : carry = b[nw-1] >> bits;
484 59400 : if (bits)
485 56400 : carry |= mpShiftRight(a, a, bits, ndigits);
486 59400 : return carry;
487 : }
488 : else
489 : {
490 148400 : bits = shift;
491 : }
492 :
493 : /* Construct mask to set low bits */
494 : /* (thanks to Jesse Chisholm for suggesting this improved technique) */
495 148400 : mask = ~(~(DIGIT_T)0 << bits);
496 :
497 148400 : y = BITS_PER_DIGIT - bits;
498 148400 : carry = 0;
499 148400 : i = ndigits;
500 1302460 : while (i--)
501 : {
502 1154060 : nextcarry = (b[i] & mask) << y;
503 1154060 : a[i] = b[i] >> bits | carry;
504 1154060 : carry = nextcarry;
505 : }
506 :
507 148400 : return carry;
508 : }
509 :
510 : /** Returns sign of (a - b) as 0, +1 or -1 (constant-time) */
511 1532 : int mpCompare_ct(const DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
512 : {
513 : /* All these vars are either 0 or 1 */
514 1532 : unsigned int gt = 0;
515 1532 : unsigned int lt = 0;
516 1532 : unsigned int mask = 1; /* Set to zero once first inequality found */
517 : unsigned int c;
518 :
519 14940 : while (ndigits--) {
520 13408 : gt |= (a[ndigits] > b[ndigits]) & mask;
521 13408 : lt |= (a[ndigits] < b[ndigits]) & mask;
522 13408 : c = (gt | lt);
523 13408 : mask &= (c-1); /* Unchanged if c==0 or mask==0, else mask=0 */
524 : }
525 :
526 1532 : return (int)gt - (int)lt; /* EQ=0 GT=+1 LT=-1 */
527 : }
528 :
529 : /** Returns sign of (a - b) as 0, +1 or -1. Not constant-time. */
530 2628 : int mpCompare(const DIGIT_T a[], const DIGIT_T b[], size_t ndigits)
531 : {
532 : /* if (ndigits == 0) return 0; // deleted [v2.5] */
533 :
534 4006 : while (ndigits--)
535 : {
536 3904 : if (a[ndigits] > b[ndigits])
537 1222 : return 1; /* GT */
538 2682 : if (a[ndigits] < b[ndigits])
539 1304 : return -1; /* LT */
540 : }
541 :
542 102 : return 0; /* EQ */
543 : }
544 :
545 281661 : static int spMultiply(DIGIT_T p[2], DIGIT_T x, DIGIT_T y)
546 : { /* Computes p = x * y */
547 : /* Ref: Arbitrary Precision Computation
548 : http://numbers.computation.free.fr/Constants/constants.html
549 :
550 : high p1 p0 low
551 : +--------+--------+--------+--------+
552 : | x1*y1 | x0*y0 |
553 : +--------+--------+--------+--------+
554 : +-+--------+--------+
555 : |1| (x0*y1 + x1*y1) |
556 : +-+--------+--------+
557 : ^carry from adding (x0*y1+x1*y1) together
558 : +-+
559 : |1|< carry from adding LOHALF t
560 : +-+ to high half of p0
561 : */
562 : DIGIT_T x0, y0, x1, y1;
563 : DIGIT_T t, u, carry;
564 :
565 : /* Split each x,y into two halves
566 : x = x0 + B*x1
567 : y = y0 + B*y1
568 : where B = 2^16, half the digit size
569 : Product is
570 : xy = x0y0 + B(x0y1 + x1y0) + B^2(x1y1)
571 : */
572 :
573 281661 : x0 = LOHALF(x);
574 281661 : x1 = HIHALF(x);
575 281661 : y0 = LOHALF(y);
576 281661 : y1 = HIHALF(y);
577 :
578 : /* Calc low part - no carry */
579 281661 : p[0] = x0 * y0;
580 :
581 : /* Calc middle part */
582 281661 : t = x0 * y1;
583 281661 : u = x1 * y0;
584 281661 : t += u;
585 281661 : if (t < u)
586 19928 : carry = 1;
587 : else
588 261733 : carry = 0;
589 :
590 : /* This carry will go to high half of p[1]
591 : + high half of t into low half of p[1] */
592 281661 : carry = TOHIGH(carry) + HIHALF(t);
593 :
594 : /* Add low half of t to high half of p[0] */
595 281661 : t = TOHIGH(t);
596 281661 : p[0] += t;
597 281661 : if (p[0] < t)
598 60273 : carry++;
599 :
600 281661 : p[1] = x1 * y1;
601 281661 : p[1] += carry;
602 :
603 :
604 281661 : return 0;
605 : }
606 :
607 : /* spDivide */
608 :
609 : #define B (MAX_HALF_DIGIT + 1)
610 :
611 148754 : static void spMultSub(DIGIT_T uu[2], DIGIT_T qhat, DIGIT_T v1, DIGIT_T v0)
612 : {
613 : /* Compute uu = uu - q(v1v0)
614 : where uu = u3u2u1u0, u3 = 0
615 : and u_n, v_n are all half-digits
616 : even though v1, v2 are passed as full digits.
617 : */
618 : DIGIT_T p0, p1, t;
619 :
620 148754 : p0 = qhat * v0;
621 148754 : p1 = qhat * v1;
622 148754 : t = p0 + TOHIGH(LOHALF(p1));
623 148754 : uu[0] -= t;
624 148754 : if (uu[0] > MAX_DIGIT - t)
625 50318 : uu[1]--; /* Borrow */
626 148754 : uu[1] -= HIHALF(p1);
627 148754 : }
628 :
629 74377 : static DIGIT_T spDivide(DIGIT_T *q, DIGIT_T *r, const DIGIT_T u[2], DIGIT_T v)
630 : { /* Computes quotient q = u / v, remainder r = u mod v
631 : where u is a double digit
632 : and q, v, r are single precision digits.
633 : Returns high digit of quotient (max value is 1)
634 : CAUTION: Assumes normalised such that v1 >= b/2
635 : where b is size of HALF_DIGIT
636 : i.e. the most significant bit of v should be one
637 :
638 : In terms of half-digits in Knuth notation:
639 : (q2q1q0) = (u4u3u2u1u0) / (v1v0)
640 : (r1r0) = (u4u3u2u1u0) mod (v1v0)
641 : for m = 2, n = 2 where u4 = 0
642 : q2 is either 0 or 1.
643 : We set q = (q1q0) and return q2 as "overflow"
644 : */
645 : DIGIT_T qhat, rhat, t, v0, v1, u0, u1, u2, u3;
646 : DIGIT_T uu[2], q2;
647 :
648 : /* Check for normalisation */
649 74377 : if (!(v & HIBITMASK))
650 : { /* Stop if assert is working, else return error */
651 : assert(v & HIBITMASK);
652 0 : *q = *r = 0;
653 0 : return MAX_DIGIT;
654 : }
655 :
656 : /* Split up into half-digits */
657 74377 : v0 = LOHALF(v);
658 74377 : v1 = HIHALF(v);
659 74377 : u0 = LOHALF(u[0]);
660 74377 : u1 = HIHALF(u[0]);
661 74377 : u2 = LOHALF(u[1]);
662 74377 : u3 = HIHALF(u[1]);
663 :
664 : /* Do three rounds of Knuth Algorithm D Vol 2 p272 */
665 :
666 : /* ROUND 1. Set j = 2 and calculate q2 */
667 : /* Estimate qhat = (u4u3)/v1 = 0 or 1
668 : then set (u4u3u2) -= qhat(v1v0)
669 : where u4 = 0.
670 : */
671 : /* [Replaced in Version 2] -->
672 : qhat = u3 / v1;
673 : if (qhat > 0)
674 : {
675 : rhat = u3 - qhat * v1;
676 : t = TOHIGH(rhat) | u2;
677 : if (qhat * v0 > t)
678 : qhat--;
679 : }
680 : <-- */
681 74377 : qhat = (u3 < v1 ? 0 : 1);
682 74377 : if (qhat > 0)
683 : { /* qhat is one, so no need to mult */
684 0 : rhat = u3 - v1;
685 : /* t = r.b + u2 */
686 0 : t = TOHIGH(rhat) | u2;
687 0 : if (v0 > t)
688 0 : qhat--;
689 : }
690 :
691 74377 : uu[1] = 0; /* (u4) */
692 74377 : uu[0] = u[1]; /* (u3u2) */
693 74377 : if (qhat > 0)
694 : {
695 : /* (u4u3u2) -= qhat(v1v0) where u4 = 0 */
696 0 : spMultSub(uu, qhat, v1, v0);
697 0 : if (HIHALF(uu[1]) != 0)
698 : { /* Add back */
699 0 : qhat--;
700 0 : uu[0] += v;
701 0 : uu[1] = 0;
702 : }
703 : }
704 74377 : q2 = qhat;
705 :
706 : /* ROUND 2. Set j = 1 and calculate q1 */
707 : /* Estimate qhat = (u3u2) / v1
708 : then set (u3u2u1) -= qhat(v1v0)
709 : */
710 74377 : t = uu[0];
711 74377 : qhat = t / v1;
712 74377 : rhat = t - qhat * v1;
713 : /* Test on v0 */
714 74377 : t = TOHIGH(rhat) | u1;
715 74377 : if ((qhat == B) || (qhat * v0 > t))
716 : {
717 19027 : qhat--;
718 19027 : rhat += v1;
719 19027 : t = TOHIGH(rhat) | u1;
720 19027 : if ((rhat < B) && (qhat * v0 > t))
721 506 : qhat--;
722 : }
723 :
724 : /* Multiply and subtract
725 : (u3u2u1)' = (u3u2u1) - qhat(v1v0)
726 : */
727 74377 : uu[1] = HIHALF(uu[0]); /* (0u3) */
728 74377 : uu[0] = TOHIGH(LOHALF(uu[0])) | u1; /* (u2u1) */
729 74377 : spMultSub(uu, qhat, v1, v0);
730 74377 : if (HIHALF(uu[1]) != 0)
731 : { /* Add back */
732 0 : qhat--;
733 0 : uu[0] += v;
734 0 : uu[1] = 0;
735 : }
736 :
737 : /* q1 = qhat */
738 74377 : *q = TOHIGH(qhat);
739 :
740 : /* ROUND 3. Set j = 0 and calculate q0 */
741 : /* Estimate qhat = (u2u1) / v1
742 : then set (u2u1u0) -= qhat(v1v0)
743 : */
744 74377 : t = uu[0];
745 74377 : qhat = t / v1;
746 74377 : rhat = t - qhat * v1;
747 : /* Test on v0 */
748 74377 : t = TOHIGH(rhat) | u0;
749 74377 : if ((qhat == B) || (qhat * v0 > t))
750 : {
751 19307 : qhat--;
752 19307 : rhat += v1;
753 19307 : t = TOHIGH(rhat) | u0;
754 19307 : if ((rhat < B) && (qhat * v0 > t))
755 527 : qhat--;
756 : }
757 :
758 : /* Multiply and subtract
759 : (u2u1u0)" = (u2u1u0)' - qhat(v1v0)
760 : */
761 74377 : uu[1] = HIHALF(uu[0]); /* (0u2) */
762 74377 : uu[0] = TOHIGH(LOHALF(uu[0])) | u0; /* (u1u0) */
763 74377 : spMultSub(uu, qhat, v1, v0);
764 74377 : if (HIHALF(uu[1]) != 0)
765 : { /* Add back */
766 0 : qhat--;
767 0 : uu[0] += v;
768 0 : uu[1] = 0;
769 : }
770 :
771 : /* q0 = qhat */
772 74377 : *q |= LOHALF(qhat);
773 :
774 : /* Remainder is in (u1u0) i.e. uu[0] */
775 74377 : *r = uu[0];
776 74377 : return q2;
777 : }
778 :
779 : #endif // __clang_analyzer__
780 :
781 : // Platform-dependent code
782 : #ifndef HSM_PLATFORM_LEDGER
783 :
784 : #include "hal/log.h"
785 :
786 2022 : void LOG_BIGD_HEX(const char *prefix,
787 : const DIGIT_T *a,
788 : size_t len,
789 : const char *suffix) {
790 2022 : if (prefix) {
791 2022 : LOG("%s", prefix);
792 : }
793 : /* Trim leading digits which are zero */
794 11728 : while (len--) {
795 11708 : if (a[len] != 0)
796 2002 : break;
797 : }
798 2022 : len++;
799 2022 : if (0 == len)
800 20 : len = 1;
801 : /* print first digit without leading zeros */
802 2022 : LOG("0x%" PRIxBIGD, a[--len]);
803 8512 : while (len--) {
804 6490 : LOG("%08" PRIxBIGD, a[len]);
805 : }
806 2022 : if (suffix) {
807 2022 : LOG("%s", suffix);
808 : }
809 2022 : }
810 :
811 : #endif
|
