Coverage Report

Created: 2024-11-21 06:47

/src/mpdecimal-4.0.0/libmpdec/basearith.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright (c) 2008-2024 Stefan Krah. All rights reserved.
3
 *
4
 * Redistribution and use in source and binary forms, with or without
5
 * modification, are permitted provided that the following conditions
6
 * are met:
7
 *
8
 * 1. Redistributions of source code must retain the above copyright
9
 *    notice, this list of conditions and the following disclaimer.
10
 * 2. Redistributions in binary form must reproduce the above copyright
11
 *    notice, this list of conditions and the following disclaimer in the
12
 *    documentation and/or other materials provided with the distribution.
13
 *
14
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24
 * SUCH DAMAGE.
25
 */
26
27
28
#include <assert.h>
29
#include <stdio.h>
30
31
#include "basearith.h"
32
#include "constants.h"
33
#include "mpdecimal.h"
34
#include "typearith.h"
35
36
37
/*********************************************************************/
38
/*                   Calculations in base MPD_RADIX                  */
39
/*********************************************************************/
40
41
/*
42
 * Knuth, TAOCP, Volume 2, 4.3.1:
43
 *     w := sum of u (len m) and v (len n)
44
 *     n > 0 and m >= n
45
 * The calling function has to handle a possible final carry.
46
 */
47
mpd_uint_t
48
_mpd_baseadd(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
49
             mpd_size_t m, mpd_size_t n)
50
42.0k
{
51
42.0k
    mpd_uint_t s;
52
42.0k
    mpd_uint_t carry = 0;
53
42.0k
    mpd_size_t i;
54
55
42.0k
    assert(n > 0 && m >= n);
56
57
    /* add n members of u and v */
58
248k
    for (i = 0; i < n; i++) {
59
206k
        s = u[i] + (v[i] + carry);
60
206k
        carry = (s < u[i]) | (s >= MPD_RADIX);
61
206k
        w[i] = carry ? s-MPD_RADIX : s;
62
206k
    }
63
    /* if there is a carry, propagate it */
64
69.0k
    for (; carry && i < m; i++) {
65
27.0k
        s = u[i] + carry;
66
27.0k
        carry = (s == MPD_RADIX);
67
27.0k
        w[i] = carry ? 0 : s;
68
27.0k
    }
69
    /* copy the rest of u */
70
57.2k
    for (; i < m; i++) {
71
15.2k
        w[i] = u[i];
72
15.2k
    }
73
74
42.0k
    return carry;
75
42.0k
}
76
77
/*
78
 * Add the contents of u to w. Carries are propagated further. The caller
79
 * has to make sure that w is big enough.
80
 */
