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