Rietveld Code Review Tool
Help | Bug tracker | Discussion group | Source code | Sign in
(10440)

Delta Between Two Patch Sets: Modules/_decimal/libmpdec/basearith.c

Issue 7652: Merge C version of decimal into py3k.
Left Patch Set: Created 7 years, 5 months ago
Right Patch Set: Created 7 years, 4 months ago
Left:
Right:
Use n/p to move between diff chunks; N/P to move between comments. Please Sign in to add in-line comments.
Jump to:
Right: Side by side diff | Download
« no previous file with change/comment | « Modules/_decimal/ISSUES.txt ('k') | Modules/_decimal/libmpdec/basearith.h » ('j') | no next file with change/comment »
Toggle Intra-line Diffs ('i') | Expand Comments ('e') | Collapse Comments ('c') | Show Comments Hide Comments ('s')
LEFTRIGHT
(no file at all)
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
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
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. Assumption: m > 0.
111 */
112 mpd_uint_t
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;
117 mpd_size_t i;
118
119 assert(m > 0);
120
121 /* add v to w */
122 s = w[0] + v;
123 carry = (s < v) | (s >= MPD_RADIX);
124 w[0] = carry ? s-MPD_RADIX : s;
125
126 /* if there is a carry, propagate it */
127 for (i = 1; carry && i < m; i++) {
128 s = w[i] + carry;
129 carry = (s == MPD_RADIX);
130 w[i] = carry ? 0 : s;
131 }
132
133 return carry;
134 }
135
136 /* Increment u. The calling function has to handle a possible carry. */
137 mpd_uint_t
138 _mpd_baseincr(mpd_uint_t *u, mpd_size_t n)
139 {
140 mpd_uint_t s;
141 mpd_uint_t carry = 1;
142 mpd_size_t i;
143
144 assert(n > 0);
145
146 /* if there is a carry, propagate it */
147 for (i = 0; carry && i < n; i++) {
148 s = u[i] + carry;
149 carry = (s == MPD_RADIX);
150 u[i] = carry ? 0 : s;
151 }
152
153 return carry;
154 }
155
156 /*
157 * Knuth, TAOCP, Volume 2, 4.3.1:
158 * w := difference of u (len m) and v (len n).
159 * number in u >= number in v;
160 */
161 void
162 _mpd_basesub(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
163 mpd_size_t m, mpd_size_t n)
164 {
165 mpd_uint_t d;
166 mpd_uint_t borrow = 0;
167 mpd_size_t i;
168
169 assert(m > 0 && n > 0);
170
171 /* subtract n members of v from u */
172 for (i = 0; i < n; i++) {
173 d = u[i] - (v[i] + borrow);
174 borrow = (u[i] < d);
175 w[i] = borrow ? d + MPD_RADIX : d;
176 }
177 /* if there is a borrow, propagate it */
178 for (; borrow && i < m; i++) {
179 d = u[i] - borrow;
180 borrow = (u[i] == 0);
181 w[i] = borrow ? MPD_RADIX-1 : d;
182 }
183 /* copy the rest of u */
184 for (; i < m; i++) {
185 w[i] = u[i];
186 }
187 }
188
189 /*
190 * Subtract the contents of u from w. w is larger than u. Borrows are
191 * propagated further, but eventually w can absorb the final borrow.
192 */
193 void
194 _mpd_basesubfrom(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n)
195 {
196 mpd_uint_t d;
197 mpd_uint_t borrow = 0;
198 mpd_size_t i;
199
200 if (n == 0) return;
201
202 /* subtract n members of u from w */
203 for (i = 0; i < n; i++) {
204 d = w[i] - (u[i] + borrow);
205 borrow = (w[i] < d);
206 w[i] = borrow ? d + MPD_RADIX : d;
207 }
208 /* if there is a borrow, propagate it */
209 for (; borrow; i++) {
210 d = w[i] - borrow;
211 borrow = (w[i] == 0);
212 w[i] = borrow ? MPD_RADIX-1 : d;
213 }
214 }
215
216 /* w := product of u (len n) and v (single word) */
217 void
218 _mpd_shortmul(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
219 {
220 mpd_uint_t hi, lo;
221 mpd_uint_t carry = 0;
222 mpd_size_t i;
223
224 assert(n > 0);
225
226 for (i=0; i < n; i++) {
227
228 _mpd_mul_words(&hi, &lo, u[i], v);
229 lo = carry + lo;
230 if (lo < carry) hi++;
231
232 _mpd_div_words_r(&carry, &w[i], hi, lo);
233 }
234 w[i] = carry;
235 }
236
237 /*
238 * Knuth, TAOCP, Volume 2, 4.3.1:
239 * w := product of u (len m) and v (len n)
240 * w must be initialized to zero
241 */
242 void
243 _mpd_basemul(mpd_uint_t *w, const mpd_uint_t *u, const mpd_uint_t *v,
244 mpd_size_t m, mpd_size_t n)
245 {
246 mpd_uint_t hi, lo;
247 mpd_uint_t carry;
248 mpd_size_t i, j;
249
250 assert(m > 0 && n > 0);
251
252 for (j=0; j < n; j++) {
253 carry = 0;
254 for (i=0; i < m; i++) {
255
256 _mpd_mul_words(&hi, &lo, u[i], v[j]);
257 lo = w[i+j] + lo;
258 if (lo < w[i+j]) hi++;
259 lo = carry + lo;
260 if (lo < carry) hi++;
261
262 _mpd_div_words_r(&carry, &w[i+j], hi, lo);
263 }
264 w[j+m] = carry;
265 }
266 }
267
268 /*
269 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
270 * w := quotient of u (len n) divided by a single word v
271 */
272 mpd_uint_t
273 _mpd_shortdiv(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n, mpd_uint_t v)
274 {
275 mpd_uint_t hi, lo;
276 mpd_uint_t rem = 0;
277 mpd_size_t i;
278
279 assert(n > 0);
280
281 for (i=n-1; i != MPD_SIZE_MAX; i--) {
282
283 _mpd_mul_words(&hi, &lo, rem, MPD_RADIX);
284 lo = u[i] + lo;
285 if (lo < u[i]) hi++;
286
287 _mpd_div_words(&w[i], &rem, hi, lo, v);
288 }
289
290 return rem;
291 }
292
293 /*
294 * Knuth, TAOCP Volume 2, 4.3.1:
295 * q, r := quotient and remainder of uconst (len nplusm)
296 * divided by vconst (len n)
297 * nplusm >= n
298 *
299 * If r is not NULL, r will contain the remainder. If r is NULL, the
300 * return value indicates if there is a remainder: 1 for true, 0 for
301 * false. A return value of -1 indicates an error.
302 */
303 int
304 _mpd_basedivmod(mpd_uint_t *q, mpd_uint_t *r,
305 const mpd_uint_t *uconst, const mpd_uint_t *vconst,
306 mpd_size_t nplusm, mpd_size_t n)
307 {
308 mpd_uint_t ustatic[MPD_MINALLOC_MAX];
309 mpd_uint_t vstatic[MPD_MINALLOC_MAX];
310 mpd_uint_t *u = ustatic;
311 mpd_uint_t *v = vstatic;
312 mpd_uint_t d, qhat, rhat, w2[2];
313 mpd_uint_t hi, lo, x;
314 mpd_uint_t carry;
315 mpd_size_t i, j, m;
316 int retval = 0;
317
318 assert(n > 1 && nplusm >= n);
319 m = sub_size_t(nplusm, n);
320
321 /* D1: normalize */
322 d = MPD_RADIX / (vconst[n-1] + 1);
323
324 if (nplusm >= MPD_MINALLOC_MAX) {
325 if ((u = mpd_alloc(nplusm+1, sizeof *u)) == NULL) {
326 return -1;
327 }
328 }
329 if (n >= MPD_MINALLOC_MAX) {
330 if ((v = mpd_alloc(n+1, sizeof *v)) == NULL) {
331 mpd_free(u);
332 return -1;
333 }
334 }
335
336 _mpd_shortmul(u, uconst, nplusm, d);
337 _mpd_shortmul(v, vconst, n, d);
338
339 /* D2: loop */
340 for (j=m; j != MPD_SIZE_MAX; j--) {
341
342 /* D3: calculate qhat and rhat */
343 rhat = _mpd_shortdiv(w2, u+j+n-1, 2, v[n-1]);
344 qhat = w2[1] * MPD_RADIX + w2[0];
345
346 while (1) {
347 if (qhat < MPD_RADIX) {
348 _mpd_singlemul(w2, qhat, v[n-2]);
349 if (w2[1] <= rhat) {
350 if (w2[1] != rhat || w2[0] <= u[j+n-2]) {
351 break;
352 }
353 }
354 }
355 qhat -= 1;
356 rhat += v[n-1];
357 if (rhat < v[n-1] || rhat >= MPD_RADIX) {
358 break;
359 }
360 }
361 /* D4: multiply and subtract */
362 carry = 0;
363 for (i=0; i <= n; i++) {
364
365 _mpd_mul_words(&hi, &lo, qhat, v[i]);
366
367 lo = carry + lo;
368 if (lo < carry) hi++;
369
370 _mpd_div_words_r(&hi, &lo, hi, lo);
371
372 x = u[i+j] - lo;
373 carry = (u[i+j] < x);
374 u[i+j] = carry ? x+MPD_RADIX : x;
375 carry += hi;
376 }
377 q[j] = qhat;
378 /* D5: test remainder */
379 if (carry) {
380 q[j] -= 1;
381 /* D6: add back */
382 (void)_mpd_baseadd(u+j, u+j, v, n+1, n);
383 }
384 }
385
386 /* D8: unnormalize */
387 if (r != NULL) {
388 _mpd_shortdiv(r, u, n, d);
389 /* we are not interested in the return value here */
390 retval = 0;
391 }
392 else {
393 retval = !_mpd_isallzero(u, n);
394 }
395
396
397 if (u != ustatic) mpd_free(u);
398 if (v != vstatic) mpd_free(v);
399 return retval;
400 }
401
402 /*
403 * Left shift of src by 'shift' digits; src may equal dest.
404 *
405 * dest := area of n mpd_uint_t with space for srcdigits+shift digits.
406 * src := coefficient with length m.
407 *
408 * The case splits in the function are non-obvious. The following
409 * equations might help:
410 *
411 * Let msdigits denote the number of digits in the most significant
412 * word of src. Then 1 <= msdigits <= rdigits.
413 *
414 * 1) shift = q * rdigits + r
415 * 2) srcdigits = qsrc * rdigits + msdigits
416 * 3) destdigits = shift + srcdigits
417 * = q * rdigits + r + qsrc * rdigits + msdigits
418 * = q * rdigits + (qsrc * rdigits + (r + msdigits))
419 *
420 * The result has q zero words, followed by the coefficient that
421 * is left-shifted by r. The case r == 0 is trivial. For r > 0, it
422 * is important to keep in mind that we always read m source words,
423 * but write m+1 destination words if r + msdigits > rdigits, m words
424 * otherwise.
425 */
426 void
427 _mpd_baseshiftl(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t n, mpd_size_t m,
428 mpd_size_t shift)
429 {
430 #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
431 /* spurious uninitialized warnings */
432 mpd_uint_t l=l, lprev=lprev, h=h;
433 #else
434 mpd_uint_t l, lprev, h;
435 #endif
436 mpd_uint_t q, r;
437 mpd_uint_t ph;
438
439 assert(m > 0 && n >= m);
440
441 _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
442
443 if (r != 0) {
444
445 ph = mpd_pow10[r];
446
447 --m; --n;
448 _mpd_divmod_pow10(&h, &lprev, src[m--], MPD_RDIGITS-r);
449 if (h != 0) { /* r + msdigits > rdigits <==> h != 0 */
450 dest[n--] = h;
451 }
452 /* write m-1 shifted words */
453 for (; m != MPD_SIZE_MAX; m--,n--) {
454 _mpd_divmod_pow10(&h, &l, src[m], MPD_RDIGITS-r);
455 dest[n] = ph * lprev + h;
456 lprev = l;
457 }
458 /* write least significant word */
459 dest[q] = ph * lprev;
460 }
461 else {
462 while (--m != MPD_SIZE_MAX) {
463 dest[m+q] = src[m];
464 }
465 }
466
467 mpd_uint_zero(dest, q);
468 }
469
470 /*
471 * Right shift of src by 'shift' digits; src may equal dest.
472 * Assumption: srcdigits-shift > 0.
473 *
474 * dest := area with space for srcdigits-shift digits.
475 * src := coefficient with length 'slen'.
476 *
477 * The case splits in the function rely on the following equations:
478 *
479 * Let msdigits denote the number of digits in the most significant
480 * word of src. Then 1 <= msdigits <= rdigits.
481 *
482 * 1) shift = q * rdigits + r
483 * 2) srcdigits = qsrc * rdigits + msdigits
484 * 3) destdigits = srcdigits - shift
485 * = qsrc * rdigits + msdigits - (q * rdigits + r)
486 * = (qsrc - q) * rdigits + msdigits - r
487 *
488 * Since destdigits > 0 and 1 <= msdigits <= rdigits:
489 *
490 * 4) qsrc >= q
491 * 5) qsrc == q ==> msdigits > r
492 *
493 * The result has slen-q words if msdigits > r, slen-q-1 words otherwise.
494 */
495 mpd_uint_t
496 _mpd_baseshiftr(mpd_uint_t *dest, mpd_uint_t *src, mpd_size_t slen,
497 mpd_size_t shift)
498 {
499 #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
500 /* spurious uninitialized warnings */
501 mpd_uint_t l=l, h=h, hprev=hprev; /* low, high, previous high */
502 #else
503 mpd_uint_t l, h, hprev; /* low, high, previous high */
504 #endif
505 mpd_uint_t rnd, rest; /* rounding digit, rest */
506 mpd_uint_t q, r;
507 mpd_size_t i, j;
508 mpd_uint_t ph;
509
510 assert(slen > 0);
511
512 _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
513
514 rnd = rest = 0;
515 if (r != 0) {
516
517 ph = mpd_pow10[MPD_RDIGITS-r];
518
519 _mpd_divmod_pow10(&hprev, &rest, src[q], r);
520 _mpd_divmod_pow10(&rnd, &rest, rest, r-1);
521
522 if (rest == 0 && q > 0) {
523 rest = !_mpd_isallzero(src, q);
524 }
525 /* write slen-q-1 words */
526 for (j=0,i=q+1; i<slen; i++,j++) {
527 _mpd_divmod_pow10(&h, &l, src[i], r);
528 dest[j] = ph * l + hprev;
529 hprev = h;
530 }
531 /* write most significant word */
532 if (hprev != 0) { /* always the case if slen==q-1 */
533 dest[j] = hprev;
534 }
535 }
536 else {
537 if (q > 0) {
538 _mpd_divmod_pow10(&rnd, &rest, src[q-1], MPD_RDIGITS-1);
539 /* is there any non-zero digit below rnd? */
540 if (rest == 0) rest = !_mpd_isallzero(src, q-1);
541 }
542 for (j = 0; j < slen-q; j++) {
543 dest[j] = src[q+j];
544 }
545 }
546
547 /* 0-4 ==> rnd+rest < 0.5 */
548 /* 5 ==> rnd+rest == 0.5 */
549 /* 6-9 ==> rnd+rest > 0.5 */
550 return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd;
551 }
552
553
554 /*********************************************************************/
555 /* Calculations in base b */
556 /*********************************************************************/
557
558 /*
559 * Add v to w (len m). The calling function has to handle a possible
560 * final carry. Assumption: m > 0.
561 */
562 mpd_uint_t
563 _mpd_shortadd_b(mpd_uint_t *w, mpd_size_t m, mpd_uint_t v, mpd_uint_t b)
564 {
565 mpd_uint_t s;
566 mpd_uint_t carry;
567 mpd_size_t i;
568
569 assert(m > 0);
570
571 /* add v to w */
572 s = w[0] + v;
573 carry = (s < v) | (s >= b);
574 w[0] = carry ? s-b : s;
575
576 /* if there is a carry, propagate it */
577 for (i = 1; carry && i < m; i++) {
578 s = w[i] + carry;
579 carry = (s == b);
580 w[i] = carry ? 0 : s;
581 }
582
583 return carry;
584 }
585
586 /* w := product of u (len n) and v (single word) */
587 void
588 _mpd_shortmul_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
589 mpd_uint_t v, mpd_uint_t b)
590 {
591 mpd_uint_t hi, lo;
592 mpd_uint_t carry = 0;
593 mpd_size_t i;
594
595 assert(n > 0);
596
597 for (i=0; i < n; i++) {
598
599 _mpd_mul_words(&hi, &lo, u[i], v);
600 lo = carry + lo;
601 if (lo < carry) hi++;
602
603 _mpd_div_words(&carry, &w[i], hi, lo, b);
604 }
605 w[i] = carry;
606 }
607
608 /*
609 * Knuth, TAOCP Volume 2, 4.3.1, exercise 16:
610 * w := quotient of u (len n) divided by a single word v
611 */
612 mpd_uint_t
613 _mpd_shortdiv_b(mpd_uint_t *w, const mpd_uint_t *u, mpd_size_t n,
614 mpd_uint_t v, mpd_uint_t b)
615 {
616 mpd_uint_t hi, lo;
617 mpd_uint_t rem = 0;
618 mpd_size_t i;
619
620 assert(n > 0);
621
622 for (i=n-1; i != MPD_SIZE_MAX; i--) {
623
624 _mpd_mul_words(&hi, &lo, rem, b);
625 lo = u[i] + lo;
626 if (lo < u[i]) hi++;
627
628 _mpd_div_words(&w[i], &rem, hi, lo, v);
629 }
630
631 return rem;
632 }
633
634
635
LEFTRIGHT

RSS Feeds Recent Issues | This issue
This is Rietveld 894c83f36cb7+