81
void
82
_mpd_baseaddto(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
83
0
{
84
0
    mpd_uint_t s;
85
0
    mpd_uint_t carry = 0;
86
0
    mpd_size_t i;
87
88
0
    if (n == 0) return;
89
90
    /* add n members of u to w */
91
0
    for (i = 0; i < n; i++) {
92
0
        s = w[i] + (u[i] + carry);
93
0
        carry = (s < w[i]) | (s >= MPD_RADIX);
94
0
        w[i] = carry ? s-MPD_RADIX : s;
95
0
    }
96
    /* if there is a carry, propagate it */
97
0
    for (; carry; i++) {
98
0
        s = w[i] + carry;
99
0
        carry = (s == MPD_RADIX);
100
0
        w[i] = carry ? 0 : s;
101
0
    }
102
0
}
103
104
/*
105
 * Add v to w (len m). The calling function has to handle a possible
106
 * final carry. Assumption: m > 0.
107
 */
108
mpd_uint_t
109
_mpd_shortadd(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v)
110
0
{
111
0
    mpd_uint_t s;
112
0
    mpd_uint_t carry;
113
0
    mpd_size_t i;
114
115
0
    assert(m > 0);
116
117
    /* add v to w */
118
0
    s = w[0] + v;
119
0
    carry = (s < v) | (s >= MPD_RADIX);
120
0
    w[0] = carry ? s-MPD_RADIX : s;
121
122
    /* if there is a carry, propagate it */
123
0
    for (i = 1; carry && i < m; i++) {
124
0
        s = w[i] + carry;
125
0
        carry = (s == MPD_RADIX);
126
0
        w[i] = carry ? 0 : s;
127
0
    }
128
129
0
    return carry;
130
0
}
131
132
/* Increment u. The calling function has to handle a possible carry. */
133
mpd_uint_t
134
_mpd_baseincr(mpd_uint_t *u, mpd_size_t n)
135
1.29k
{
136
1.29k
    mpd_uint_t s;
137
1.29k
    mpd_uint_t carry = 1;
138
1.29k
    mpd_size_t i;
139
140
1.29k
    assert(n > 0);
141
142
    /* if there is a carry, propagate it */
143
2.62k
    for (i = 0; carry && i < n; i++) {
144
1.32k
        s = u[i] + carry;
145
1.32k
        carry = (s == MPD_RADIX);
146
1.32k
        u[i] = carry ? 0 : s;
147
1.32k
    }
148
149
1.29k
    return carry;
150
1.29k
}
151
152
/*
153
 * Knuth, TAOCP, Volume 2, 4.3.1:
154
 *     w := difference of u (len m) and v (len n).
155
 *     number in u >= number in v;
156
 */
157
void
158
_mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
159
             mpd_size_t m, mpd_size_t n)
160
9.30k
{
161
9.30k
    mpd_uint_t d;
162
9.30k
    mpd_uint_t borrow = 0;
163
9.30k
    mpd_size_t i;
164
165
9.30k
    assert(m > 0 && n > 0);
166
167
    /* subtract n members of v from u */
168
43.6k
    for (i = 0; i < n; i++) {
169
34.3k
        d = u[i] - (v[i] + borrow);
170
34.3k
        borrow = (u[i] < d);
171
34.3k
        w[i] = borrow ? d + MPD_RADIX : d;
172
34.3k
    }
173
    /* if there is a borrow, propagate it */
174
11.7k
    for (; borrow && i < m; i++) {
175
2.44k
        d = u[i] - borrow;
176
2.44k
        borrow = (u[i] == 0);
177
2.44k
        w[i] = borrow ? MPD_RADIX-1 : d;
178
2.44k
    }
179
    /* copy the rest of u */
180
18.9k
    for (; i < m; i++) {
181
9.59k
        w[i] = u[i];
182
9.59k
    }
183
9.30k
}
184
185
/*
186
 * Subtract the contents of u from w. w is larger than u. Borrows are
187
 * propagated further, but eventually w can absorb the final borrow.
188
 */
189
void
190
_mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
191
0
{
192
0
    mpd_uint_t d;
193
0
    mpd_uint_t borrow = 0;
194
0
    mpd_size_t i;
195
196
0
    if (n == 0) return;
197
198
    /* subtract n members of u from w */
199
0
    for (i = 0; i < n; i++) {
200
0
        d = w[i] - (u[i] + borrow);
201
0
        borrow = (w[i] < d);
202
0
        w[i] = borrow ? d + MPD_RADIX : d;
203
0
    }
204
    /* if there is a borrow, propagate it */
205
0
    for (; borrow; i++) {
206
0
        d = w[i] - borrow;
207
0
        borrow = (w[i] == 0);
208
0
        w[i] = borrow ? MPD_RADIX-1 : d;
209
0
    }
210
0
}
211
212
/* w := product of u (len n) and v (single word) */
213
void
214
_mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
215
1.92M
{
216
1.92M
    mpd_uint_t hi, lo;
217
1.92M
    mpd_uint_t carry = 0;
218
1.92M
    mpd_size_t i;
219
220
1.92M
    assert(n > 0);
221
222
28.5M
    for (i=0; i < n; i++) {
223
224
26.6M
        _mpd_mul_words(&hi, &lo, u[i], v);
225
26.6M
        lo = carry + lo;
226
26.6M
        if (lo < carry) hi++;
227
228
26.6M
        _mpd_div_words_r(&carry, &w[i], hi, lo);
229
26.6M
    }
230
1.92M
    w[i] = carry;
231
1.92M
}
232
233
/*
234
 * Knuth, TAOCP, Volume 2, 4.3.1:
235
 *     w := product of u (len m) and v (len n)
236
 *     w must be initialized to zero
237
 */
238
void
239
_mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
240
             mpd_size_t m, mpd_size_t n)
241
822k
{
242
822k
    mpd_uint_t hi, lo;
243
822k
    mpd_uint_t carry;
244
822k
    mpd_size_t i, j;
245
246
822k
    assert(m > 0 && n > 0);
247
248
5.22M
    for (j=0; j < n; j++) {
249
4.40M
        carry = 0;
250
28.4M
        for (i=0; i < m; i++) {
251
252
24.0M
            _mpd_mul_words(&hi, &lo, u[i], v[j]);
253
24.0M
            lo = w[i+j] + lo;
254
24.0M
            if (lo < w[i+j]) hi++;
255
24.0M
            lo = carry + lo;
256
24.0M
            if (lo < carry) hi++;
257
258
24.0M
            _mpd_div_words_r(&carry, &w[i+j], hi, lo);
259
24.0M
        }
260
4.40M
        w[j+m] = carry;
261
4.40M
    }
262
822k
}
263
264
/*
265
 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
266
 *     w := quotient of u (len n) divided by a single word v
267
 */
268
mpd_uint_t
269
_mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
270
19.2M
{
271
19.2M
    mpd_uint_t hi, lo;
272
19.2M
    mpd_uint_t rem = 0;
273
19.2M
    mpd_size_t i;
274
275
19.2M
    assert(n > 0);
276
277
62.8M
    for (i=n-1; i != MPD_SIZE_MAX; i--) {
278
279
43.6M
        _mpd_mul_words(&hi, &lo, rem, MPD_RADIX);
280
43.6M
        lo = u[i] + lo;
281
43.6M
        if (lo < u[i]) hi++;
282
283
43.6M
        _mpd_div_words(&w[i], &rem, hi, lo, v);
284
43.6M
    }
285
286
19.2M
    return rem;
287
19.2M
}
288
289
/*
290
 * Knuth, TAOCP Volume 2, 4.3.1:
291
 *     q, r := quotient and remainder of uconst (len nplusm)
292
 *             divided by vconst (len n)
293
 *     nplusm >= n
294
 *
295
 * If r is not NULL, r will contain the remainder. If r is NULL, the
296
 * return value indicates if there is a remainder: 1 for true, 0 for
297
 * false.  A return value of -1 indicates an error.
298
 */
299
int
300
_mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r,
301
                const mpd_uint_t *uconst, const mpd_uint_t *vconst,
302
                mpd_size_t nplusm, mpd_size_t n)
303
936k
{
304
936k
    mpd_uint_t ustatic[MPD_MINALLOC_MAX];
305
936k
    mpd_uint_t vstatic[MPD_MINALLOC_MAX];
306
936k
    mpd_uint_t *u = ustatic;
307
936k
    mpd_uint_t *v = vstatic;
308
936k
    mpd_uint_t d, qhat, rhat, w2[2];
309
936k
    mpd_uint_t hi, lo, x;
310
936k
    mpd_uint_t carry;
311
936k
    mpd_size_t i, j, m;
312
936k
    int retval = 0;
313
314
936k
    assert(n > 1 && nplusm >= n);
315
936k
    m = sub_size_t(nplusm, n);
316
317
    /* D1: normalize */
318
936k
    d = MPD_RADIX / (vconst[n-1] + 1);
319
320
936k
    if (nplusm >= MPD_MINALLOC_MAX) {
321
2.46k
        if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {
322
0
            return -1;
323
0
        }
324
2.46k
    }
325
936k
    if (n >= MPD_MINALLOC_MAX) {
326
0
        if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {
327
0
            mpd_free(u);
328
0
            return -1;
329
0
        }
330
0
    }
331
332
936k
    _mpd_shortmul(u, uconst, nplusm, d);
333
936k
    _mpd_shortmul(v, vconst, n, d);
334
335
    /* D2: loop */
336
18.7M
    for (j=m; j != MPD_SIZE_MAX; j--) {
337
338
        /* D3: calculate qhat and rhat */
339
17.7M
        rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);
340
17.7M
        qhat = w2[1] * MPD_RADIX + w2[0];
341
342
18.0M
        while (1) {
343
18.0M
            if (qhat < MPD_RADIX) {
344
17.9M
                _mpd_singlemul(w2, qhat, v[n-2]);
345
17.9M
                if (w2[1] <= rhat) {
346
14.9M
                    if (w2[1] != rhat || w2[0] <= u[j+n-2]) {
347
14.9M
                        break;
348
14.9M
                    }
349
14.9M
                }
350
17.9M
            }
351
3.04M
            qhat -= 1;
352
3.04M
            rhat += v[n-1];
353
3.04M
            if (rhat < v[n-1] || rhat >= MPD_RADIX) {
354
2.80M
                break;
355
2.80M
            }
356
3.04M
        }
357
        /* D4: multiply and subtract */
358
17.7M
        carry = 0;
359
126M
        for (i=0; i <= n; i++) {
360
361
108M
            _mpd_mul_words(&hi, &lo, qhat, v[i]);
362
363
108M
            lo = carry + lo;
364
108M
            if (lo < carry) hi++;
365
366
108M
            _mpd_div_words_r(&hi, &lo, hi, lo);
367
368
108M
            x = u[i+j] - lo;
369
108M
            carry = (u[i+j] < x);
370
108M
            u[i+j] = carry ? x+MPD_RADIX : x;
371
108M
            carry += hi;
372
108M
        }
373
17.7M
        q[j] = qhat;
374
        /* D5: test remainder */
375
17.7M
        if (carry) {
376
26.6k
            q[j] -= 1;
377
            /* D6: add back */
378
26.6k
            (void)_mpd_baseadd(u+j, u+j, v, n+1, n);
379
26.6k
        }
380
17.7M
    }
381
382
    /* D8: unnormalize */
383
936k
    if (r != NULL) {
384
933k
        _mpd_shortdiv(r, u, n, d);
385
        /* we are not interested in the return value here */
386
933k
        retval = 0;
387
933k
    }
388
2.46k
    else {
389
2.46k
        retval = !_mpd_isallzero(u, n);
390
2.46k
    }
391
392
393
936k
if (u != ustatic) mpd_free(u);
394
936k
if (v != vstatic) mpd_free(v);
395
936k
return retval;
396
936k
}
397
398
/*
399
 * Left shift of src by 'shift' digits; src may equal dest.
400
 *
401
 *  dest := area of n mpd_uint_t with space for srcdigits+shift digits.
402
 *  src  := coefficient with length m.
403
 *
404
 * The case splits in the function are non-obvious. The following
405
 * equations might help:
406
 *
407
 *  Let msdigits denote the number of digits in the most significant
408
 *  word of src. Then 1 <= msdigits <= rdigits.
409
 *
410
 *   1) shift = q * rdigits + r
411
 *   2) srcdigits = qsrc * rdigits + msdigits
412
 *   3) destdigits = shift + srcdigits
413
 *                 = q * rdigits + r + qsrc * rdigits + msdigits
414
 *                 = q * rdigits + (qsrc * rdigits + (r + msdigits))
415
 *
416
 *  The result has q zero words, followed by the coefficient that is left-
417
 *  shifted by r.  The case r == 0 is trivial. For r > 0, it is important
418
 *  to keep in mind that we always read m source words, but write m+1
419
 *  destination words if r + msdigits > rdigits, m words otherwise.
420
 */
421
void
422
_mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m,
423
                mpd_size_t shift)
424
2.78k
{
425
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
426
    /* spurious uninitialized warnings */
427
    mpd_uint_t l=l, lprev=lprev, h=h;
428
#else
429
2.78k
    mpd_uint_t l, lprev, h;
430
2.78k
#endif
431
2.78k
    mpd_uint_t q, r;
432
2.78k
    mpd_uint_t ph;
433
434
2.78k
    assert(m > 0 && n >= m);
435
436
2.78k
    _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
437
438
2.78k
    if (r != 0) {
439
440
2.64k
        ph = mpd_pow10[r];
441
442
2.64k
        --m; --n;
443
2.64k
        _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);
444
2.64k
        if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */
445
1.17k
            dest[n--] = h;
446
1.17k
        }
447
        /* write m-1 shifted words */
448
11.1k
        for (; m != MPD_SIZE_MAX; m--,n--) {
449
8.55k
            _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
450
8.55k
            dest[n] = ph * lprev + h;
451
8.55k
            lprev = l;
452
8.55k
        }
453
        /* write least significant word */
454
2.64k
        dest[q] = ph * lprev;
455
2.64k
    }
456
138
    else {
457
752
        while (--m != MPD_SIZE_MAX) {
458
614
            dest[m+q] = src[m];
459
614
        }
460
138
    }
461
462
2.78k
    mpd_uint_zero(dest, q);
463
2.78k
}
464
465
/*
466
 * Right shift of src by 'shift' digits; src may equal dest.
467
 * Assumption: srcdigits-shift > 0.
468
 *
469
 *  dest := area with space for srcdigits-shift digits.
470
 *  src  := coefficient with length 'slen'.
471
 *
472
 * The case splits in the function rely on the following equations:
473
 *
474
 *  Let msdigits denote the number of digits in the most significant
475
 *  word of src. Then 1 <= msdigits <= rdigits.
476
 *
477
 *  1) shift = q * rdigits + r
478
 *  2) srcdigits = qsrc * rdigits + msdigits
479
 *  3) destdigits = srcdigits - shift
480
 *                = qsrc * rdigits + msdigits - (q * rdigits + r)
481
 *                = (qsrc - q) * rdigits + msdigits - r
482
 *
483
 * Since destdigits > 0 and 1 <= msdigits <= rdigits:
484
 *
485
 *  4) qsrc >= q
486
 *  5) qsrc == q  ==>  msdigits > r
487
 *
488
 * The result has slen-q words if msdigits > r, slen-q-1 words otherwise.
489
 */
490
mpd_uint_t
491
_mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,
492
                mpd_size_t shift)
493
2.78k
{
494
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
495
    /* spurious uninitialized warnings */
496
    mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */
497
#else
498
2.78k
    mpd_uint_t l, h, hprev; /* low, high, previous high */
499
2.78k
#endif
500
2.78k
    mpd_uint_t rnd, rest;   /* rounding digit, rest */
501
2.78k
    mpd_uint_t q, r;
502
2.78k
    mpd_size_t i, j;
503
2.78k
    mpd_uint_t ph;
504
505
2.78k
    assert(slen > 0);
506
507
2.78k
    _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
508
509
2.78k
    rnd = rest = 0;
510
2.78k
    if (r != 0) {
511
512
2.76k
        ph = mpd_pow10[MPD_RDIGITS-r];
513
514
2.76k
        _mpd_divmod_pow10(&hprev, &rest, src[q], r);
515
2.76k
        _mpd_divmod_pow10(&rnd, &rest, rest, r-1);
516
517
2.76k
        if (rest == 0 && q > 0) {
518
265
            rest = !_mpd_isallzero(src, q);
519
265
        }
520
        /* write slen-q-1 words */
521
13.1M
        for (j=0,i=q+1; i<slen; i++,j++) {
522
13.1M
            _mpd_divmod_pow10(&h, &l, src[i], r);
523
13.1M
            dest[j] = ph * l + hprev;
524
13.1M
            hprev = h;
525
13.1M
        }
526
        /* write most significant word */
527
2.76k
        if (hprev != 0) { /* always the case if slen==q-1 */
528
2.65k
            dest[j] = hprev;
529
2.65k
        }
530
2.76k
    }
531
13
    else {
532
13
        if (q > 0) {
533
13
            _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);
534
            /* is there any non-zero digit below rnd? */
535
13
            if (rest == 0) rest = !_mpd_isallzero(src, q-1);
536
13
        }
537
72
        for (j = 0; j < slen-q; j++) {
538
59
            dest[j] = src[q+j];
539
59
        }
540
13
    }
541
542
    /* 0-4  ==> rnd+rest < 0.5   */
543
    /* 5    ==> rnd+rest == 0.5  */
544
    /* 6-9  ==> rnd+rest > 0.5   */
545
2.78k
    return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd;
546
2.78k
}
547
548
/*********************************************************************/
549
/*                      Calculations in base b                       */
550
/*********************************************************************/
551
552
/*
553
 * Add v to w (len m). The calling function has to handle a possible
554
 * final carry. Assumption: m > 0.
555
 */
556
mpd_uint_t
557
_mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b)
558
0
{
559
0
    mpd_uint_t s;
560
0
    mpd_uint_t carry;
561
0
    mpd_size_t i;
562
563
0
    assert(m > 0);
564
565
    /* add v to w */
566
0
    s = w[0] + v;
567
0
    carry = (s < v) | (s >= b);
568
0
    w[0] = carry ? s-b : s;
569
570
    /* if there is a carry, propagate it */
571
0
    for (i = 1; carry && i < m; i++) {
572
0
        s = w[i] + carry;
573
0
        carry = (s == b);
574
0
        w[i] = carry ? 0 : s;
575
0
    }
576
577
0
    return carry;
578
0
}
579
580
/* w := product of u (len n) and v (single word). Return carry. */
581
mpd_uint_t
582
_mpd_shortmul_c(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
583
0
{
584
0
    mpd_uint_t hi, lo;
585
0
    mpd_uint_t carry = 0;
586
0
    mpd_size_t i;
587
588
0
    assert(n > 0);
589
590
0
    for (i=0; i < n; i++) {
591
592
0
        _mpd_mul_words(&hi, &lo, u[i], v);
593
0
        lo = carry + lo;
594
0
        if (lo < carry) hi++;
595
596
0
        _mpd_div_words_r(&carry, &w[i], hi, lo);
597
0
    }
598
599
0
    return carry;
600
0
}
601
602
/* w := product of u (len n) and v (single word) */
603
mpd_uint_t
604
_mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
605
                mpd_uint_t v, mpd_uint_t b)
606
0
{
607
0
    mpd_uint_t hi, lo;
608
0
    mpd_uint_t carry = 0;
609
0
    mpd_size_t i;
610
611
0
    assert(n > 0);
612
613
0
    for (i=0; i < n; i++) {
614
615
0
        _mpd_mul_words(&hi, &lo, u[i], v);
616
0
        lo = carry + lo;
617
0
        if (lo < carry) hi++;
618
619
0
        _mpd_div_words(&carry, &w[i], hi, lo, b);
620
0
    }
621
622
0
    return carry;
623
0
}
624
625
/*
626
 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
627
 *     w := quotient of u (len n) divided by a single word v
628
 */
629
mpd_uint_t
630
_mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
631
                mpd_uint_t v, mpd_uint_t b)
632
0
{
633
0
    mpd_uint_t hi, lo;
634
0
    mpd_uint_t rem = 0;
635
0
    mpd_size_t i;
636
637
0
    assert(n > 0);
638
639
0
    for (i=n-1; i != MPD_SIZE_MAX; i--) {
640
641
0
        _mpd_mul_words(&hi, &lo, rem, b);
642
0
        lo = u[i] + lo;
643
0
        if (lo < u[i]) hi++;
644
645
0
        _mpd_div_words(&w[i], &rem, hi, lo, v);
646
0
    }
647
648
0
    return rem;
649
0
}