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

Side by Side Diff: Modules/_decimal/libmpdec/mpdecimal.c

Issue 7652: Merge C version of decimal into py3k.
Patch Set: Created 7 years, 8 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:
View unified diff | Download patch
« no previous file with comments | « Modules/_decimal/libmpdec/memory.h ('k') | Modules/_decimal/libmpdec/mpdecimal.h » ('j') | no next file with comments »
Toggle Intra-line Diffs ('i') | Expand Comments ('e') | Collapse Comments ('c') | Show Comments Hide Comments ('s')
OLDNEW
(Empty)
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 <limits.h>
34 #include <math.h>
35 #include "basearith.h"
36 #include "bits.h"
37 #include "convolute.h"
38 #include "crt.h"
39 #include "errno.h"
40 #include "memory.h"
41 #include "typearith.h"
42 #include "umodarith.h"
43
44 #ifdef PPRO
45 #if defined(_MSC_VER)
46 #include <float.h>
47 #pragma fenv_access(on)
48 #elif !defined(__OpenBSD__) && !defined(__NetBSD__)
49 /* C99 */
50 #include <fenv.h>
51 #pragma STDC FENV_ACCESS ON
52 #endif
53 #endif
54
55 #if defined(__x86_64__) && defined(__GLIBC__) && !defined(__INTEL_COMPILER)
56 #define USE_80BIT_LONG_DOUBLE
57 #endif
58
59 #if defined(_MSC_VER)
60 #define ALWAYS_INLINE __forceinline
61 #elif defined(LEGACY_COMPILER)
62 #define ALWAYS_INLINE
63 #undef inline
64 #define inline
65 #else
66 #ifdef TEST_COVERAGE
67 #define ALWAYS_INLINE
68 #else
69 #define ALWAYS_INLINE inline __attribute__ ((always_inline))
70 #endif
71 #endif
72
73
74 #define MPD_NEWTONDIV_CUTOFF 1024L
75
76 #define MPD_NEW_STATIC(name, flags, exp, digits, len) \
77 mpd_uint_t name##_data[MPD_MINALLOC_MAX]; \
78 mpd_t name = {flags|MPD_STATIC|MPD_STATIC_DATA, exp, digits, \
79 len, MPD_MINALLOC_MAX, name##_data}
80
81 #define MPD_NEW_CONST(name, flags, exp, digits, len, alloc, initval) \
82 mpd_uint_t name##_data[alloc] = {initval}; \
83 mpd_t name = {flags|MPD_STATIC|MPD_CONST_DATA, exp, digits, \
84 len, alloc, name##_data}
85
86 #define MPD_NEW_SHARED(name, a) \
87 mpd_t name = {(a->flags&~MPD_DATAFLAGS)|MPD_STATIC|MPD_SHARED_DATA, \
88 a->exp, a->digits, a->len, a->alloc, a->data}
89
90
91 static mpd_uint_t data_one[1] = {1};
92 static mpd_uint_t data_zero[1] = {0};
93 static const mpd_t one = {MPD_STATIC|MPD_CONST_DATA, 0, 1, 1, 1, data_one};
94 static const mpd_t minus_one = {MPD_NEG|MPD_STATIC|MPD_CONST_DATA, 0, 1, 1, 1,
95 data_one};
96 static const mpd_t zero = {MPD_STATIC|MPD_CONST_DATA, 0, 1, 1, 1, data_zero};
97
98 static inline void _mpd_check_exp(mpd_t *dec, const mpd_context_t *ctx,
99 uint32_t *status);
100 static void _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a,
101 mpd_ssize_t exp);
102 static inline mpd_ssize_t _mpd_real_size(mpd_uint_t *data, mpd_ssize_t size);
103
104 static void _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
105 const mpd_context_t *ctx, uint32_t *status);
106 static inline void _mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b,
107 const mpd_context_t *ctx, uint32_t *status);
108 static void _mpd_qbarrett_divmod(mpd_t *q, mpd_t *r, const mpd_t *a,
109 const mpd_t *b, uint32_t *status);
110 static inline void _mpd_qpow_uint(mpd_t *result, mpd_t *base, mpd_uint_t exp,
111 uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status);
112
113 mpd_uint_t mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n);
114
115
116 /******************************************************************************/
117 /* Performance critical inline functions */
118 /******************************************************************************/
119
120 #ifdef CONFIG_64
121 /* Digits in a word, primarily useful for the most significant word. */
122 ALWAYS_INLINE int
123 mpd_word_digits(mpd_uint_t word)
124 {
125 if (word < mpd_pow10[9]) {
126 if (word < mpd_pow10[4]) {
127 if (word < mpd_pow10[2]) {
128 return (word < mpd_pow10[1]) ? 1 : 2;
129 }
130 return (word < mpd_pow10[3]) ? 3 : 4;
131 }
132 if (word < mpd_pow10[6]) {
133 return (word < mpd_pow10[5]) ? 5 : 6;
134 }
135 if (word < mpd_pow10[8]) {
136 return (word < mpd_pow10[7]) ? 7 : 8;
137 }
138 return 9;
139 }
140 if (word < mpd_pow10[14]) {
141 if (word < mpd_pow10[11]) {
142 return (word < mpd_pow10[10]) ? 10 : 11;
143 }
144 if (word < mpd_pow10[13]) {
145 return (word < mpd_pow10[12]) ? 12 : 13;
146 }
147 return 14;
148 }
149 if (word < mpd_pow10[18]) {
150 if (word < mpd_pow10[16]) {
151 return (word < mpd_pow10[15]) ? 15 : 16;
152 }
153 return (word < mpd_pow10[17]) ? 17 : 18;
154 }
155
156 return (word < mpd_pow10[19]) ? 19 : 20;
157 }
158 #else
159 ALWAYS_INLINE int
160 mpd_word_digits(mpd_uint_t word)
161 {
162 if (word < mpd_pow10[4]) {
163 if (word < mpd_pow10[2]) {
164 return (word < mpd_pow10[1]) ? 1 : 2;
165 }
166 return (word < mpd_pow10[3]) ? 3 : 4;
167 }
168 if (word < mpd_pow10[6]) {
169 return (word < mpd_pow10[5]) ? 5 : 6;
170 }
171 if (word < mpd_pow10[8]) {
172 return (word < mpd_pow10[7]) ? 7 : 8;
173 }
174
175 return (word < mpd_pow10[9]) ? 9 : 10;
176 }
177 #endif
178
179
180 /* Adjusted exponent */
181 ALWAYS_INLINE mpd_ssize_t
182 mpd_adjexp(const mpd_t *dec)
183 {
184 return (dec->exp + dec->digits) - 1;
185 }
186
187 /* Etiny */
188 ALWAYS_INLINE mpd_ssize_t
189 mpd_etiny(const mpd_context_t *ctx)
190 {
191 return ctx->emin - (ctx->prec - 1);
192 }
193
194 /* Etop: used for folding down in IEEE clamping */
195 ALWAYS_INLINE mpd_ssize_t
196 mpd_etop(const mpd_context_t *ctx)
197 {
198 return ctx->emax - (ctx->prec - 1);
199 }
200
201 /* Most significant word */
202 ALWAYS_INLINE mpd_uint_t
203 mpd_msword(const mpd_t *dec)
204 {
205 assert(dec->len > 0);
206 return dec->data[dec->len-1];
207 }
208
209 /* Most significant digit of a word */
210 inline mpd_uint_t
211 mpd_msd(mpd_uint_t word)
212 {
213 int n;
214
215 n = mpd_word_digits(word);
216 return word / mpd_pow10[n-1];
217 }
218
219 /* Least significant digit of a word */
220 ALWAYS_INLINE mpd_uint_t
221 mpd_lsd(mpd_uint_t word)
222 {
223 return word % 10;
224 }
225
226 /* Coefficient size needed to store 'digits' */
227 ALWAYS_INLINE mpd_ssize_t
228 mpd_digits_to_size(mpd_ssize_t digits)
229 {
230 mpd_ssize_t q, r;
231
232 _mpd_idiv_word(&q, &r, digits, MPD_RDIGITS);
233 return (r == 0) ? q : q+1;
234 }
235
236 /* Number of digits in the exponent. Not defined for MPD_SSIZE_MIN. */
237 inline int
238 mpd_exp_digits(mpd_ssize_t exp)
239 {
240 exp = (exp < 0) ? -exp : exp;
241 return mpd_word_digits(exp);
242 }
243
244 /* Canonical */
245 ALWAYS_INLINE int
246 mpd_iscanonical(const mpd_t *dec UNUSED)
247 {
248 return 1;
249 }
250
251 /* Finite */
252 ALWAYS_INLINE int
253 mpd_isfinite(const mpd_t *dec)
254 {
255 return !(dec->flags & MPD_SPECIAL);
256 }
257
258 /* Infinite */
259 ALWAYS_INLINE int
260 mpd_isinfinite(const mpd_t *dec)
261 {
262 return dec->flags & MPD_INF;
263 }
264
265 /* NaN */
266 ALWAYS_INLINE int
267 mpd_isnan(const mpd_t *dec)
268 {
269 return dec->flags & (MPD_NAN|MPD_SNAN);
270 }
271
272 /* Negative */
273 ALWAYS_INLINE int
274 mpd_isnegative(const mpd_t *dec)
275 {
276 return dec->flags & MPD_NEG;
277 }
278
279 /* Positive */
280 ALWAYS_INLINE int
281 mpd_ispositive(const mpd_t *dec)
282 {
283 return !(dec->flags & MPD_NEG);
284 }
285
286 /* qNaN */
287 ALWAYS_INLINE int
288 mpd_isqnan(const mpd_t *dec)
289 {
290 return dec->flags & MPD_NAN;
291 }
292
293 /* Signed */
294 ALWAYS_INLINE int
295 mpd_issigned(const mpd_t *dec)
296 {
297 return dec->flags & MPD_NEG;
298 }
299
300 /* sNaN */
301 ALWAYS_INLINE int
302 mpd_issnan(const mpd_t *dec)
303 {
304 return dec->flags & MPD_SNAN;
305 }
306
307 /* Special */
308 ALWAYS_INLINE int
309 mpd_isspecial(const mpd_t *dec)
310 {
311 return dec->flags & MPD_SPECIAL;
312 }
313
314 /* Zero */
315 ALWAYS_INLINE int
316 mpd_iszero(const mpd_t *dec)
317 {
318 return !mpd_isspecial(dec) && mpd_msword(dec) == 0;
319 }
320
321 /* Test for zero when specials have been ruled out already */
322 ALWAYS_INLINE int
323 mpd_iszerocoeff(const mpd_t *dec)
324 {
325 return mpd_msword(dec) == 0;
326 }
327
328 /* Normal */
329 inline int
330 mpd_isnormal(const mpd_t *dec, const mpd_context_t *ctx)
331 {
332 if (mpd_isspecial(dec)) return 0;
333 if (mpd_iszerocoeff(dec)) return 0;
334
335 return mpd_adjexp(dec) >= ctx->emin;
336 }
337
338 /* Subnormal */
339 inline int
340 mpd_issubnormal(const mpd_t *dec, const mpd_context_t *ctx)
341 {
342 if (mpd_isspecial(dec)) return 0;
343 if (mpd_iszerocoeff(dec)) return 0;
344
345 return mpd_adjexp(dec) < ctx->emin;
346 }
347
348 /* Odd word */
349 ALWAYS_INLINE int
350 mpd_isoddword(mpd_uint_t word)
351 {
352 return word & 1;
353 }
354
355 /* Odd coefficient */
356 ALWAYS_INLINE int
357 mpd_isoddcoeff(const mpd_t *dec)
358 {
359 return mpd_isoddword(dec->data[0]);
360 }
361
362 /* 0 if dec is positive, 1 if dec is negative */
363 ALWAYS_INLINE uint8_t
364 mpd_sign(const mpd_t *dec)
365 {
366 return dec->flags & MPD_NEG;
367 }
368
369 /* 1 if dec is positive, -1 if dec is negative */
370 ALWAYS_INLINE int
371 mpd_arith_sign(const mpd_t *dec)
372 {
373 return 1 - 2 * mpd_isnegative(dec);
374 }
375
376 /* Radix */
377 ALWAYS_INLINE long
378 mpd_radix(void)
379 {
380 return 10;
381 }
382
383 /* Dynamic decimal */
384 ALWAYS_INLINE int
385 mpd_isdynamic(mpd_t *dec)
386 {
387 return !(dec->flags & MPD_STATIC);
388 }
389
390 /* Static decimal */
391 ALWAYS_INLINE int
392 mpd_isstatic(mpd_t *dec)
393 {
394 return dec->flags & MPD_STATIC;
395 }
396
397 /* Data of decimal is dynamic */
398 ALWAYS_INLINE int
399 mpd_isdynamic_data(mpd_t *dec)
400 {
401 return !(dec->flags & MPD_DATAFLAGS);
402 }
403
404 /* Data of decimal is static */
405 ALWAYS_INLINE int
406 mpd_isstatic_data(mpd_t *dec)
407 {
408 return dec->flags & MPD_STATIC_DATA;
409 }
410
411 /* Data of decimal is shared */
412 ALWAYS_INLINE int
413 mpd_isshared_data(mpd_t *dec)
414 {
415 return dec->flags & MPD_SHARED_DATA;
416 }
417
418 /* Data of decimal is const */
419 ALWAYS_INLINE int
420 mpd_isconst_data(mpd_t *dec)
421 {
422 return dec->flags & MPD_CONST_DATA;
423 }
424
425
426 /******************************************************************************/
427 /* Inline memory handling */
428 /******************************************************************************/
429
430 /* Fill destination with zeros */
431 ALWAYS_INLINE void
432 mpd_uint_zero(mpd_uint_t *dest, mpd_size_t len)
433 {
434 mpd_size_t i;
435
436 for (i = 0; i < len; i++) {
437 dest[i] = 0;
438 }
439 }
440
441 /* Free a decimal */
442 ALWAYS_INLINE void
443 mpd_del(mpd_t *dec)
444 {
445 if (mpd_isdynamic_data(dec)) {
446 mpd_free(dec->data);
447 }
448 if (mpd_isdynamic(dec)) {
449 mpd_free(dec);
450 }
451 }
452
453 /*
454 * Update the memory size for the coefficient. Existing data up to size is
455 * left untouched.
456 *
457 * Error handling: When relloc fails, result->data will still be a valid pointer
458 * to the old memory area of size result->len. If the requested size is less tha n
459 * result->len, we can continue normally, so we treat the failure as a soft erro r.
460 * If the requested size is greater than the old area, MPD_Malloc_error is
461 * set and the result will be a NaN.
462 */
463 ALWAYS_INLINE int
464 mpd_qresize(mpd_t *result, mpd_ssize_t size, uint32_t *status)
465 {
466 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
467 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
468
469 if (mpd_isstatic_data(result)) {
470 if (size > result->alloc) {
471 return mpd_switch_to_dyn(result, size, status);
472 }
473 }
474 else if (size != result->alloc && size >= MPD_MINALLOC) {
475 return mpd_realloc_dyn(result, size, status);
476 }
477
478 return 1;
479 }
480
481 /* Same as mpd_qresize, but the complete coefficient (including the old
482 * memory area!) is initialized to zero. */
483 ALWAYS_INLINE int
484 mpd_qresize_zero(mpd_t *result, mpd_ssize_t size, uint32_t *status)
485 {
486 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
487 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
488
489 if (mpd_isstatic_data(result)) {
490 if (size > result->alloc) {
491 return mpd_switch_to_dyn_zero(result, size, status);
492 }
493 }
494 else if (size != result->alloc && size >= MPD_MINALLOC) {
495 if (!mpd_realloc_dyn(result, size, status)) {
496 return 0;
497 }
498 }
499
500 mpd_uint_zero(result->data, size);
501
502 return 1;
503 }
504
505 /*
506 * Reduce memory size for the coefficient to MPD_MINALLOC. In theory,
507 * realloc may fail even when reducing the memory size. But in that case
508 * the old memory area is always big enough, so checking for MPD_Malloc_error
509 * is not imperative.
510 */
511 ALWAYS_INLINE void
512 mpd_minalloc(mpd_t *result)
513 {
514 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
515 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
516
517 if (!mpd_isstatic_data(result) && result->alloc > MPD_MINALLOC) {
518 uint8_t err = 0;
519 result->data = mpd_realloc(result->data, MPD_MINALLOC,
520 sizeof *result->data, &err);
521 if (!err) {
522 result->alloc = MPD_MINALLOC;
523 }
524 }
525 }
526
527 int
528 mpd_resize(mpd_t *result, mpd_ssize_t size, mpd_context_t *ctx)
529 {
530 uint32_t status = 0;
531 if (!mpd_qresize(result, size, &status)) {
532 mpd_addstatus_raise(ctx, status);
533 return 0;
534 }
535 return 1;
536 }
537
538 int
539 mpd_resize_zero(mpd_t *result, mpd_ssize_t size, mpd_context_t *ctx)
540 {
541 uint32_t status = 0;
542 if (!mpd_qresize_zero(result, size, &status)) {
543 mpd_addstatus_raise(ctx, status);
544 return 0;
545 }
546 return 1;
547 }
548
549
550 /******************************************************************************/
551 /* Set attributes of a decimal */
552 /******************************************************************************/
553
554 /* Set digits. result->len is assumed to be correct. */
555 inline void
556 mpd_setdigits(mpd_t *result)
557 {
558 mpd_ssize_t wdigits = mpd_word_digits(mpd_msword(result));
559 result->digits = wdigits + (result->len-1) * MPD_RDIGITS;
560 }
561
562 /* Set sign */
563 ALWAYS_INLINE void
564 mpd_set_sign(mpd_t *result, uint8_t sign)
565 {
566 result->flags &= ~MPD_NEG;
567 result->flags |= sign;
568 }
569
570 /* Copy sign from another decimal */
571 ALWAYS_INLINE void
572 mpd_signcpy(mpd_t *result, mpd_t *a)
573 {
574 uint8_t sign = a->flags&MPD_NEG;
575
576 result->flags &= ~MPD_NEG;
577 result->flags |= sign;
578 }
579
580 /* Set infinity */
581 ALWAYS_INLINE void
582 mpd_set_infinity(mpd_t *result)
583 {
584 result->flags &= ~MPD_SPECIAL;
585 result->flags |= MPD_INF;
586 }
587
588 /* Set qNaN */
589 ALWAYS_INLINE void
590 mpd_set_qnan(mpd_t *result)
591 {
592 result->flags &= ~MPD_SPECIAL;
593 result->flags |= MPD_NAN;
594 }
595
596 /* Set sNaN */
597 ALWAYS_INLINE void
598 mpd_set_snan(mpd_t *result)
599 {
600 result->flags &= ~MPD_SPECIAL;
601 result->flags |= MPD_SNAN;
602 }
603
604 /* Set to negative */
605 ALWAYS_INLINE void
606 mpd_set_negative(mpd_t *result)
607 {
608 result->flags |= MPD_NEG;
609 }
610
611 /* Set to positive */
612 ALWAYS_INLINE void
613 mpd_set_positive(mpd_t *result)
614 {
615 result->flags &= ~MPD_NEG;
616 }
617
618 /* Set to dynamic */
619 ALWAYS_INLINE void
620 mpd_set_dynamic(mpd_t *result)
621 {
622 result->flags &= ~MPD_STATIC;
623 }
624
625 /* Set to static */
626 ALWAYS_INLINE void
627 mpd_set_static(mpd_t *result)
628 {
629 result->flags |= MPD_STATIC;
630 }
631
632 /* Set data to dynamic */
633 ALWAYS_INLINE void
634 mpd_set_dynamic_data(mpd_t *result)
635 {
636 result->flags &= ~MPD_DATAFLAGS;
637 }
638
639 /* Set data to static */
640 ALWAYS_INLINE void
641 mpd_set_static_data(mpd_t *result)
642 {
643 result->flags &= ~MPD_DATAFLAGS;
644 result->flags |= MPD_STATIC_DATA;
645 }
646
647 /* Set data to shared */
648 ALWAYS_INLINE void
649 mpd_set_shared_data(mpd_t *result)
650 {
651 result->flags &= ~MPD_DATAFLAGS;
652 result->flags |= MPD_SHARED_DATA;
653 }
654
655 /* Set data to const */
656 ALWAYS_INLINE void
657 mpd_set_const_data(mpd_t *result)
658 {
659 result->flags &= ~MPD_DATAFLAGS;
660 result->flags |= MPD_CONST_DATA;
661 }
662
663 /* Clear flags, preserving memory attributes. */
664 ALWAYS_INLINE void
665 mpd_clear_flags(mpd_t *result)
666 {
667 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
668 }
669
670 /* Set flags, preserving memory attributes. */
671 ALWAYS_INLINE void
672 mpd_set_flags(mpd_t *result, uint8_t flags)
673 {
674 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
675 result->flags |= flags;
676 }
677
678 /* Copy flags, preserving memory attributes of result. */
679 ALWAYS_INLINE void
680 mpd_copy_flags(mpd_t *result, const mpd_t *a)
681 {
682 uint8_t aflags = a->flags;
683 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
684 result->flags |= (aflags & ~(MPD_STATIC|MPD_DATAFLAGS));
685 }
686
687 /* Make a work context */
688 static inline void
689 mpd_workcontext(mpd_context_t *workctx, const mpd_context_t *ctx)
690 {
691 workctx->prec = ctx->prec;
692 workctx->emax = ctx->emax;
693 workctx->emin = ctx->emin;
694 workctx->round = ctx->round;
695 workctx->traps = 0;
696 workctx->status= 0;
697 workctx->newtrap= 0;
698 workctx->clamp = ctx->clamp;
699 workctx->allcr = ctx->allcr;
700 }
701
702
703 /******************************************************************************/
704 /* Getting and setting parts of decimals */
705 /******************************************************************************/
706
707 /* Flip the sign of a decimal */
708 static inline void
709 _mpd_negate(mpd_t *dec)
710 {
711 dec->flags ^= MPD_NEG;
712 }
713
714 /* Set coefficient to zero */
715 void
716 mpd_zerocoeff(mpd_t *result)
717 {
718 mpd_minalloc(result);
719 result->digits = 1;
720 result->len = 1;
721 result->data[0] = 0;
722 }
723
724 /* Set the coefficient to all nines. */
725 void
726 mpd_qmaxcoeff(mpd_t *result, const mpd_context_t *ctx, uint32_t *status)
727 {
728 mpd_ssize_t len, r;
729
730 _mpd_idiv_word(&len, &r, ctx->prec, MPD_RDIGITS);
731 len = (r == 0) ? len : len+1;
732
733 if (!mpd_qresize(result, len, status)) {
734 return;
735 }
736
737 result->len = len;
738 result->digits = ctx->prec;
739
740 --len;
741 if (r > 0) {
742 result->data[len--] = mpd_pow10[r]-1;
743 }
744 for (; len >= 0; --len) {
745 result->data[len] = MPD_RADIX-1;
746 }
747 }
748
749 /*
750 * Cut off the most significant digits so that the rest fits in ctx->prec.
751 * Cannot fail.
752 */
753 static void
754 _mpd_cap(mpd_t *result, const mpd_context_t *ctx)
755 {
756 uint32_t dummy;
757 mpd_ssize_t len, r;
758
759 if (result->len > 0 && result->digits > ctx->prec) {
760 _mpd_idiv_word(&len, &r, ctx->prec, MPD_RDIGITS);
761 len = (r == 0) ? len : len+1;
762
763 if (r != 0) {
764 result->data[len-1] %= mpd_pow10[r];
765 }
766
767 len = _mpd_real_size(result->data, len);
768 /* resize to fewer words cannot fail */
769 mpd_qresize(result, len, &dummy);
770 result->len = len;
771 mpd_setdigits(result);
772 }
773 if (mpd_iszero(result)) {
774 _settriple(result, mpd_sign(result), 0, result->exp);
775 }
776 }
777
778 /*
779 * Cut off the most significant digits of a NaN payload so that the rest
780 * fits in ctx->prec - ctx->clamp. Cannot fail.
781 */
782 static void
783 _mpd_fix_nan(mpd_t *result, const mpd_context_t *ctx)
784 {
785 uint32_t dummy;
786 mpd_ssize_t prec;
787 mpd_ssize_t len, r;
788
789 prec = ctx->prec - ctx->clamp;
790 if (result->len > 0 && result->digits > prec) {
791 if (prec == 0) {
792 mpd_minalloc(result);
793 result->len = result->digits = 0;
794 }
795 else {
796 _mpd_idiv_word(&len, &r, prec, MPD_RDIGITS);
797 len = (r == 0) ? len : len+1;
798
799 if (r != 0) {
800 result->data[len-1] %= mpd_pow10[r];
801 }
802
803 len = _mpd_real_size(result->data, len);
804 /* resize to fewer words cannot fail */
805 mpd_qresize(result, len, &dummy);
806 result->len = len;
807 mpd_setdigits(result);
808 if (mpd_iszerocoeff(result)) {
809 /* NaN0 is not a valid representation */
810 result->len = result->digits = 0;
811 }
812 }
813 }
814 }
815
816 /*
817 * Get n most significant digits from a decimal, where 0 < n <= MPD_UINT_DIGITS.
818 * Assumes MPD_UINT_DIGITS == MPD_RDIGITS+1, which is true for 32 and 64 bit
819 * machines.
820 *
821 * The result of the operation will be in lo. If the operation is impossible,
822 * hi will be nonzero. This is used to indicate an error.
823 */
824 static inline void
825 _mpd_get_msdigits(mpd_uint_t *hi, mpd_uint_t *lo, const mpd_t *dec,
826 unsigned int n)
827 {
828 mpd_uint_t r, tmp;
829
830 assert(0 < n && n <= MPD_RDIGITS+1);
831
832 _mpd_div_word(&tmp, &r, dec->digits, MPD_RDIGITS);
833 r = (r == 0) ? MPD_RDIGITS : r; /* digits in the most significant word * /
834
835 *hi = 0;
836 *lo = dec->data[dec->len-1];
837 if (n <= r) {
838 *lo /= mpd_pow10[r-n];
839 }
840 else if (dec->len > 1) {
841 /* at this point 1 <= r < n <= MPD_RDIGITS+1 */
842 _mpd_mul_words(hi, lo, *lo, mpd_pow10[n-r]);
843 tmp = dec->data[dec->len-2] / mpd_pow10[MPD_RDIGITS-(n-r)];
844 *lo = *lo + tmp;
845 if (*lo < tmp) (*hi)++;
846 }
847 }
848
849
850 /******************************************************************************/
851 /* Gathering information about a decimal */
852 /******************************************************************************/
853
854 /* The real size of the coefficient without leading zero words. */
855 static inline mpd_ssize_t
856 _mpd_real_size(mpd_uint_t *data, mpd_ssize_t size)
857 {
858 while (size > 1 && data[size-1] == 0) {
859 size--;
860 }
861
862 return size;
863 }
864
865 /* Return number of trailing zeros. No errors are possible. */
866 mpd_ssize_t
867 mpd_trail_zeros(const mpd_t *dec)
868 {
869 mpd_uint_t word;
870 mpd_ssize_t i, tz = 0;
871
872 for (i=0; i < dec->len; ++i) {
873 if (dec->data[i] != 0) {
874 word = dec->data[i];
875 tz = i * MPD_RDIGITS;
876 while (word % 10 == 0) {
877 word /= 10;
878 tz++;
879 }
880 break;
881 }
882 }
883
884 return tz;
885 }
886
887 /* Integer: Undefined for specials */
888 static int
889 _mpd_isint(const mpd_t *dec)
890 {
891 mpd_ssize_t tz;
892
893 if (mpd_iszerocoeff(dec)) {
894 return 1;
895 }
896
897 tz = mpd_trail_zeros(dec);
898 return (dec->exp + tz >= 0);
899 }
900
901 /* Integer */
902 int
903 mpd_isinteger(const mpd_t *dec)
904 {
905 if (mpd_isspecial(dec)) {
906 return 0;
907 }
908 return _mpd_isint(dec);
909 }
910
911 /* Word is a power of 10 */
912 static int
913 mpd_word_ispow10(mpd_uint_t word)
914 {
915 int n;
916
917 n = mpd_word_digits(word);
918 if (word == mpd_pow10[n-1]) {
919 return 1;
920 }
921
922 return 0;
923 }
924
925 /* Coefficient is a power of 10 */
926 static int
927 mpd_coeff_ispow10(const mpd_t *dec)
928 {
929 if (mpd_word_ispow10(mpd_msword(dec))) {
930 if (_mpd_isallzero(dec->data, dec->len-1)) {
931 return 1;
932 }
933 }
934
935 return 0;
936 }
937
938 /* All digits of a word are nines */
939 static int
940 mpd_word_isallnine(mpd_uint_t word)
941 {
942 int n;
943
944 n = mpd_word_digits(word);
945 if (word == mpd_pow10[n]-1) {
946 return 1;
947 }
948
949 return 0;
950 }
951
952 /* All digits of the coefficient are nines */
953 static int
954 mpd_coeff_isallnine(const mpd_t *dec)
955 {
956 if (mpd_word_isallnine(mpd_msword(dec))) {
957 if (_mpd_isallnine(dec->data, dec->len-1)) {
958 return 1;
959 }
960 }
961
962 return 0;
963 }
964
965 /* Odd decimal: Undefined for non-integers! */
966 int
967 mpd_isodd(const mpd_t *dec)
968 {
969 mpd_uint_t q, r;
970 assert(mpd_isinteger(dec));
971 if (mpd_iszerocoeff(dec)) return 0;
972 if (dec->exp < 0) {
973 _mpd_div_word(&q, &r, -dec->exp, MPD_RDIGITS);
974 q = dec->data[q] / mpd_pow10[r];
975 return mpd_isoddword(q);
976 }
977 return dec->exp == 0 && mpd_isoddword(dec->data[0]);
978 }
979
980 /* Even: Undefined for non-integers! */
981 int
982 mpd_iseven(const mpd_t *dec)
983 {
984 return !mpd_isodd(dec);
985 }
986
987 /******************************************************************************/
988 /* Getting and setting decimals */
989 /******************************************************************************/
990
991 /* Internal function: Set a static decimal from a triple, no error checking. */
992 static void
993 _ssettriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp)
994 {
995 mpd_set_flags(result, sign);
996 result->exp = exp;
997 _mpd_div_word(&result->data[1], &result->data[0], a, MPD_RADIX);
998 result->len = (result->data[1] == 0) ? 1 : 2;
999 mpd_setdigits(result);
1000 }
1001
1002 /* Internal function: Set a decimal from a triple, no error checking. */
1003 static void
1004 _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp)
1005 {
1006 mpd_minalloc(result);
1007 mpd_set_flags(result, sign);
1008 result->exp = exp;
1009 _mpd_div_word(&result->data[1], &result->data[0], a, MPD_RADIX);
1010 result->len = (result->data[1] == 0) ? 1 : 2;
1011 mpd_setdigits(result);
1012 }
1013
1014 /* Set a special number from a triple */
1015 void
1016 mpd_setspecial(mpd_t *result, uint8_t sign, uint8_t type)
1017 {
1018 mpd_minalloc(result);
1019 result->flags &= ~(MPD_NEG|MPD_SPECIAL);
1020 result->flags |= (sign|type);
1021 result->exp = result->digits = result->len = 0;
1022 }
1023
1024 /* Set result of NaN with an error status */
1025 void
1026 mpd_seterror(mpd_t *result, uint32_t flags, uint32_t *status)
1027 {
1028 mpd_minalloc(result);
1029 mpd_set_qnan(result);
1030 mpd_set_positive(result);
1031 result->exp = result->digits = result->len = 0;
1032 *status |= flags;
1033 }
1034
1035 /* quietly set a static decimal from an mpd_ssize_t */
1036 void
1037 mpd_qsset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx,
1038 uint32_t *status)
1039 {
1040 mpd_uint_t u;
1041 uint8_t sign = MPD_POS;
1042
1043 if (a < 0) {
1044 if (a == MPD_SSIZE_MIN) {
1045 u = (mpd_uint_t)MPD_SSIZE_MAX +
1046 (-(MPD_SSIZE_MIN+MPD_SSIZE_MAX));
1047 }
1048 else {
1049 u = -a;
1050 }
1051 sign = MPD_NEG;
1052 }
1053 else {
1054 u = a;
1055 }
1056 _ssettriple(result, sign, u, 0);
1057 mpd_qfinalize(result, ctx, status);
1058 }
1059
1060 /* quietly set a static decimal from an mpd_uint_t */
1061 void
1062 mpd_qsset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx,
1063 uint32_t *status)
1064 {
1065 _ssettriple(result, MPD_POS, a, 0);
1066 mpd_qfinalize(result, ctx, status);
1067 }
1068
1069 /* quietly set a static decimal from an int32_t */
1070 void
1071 mpd_qsset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx,
1072 uint32_t *status)
1073 {
1074 mpd_qsset_ssize(result, a, ctx, status);
1075 }
1076
1077 /* quietly set a static decimal from a uint32_t */
1078 void
1079 mpd_qsset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx,
1080 uint32_t *status)
1081 {
1082 mpd_qsset_uint(result, a, ctx, status);
1083 }
1084
1085 #ifdef CONFIG_64
1086 /* quietly set a static decimal from an int64_t */
1087 void
1088 mpd_qsset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1089 uint32_t *status)
1090 {
1091 mpd_qsset_ssize(result, a, ctx, status);
1092 }
1093
1094 /* quietly set a static decimal from a uint64_t */
1095 void
1096 mpd_qsset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1097 uint32_t *status)
1098 {
1099 mpd_qsset_uint(result, a, ctx, status);
1100 }
1101 #endif
1102
1103 /* quietly set a decimal from an mpd_ssize_t */
1104 void
1105 mpd_qset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx,
1106 uint32_t *status)
1107 {
1108 mpd_minalloc(result);
1109 mpd_qsset_ssize(result, a, ctx, status);
1110 }
1111
1112 /* quietly set a decimal from an mpd_uint_t */
1113 void
1114 mpd_qset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx,
1115 uint32_t *status)
1116 {
1117 _settriple(result, MPD_POS, a, 0);
1118 mpd_qfinalize(result, ctx, status);
1119 }
1120
1121 /* quietly set a decimal from an int32_t */
1122 void
1123 mpd_qset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx,
1124 uint32_t *status)
1125 {
1126 mpd_qset_ssize(result, a, ctx, status);
1127 }
1128
1129 /* quietly set a decimal from a uint32_t */
1130 void
1131 mpd_qset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx,
1132 uint32_t *status)
1133 {
1134 mpd_qset_uint(result, a, ctx, status);
1135 }
1136
1137 #if defined(CONFIG_32) && !defined(LEGACY_COMPILER)
1138 /* set a decimal from a uint64_t */
1139 static void
1140 _c32setu64(mpd_t *result, uint64_t u, uint8_t sign, uint32_t *status)
1141 {
1142 mpd_uint_t w[3];
1143 uint64_t q;
1144 int i, len;
1145
1146 len = 0;
1147 do {
1148 q = u / MPD_RADIX;
1149 w[len] = (mpd_uint_t)(u - q * MPD_RADIX);
1150 u = q; len++;
1151 } while (u != 0);
1152
1153 if (!mpd_qresize(result, len, status)) {
1154 return;
1155 }
1156 for (i = 0; i < len; i++) {
1157 result->data[i] = w[i];
1158 }
1159
1160 mpd_set_sign(result, sign);
1161 result->exp = 0;
1162 result->len = len;
1163 mpd_setdigits(result);
1164 }
1165
1166 static void
1167 _c32_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1168 uint32_t *status)
1169 {
1170 _c32setu64(result, a, MPD_POS, status);
1171 mpd_qfinalize(result, ctx, status);
1172 }
1173
1174 /* set a decimal from an int64_t */
1175 static void
1176 _c32_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1177 uint32_t *status)
1178 {
1179 uint64_t u;
1180 uint8_t sign = MPD_POS;
1181
1182 if (a < 0) {
1183 if (a == INT64_MIN) {
1184 u = (uint64_t)INT64_MAX + (-(INT64_MIN+INT64_MAX));
1185 }
1186 else {
1187 u = -a;
1188 }
1189 sign = MPD_NEG;
1190 }
1191 else {
1192 u = a;
1193 }
1194 _c32setu64(result, u, sign, status);
1195 mpd_qfinalize(result, ctx, status);
1196 }
1197 #endif /* CONFIG_32 && !LEGACY_COMPILER */
1198
1199 #ifndef LEGACY_COMPILER
1200 /* quietly set a decimal from an int64_t */
1201 void
1202 mpd_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1203 uint32_t *status)
1204 {
1205 #ifdef CONFIG_64
1206 mpd_qset_ssize(result, a, ctx, status);
1207 #else
1208 _c32_qset_i64(result, a, ctx, status);
1209 #endif
1210 }
1211
1212 /* quietly set a decimal from a uint64_t */
1213 void
1214 mpd_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1215 uint32_t *status)
1216 {
1217 #ifdef CONFIG_64
1218 mpd_qset_uint(result, a, ctx, status);
1219 #else
1220 _c32_qset_u64(result, a, ctx, status);
1221 #endif
1222 }
1223 #endif /* !LEGACY_COMPILER */
1224
1225
1226 /*
1227 * Quietly get an mpd_uint_t from a decimal. Assumes
1228 * MPD_UINT_DIGITS == MPD_RDIGITS+1, which is true for
1229 * 32 and 64 bit machines.
1230 *
1231 * If the operation is impossible, MPD_Invalid_operation is set.
1232 */
1233 static mpd_uint_t
1234 _mpd_qget_uint(int use_sign, const mpd_t *a, uint32_t *status)
1235 {
1236 mpd_t tmp;
1237 mpd_uint_t tmp_data[2];
1238 mpd_uint_t lo, hi;
1239
1240 if (mpd_isspecial(a)) {
1241 *status |= MPD_Invalid_operation;
1242 return MPD_UINT_MAX;
1243 }
1244 if (mpd_iszero(a)) {
1245 return 0;
1246 }
1247 if (use_sign && mpd_isnegative(a)) {
1248 *status |= MPD_Invalid_operation;
1249 return MPD_UINT_MAX;
1250 }
1251
1252 if (a->digits+a->exp > MPD_RDIGITS+1) {
1253 *status |= MPD_Invalid_operation;
1254 return MPD_UINT_MAX;
1255 }
1256
1257 if (a->exp < 0) {
1258 if (!_mpd_isint(a)) {
1259 *status |= MPD_Invalid_operation;
1260 return MPD_UINT_MAX;
1261 }
1262 /* At this point a->digits+a->exp <= MPD_RDIGITS+1,
1263 * so the shift fits. */
1264 tmp.data = tmp_data;
1265 tmp.flags = MPD_STATIC|MPD_CONST_DATA;
1266 mpd_qsshiftr(&tmp, a, -a->exp);
1267 tmp.exp = 0;
1268 a = &tmp;
1269 }
1270
1271 _mpd_get_msdigits(&hi, &lo, a, MPD_RDIGITS+1);
1272 if (hi) {
1273 *status |= MPD_Invalid_operation;
1274 return MPD_UINT_MAX;
1275 }
1276
1277 if (a->exp > 0) {
1278 _mpd_mul_words(&hi, &lo, lo, mpd_pow10[a->exp]);
1279 if (hi) {
1280 *status |= MPD_Invalid_operation;
1281 return MPD_UINT_MAX;
1282 }
1283 }
1284
1285 return lo;
1286 }
1287
1288 /*
1289 * Sets Invalid_operation for:
1290 * - specials
1291 * - negative numbers (except negative zero)
1292 * - non-integers
1293 * - overflow
1294 */
1295 mpd_uint_t
1296 mpd_qget_uint(const mpd_t *a, uint32_t *status)
1297 {
1298 return _mpd_qget_uint(1, a, status);
1299 }
1300
1301 /* Same as above, but gets the absolute value, i.e. the sign is ignored. */
1302 mpd_uint_t
1303 mpd_qabs_uint(const mpd_t *a, uint32_t *status)
1304 {
1305 return _mpd_qget_uint(0, a, status);
1306 }
1307
1308 /* quietly get an mpd_ssize_t from a decimal */
1309 mpd_ssize_t
1310 mpd_qget_ssize(const mpd_t *a, uint32_t *status)
1311 {
1312 mpd_uint_t u;
1313 int isneg;
1314
1315 u = mpd_qabs_uint(a, status);
1316 if (*status&MPD_Invalid_operation) {
1317 return MPD_SSIZE_MAX;
1318 }
1319
1320 isneg = mpd_isnegative(a);
1321 if (u <= MPD_SSIZE_MAX) {
1322 return isneg ? -((mpd_ssize_t)u) : (mpd_ssize_t)u;
1323 }
1324 else if (isneg && u-1 == MPD_SSIZE_MAX) {
1325 return MPD_SSIZE_MIN;
1326 }
1327
1328 *status |= MPD_Invalid_operation;
1329 return MPD_SSIZE_MAX;
1330 }
1331
1332 #ifdef CONFIG_64
1333 /* quietly get a uint64_t from a decimal */
1334 uint64_t
1335 mpd_qget_u64(const mpd_t *a, uint32_t *status)
1336 {
1337 return mpd_qget_uint(a, status);
1338 }
1339
1340 /* quietly get an int64_t from a decimal */
1341 int64_t
1342 mpd_qget_i64(const mpd_t *a, uint32_t *status)
1343 {
1344 return mpd_qget_ssize(a, status);
1345 }
1346 #else
1347 /* quietly get a uint32_t from a decimal */
1348 uint32_t
1349 mpd_qget_u32(const mpd_t *a, uint32_t *status)
1350 {
1351 return mpd_qget_uint(a, status);
1352 }
1353
1354 /* quietly get an int32_t from a decimal */
1355 int32_t
1356 mpd_qget_i32(const mpd_t *a, uint32_t *status)
1357 {
1358 return mpd_qget_ssize(a, status);
1359 }
1360 #endif
1361
1362
1363 /******************************************************************************/
1364 /* Filtering input of functions, finalizing output of functions */
1365 /******************************************************************************/
1366
1367 /*
1368 * Check if the operand is NaN, copy to result and return 1 if this is
1369 * the case. Copying can fail since NaNs are allowed to have a payload that
1370 * does not fit in MPD_MINALLOC.
1371 */
1372 int
1373 mpd_qcheck_nan(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
1374 uint32_t *status)
1375 {
1376 if (mpd_isnan(a)) {
1377 *status |= mpd_issnan(a) ? MPD_Invalid_operation : 0;
1378 mpd_qcopy(result, a, status);
1379 mpd_set_qnan(result);
1380 _mpd_fix_nan(result, ctx);
1381 return 1;
1382 }
1383 return 0;
1384 }
1385
1386 /*
1387 * Check if either operand is NaN, copy to result and return 1 if this
1388 * is the case. Copying can fail since NaNs are allowed to have a payload
1389 * that does not fit in MPD_MINALLOC.
1390 */
1391 int
1392 mpd_qcheck_nans(mpd_t *result, const mpd_t *a, const mpd_t *b,
1393 const mpd_context_t *ctx, uint32_t *status)
1394 {
1395 if ((a->flags|b->flags)&(MPD_NAN|MPD_SNAN)) {
1396 const mpd_t *choice = b;
1397 if (mpd_issnan(a)) {
1398 choice = a;
1399 *status |= MPD_Invalid_operation;
1400 }
1401 else if (mpd_issnan(b)) {
1402 *status |= MPD_Invalid_operation;
1403 }
1404 else if (mpd_isqnan(a)) {
1405 choice = a;
1406 }
1407 mpd_qcopy(result, choice, status);
1408 mpd_set_qnan(result);
1409 _mpd_fix_nan(result, ctx);
1410 return 1;
1411 }
1412 return 0;
1413 }
1414
1415 /*
1416 * Check if one of the operands is NaN, copy to result and return 1 if this
1417 * is the case. Copying can fail since NaNs are allowed to have a payload
1418 * that does not fit in MPD_MINALLOC.
1419 */
1420 static int
1421 mpd_qcheck_3nans(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c,
1422 const mpd_context_t *ctx, uint32_t *status)
1423 {
1424 if ((a->flags|b->flags|c->flags)&(MPD_NAN|MPD_SNAN)) {
1425 const mpd_t *choice = c;
1426 if (mpd_issnan(a)) {
1427 choice = a;
1428 *status |= MPD_Invalid_operation;
1429 }
1430 else if (mpd_issnan(b)) {
1431 choice = b;
1432 *status |= MPD_Invalid_operation;
1433 }
1434 else if (mpd_issnan(c)) {
1435 *status |= MPD_Invalid_operation;
1436 }
1437 else if (mpd_isqnan(a)) {
1438 choice = a;
1439 }
1440 else if (mpd_isqnan(b)) {
1441 choice = b;
1442 }
1443 mpd_qcopy(result, choice, status);
1444 mpd_set_qnan(result);
1445 _mpd_fix_nan(result, ctx);
1446 return 1;
1447 }
1448 return 0;
1449 }
1450
1451 /* Check if rounding digit 'rnd' leads to an increment. */
1452 static inline int
1453 _mpd_rnd_incr(const mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx)
1454 {
1455 int ld;
1456
1457 switch (ctx->round) {
1458 case MPD_ROUND_DOWN: case MPD_ROUND_TRUNC:
1459 return 0;
1460 case MPD_ROUND_HALF_UP:
1461 return (rnd >= 5);
1462 case MPD_ROUND_HALF_EVEN:
1463 return (rnd > 5) || ((rnd == 5) && mpd_isoddcoeff(dec));
1464 case MPD_ROUND_CEILING:
1465 return !(rnd == 0 || mpd_isnegative(dec));
1466 case MPD_ROUND_FLOOR:
1467 return !(rnd == 0 || mpd_ispositive(dec));
1468 case MPD_ROUND_HALF_DOWN:
1469 return (rnd > 5);
1470 case MPD_ROUND_UP:
1471 return !(rnd == 0);
1472 case MPD_ROUND_05UP:
1473 ld = (int)mpd_lsd(dec->data[0]);
1474 return (!(rnd == 0) && (ld == 0 || ld == 5));
1475 default:
1476 /* Without a valid context, further results will be undefined. * /
1477 return 0; /* GCOV_NOT_REACHED */
1478 }
1479 }
1480
1481 /*
1482 * Apply rounding to a decimal that has been right-shifted into a full
1483 * precision decimal. If an increment leads to an overflow of the precision,
1484 * adjust the coefficient and the exponent and check the new exponent for
1485 * overflow.
1486 */
1487 static inline void
1488 _mpd_apply_round(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx,
1489 uint32_t *status)
1490 {
1491 if (_mpd_rnd_incr(dec, rnd, ctx)) {
1492 /* We have a number with exactly ctx->prec digits. The increment
1493 * can only lead to an overflow if the decimal is all nines. In
1494 * that case, the result is a power of ten with prec+1 digits.
1495 *
1496 * If the precision is a multiple of MPD_RDIGITS, this situation is
1497 * detected by _mpd_baseincr returning a carry.
1498 * If the precision is not a multiple of MPD_RDIGITS, we have to
1499 * check if the result has one digit too many.
1500 */
1501 mpd_uint_t carry = _mpd_baseincr(dec->data, dec->len);
1502 if (carry) {
1503 dec->data[dec->len-1] = mpd_pow10[MPD_RDIGITS-1];
1504 dec->exp += 1;
1505 _mpd_check_exp(dec, ctx, status);
1506 return;
1507 }
1508 mpd_setdigits(dec);
1509 if (dec->digits > ctx->prec) {
1510 mpd_qshiftr_inplace(dec, 1);
1511 dec->exp += 1;
1512 dec->digits = ctx->prec;
1513 _mpd_check_exp(dec, ctx, status);
1514 }
1515 }
1516 }
1517
1518 /*
1519 * Apply rounding to a decimal. Allow overflow of the precision.
1520 */
1521 static inline void
1522 _mpd_apply_round_excess(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx,
1523 uint32_t *status)
1524 {
1525 if (_mpd_rnd_incr(dec, rnd, ctx)) {
1526 mpd_uint_t carry = _mpd_baseincr(dec->data, dec->len);
1527 if (carry) {
1528 if (!mpd_qresize(dec, dec->len+1, status)) {
1529 return;
1530 }
1531 dec->data[dec->len] = 1;
1532 dec->len += 1;
1533 }
1534 mpd_setdigits(dec);
1535 }
1536 }
1537
1538 /*
1539 * Apply rounding to a decimal that has been right-shifted into a decimal
1540 * with full precision or less. Return failure if an increment would
1541 * overflow the precision.
1542 */
1543 static inline int
1544 _mpd_apply_round_fit(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx,
1545 uint32_t *status)
1546 {
1547 if (_mpd_rnd_incr(dec, rnd, ctx)) {
1548 mpd_uint_t carry = _mpd_baseincr(dec->data, dec->len);
1549 if (carry) {
1550 if (!mpd_qresize(dec, dec->len+1, status)) {
1551 return 0;
1552 }
1553 dec->data[dec->len] = 1;
1554 dec->len += 1;
1555 }
1556 mpd_setdigits(dec);
1557 if (dec->digits > ctx->prec) {
1558 mpd_seterror(dec, MPD_Invalid_operation, status);
1559 return 0;
1560 }
1561 }
1562 return 1;
1563 }
1564
1565 /* Check a normal number for overflow, underflow, clamping. If the operand
1566 is modified, it will be zero, special or (sub)normal with a coefficient
1567 that fits into the current context precision. */
1568 static inline void
1569 _mpd_check_exp(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status)
1570 {
1571 mpd_ssize_t adjexp, etiny, shift;
1572 int rnd;
1573
1574 adjexp = mpd_adjexp(dec);
1575 if (adjexp > ctx->emax) {
1576
1577 if (mpd_iszerocoeff(dec)) {
1578 dec->exp = ctx->emax;
1579 if (ctx->clamp) {
1580 dec->exp -= (ctx->prec-1);
1581 }
1582 mpd_zerocoeff(dec);
1583 *status |= MPD_Clamped;
1584 return;
1585 }
1586
1587 switch (ctx->round) {
1588 case MPD_ROUND_HALF_UP: case MPD_ROUND_HALF_EVEN:
1589 case MPD_ROUND_HALF_DOWN: case MPD_ROUND_UP:
1590 case MPD_ROUND_TRUNC:
1591 mpd_setspecial(dec, mpd_sign(dec), MPD_INF);
1592 break;
1593 case MPD_ROUND_DOWN: case MPD_ROUND_05UP:
1594 mpd_qmaxcoeff(dec, ctx, status);
1595 dec->exp = ctx->emax - ctx->prec + 1;
1596 break;
1597 case MPD_ROUND_CEILING:
1598 if (mpd_isnegative(dec)) {
1599 mpd_qmaxcoeff(dec, ctx, status);
1600 dec->exp = ctx->emax - ctx->prec + 1;
1601 }
1602 else {
1603 mpd_setspecial(dec, MPD_POS, MPD_INF);
1604 }
1605 break;
1606 case MPD_ROUND_FLOOR:
1607 if (mpd_ispositive(dec)) {
1608 mpd_qmaxcoeff(dec, ctx, status);
1609 dec->exp = ctx->emax - ctx->prec + 1;
1610 }
1611 else {
1612 mpd_setspecial(dec, MPD_NEG, MPD_INF);
1613 }
1614 break;
1615 default: /* debug */
1616 abort(); /* GCOV_NOT_REACHED */
1617 }
1618
1619 *status |= MPD_Overflow|MPD_Inexact|MPD_Rounded;
1620
1621 } /* fold down */
1622 else if (ctx->clamp && dec->exp > mpd_etop(ctx)) {
1623 /* At this point adjexp=exp+digits-1 <= emax and exp > etop=emax -prec+1:
1624 * (1) shift = exp -emax+prec-1 > 0
1625 * (2) digits+shift = exp+digits-1 - emax + prec <= prec */
1626 shift = dec->exp - mpd_etop(ctx);
1627 if (!mpd_qshiftl(dec, dec, shift, status)) {
1628 return;
1629 }
1630 dec->exp -= shift;
1631 *status |= MPD_Clamped;
1632 if (!mpd_iszerocoeff(dec) && adjexp < ctx->emin) {
1633 /* Underflow is impossible, since exp < etiny=emin-prec+ 1
1634 * and exp > etop=emax-prec+1 would imply emax < emin. * /
1635 *status |= MPD_Subnormal;
1636 }
1637 }
1638 else if (adjexp < ctx->emin) {
1639
1640 etiny = mpd_etiny(ctx);
1641
1642 if (mpd_iszerocoeff(dec)) {
1643 if (dec->exp < etiny) {
1644 dec->exp = etiny;
1645 mpd_zerocoeff(dec);
1646 *status |= MPD_Clamped;
1647 }
1648 return;
1649 }
1650
1651 *status |= MPD_Subnormal;
1652 if (dec->exp < etiny) {
1653 /* At this point adjexp=exp+digits-1 < emin and exp < et iny=emin-prec+1:
1654 * (1) shift = emin-prec+1 - exp > 0
1655 * (2) digits-shift = exp+digits-1 - emin + prec < pre c */
1656 shift = etiny - dec->exp;
1657 rnd = (int)mpd_qshiftr_inplace(dec, shift);
1658 dec->exp = etiny;
1659 /* We always have a spare digit in case of an increment. */
1660 _mpd_apply_round_excess(dec, rnd, ctx, status);
1661 *status |= MPD_Rounded;
1662 if (rnd) {
1663 *status |= (MPD_Inexact|MPD_Underflow);
1664 if (mpd_iszerocoeff(dec)) {
1665 mpd_zerocoeff(dec);
1666 *status |= MPD_Clamped;
1667 }
1668 }
1669 }
1670 /* Case exp >= etiny=emin-prec+1:
1671 * (1) adjexp=exp+digits-1 < emin
1672 * (2) digits < emin-exp+1 <= prec */
1673 }
1674 }
1675
1676 /* Transcendental functions do not always set Underflow reliably,
1677 * since they only use as much precision as is necessary for correct
1678 * rounding. If a result like 1.0000000000e-101 is finalized, there
1679 * is no rounding digit that would trigger Underflow. But we can
1680 * assume Inexact, so a short check suffices. */
1681 static inline void
1682 mpd_check_underflow(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status)
1683 {
1684 if (mpd_adjexp(dec) < ctx->emin && !mpd_iszero(dec) &&
1685 dec->exp < mpd_etiny(ctx)) {
1686 *status |= MPD_Underflow;
1687 }
1688 }
1689
1690 /* Check if a normal number must be rounded after the exponent has been checked. */
1691 static inline void
1692 _mpd_check_round(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status)
1693 {
1694 mpd_uint_t rnd;
1695 mpd_ssize_t shift;
1696
1697 /* must handle specials: _mpd_check_exp() can produce infinities or NaNs */
1698 if (mpd_isspecial(dec)) {
1699 return;
1700 }
1701
1702 if (dec->digits > ctx->prec) {
1703 shift = dec->digits - ctx->prec;
1704 rnd = mpd_qshiftr_inplace(dec, shift);
1705 dec->exp += shift;
1706 _mpd_apply_round(dec, rnd, ctx, status);
1707 *status |= MPD_Rounded;
1708 if (rnd) {
1709 *status |= MPD_Inexact;
1710 }
1711 }
1712 }
1713
1714 /* Finalize all operations. */
1715 void
1716 mpd_qfinalize(mpd_t *result, const mpd_context_t *ctx, uint32_t *status)
1717 {
1718 if (mpd_isspecial(result)) {
1719 if (mpd_isnan(result)) {
1720 _mpd_fix_nan(result, ctx);
1721 }
1722 return;
1723 }
1724
1725 _mpd_check_exp(result, ctx, status);
1726 _mpd_check_round(result, ctx, status);
1727 }
1728
1729
1730 /******************************************************************************/
1731 /* Copying */
1732 /******************************************************************************/
1733
1734 /* Internal function: Copy a decimal, share data with src: USE WITH CARE! */
1735 static inline void
1736 _mpd_copy_shared(mpd_t *dest, const mpd_t *src)
1737 {
1738 dest->flags = src->flags;
1739 dest->exp = src->exp;
1740 dest->digits = src->digits;
1741 dest->len = src->len;
1742 dest->alloc = src->alloc;
1743 dest->data = src->data;
1744
1745 mpd_set_shared_data(dest);
1746 }
1747
1748 /*
1749 * Copy a decimal. In case of an error, status is set to MPD_Malloc_error.
1750 */
1751 int
1752 mpd_qcopy(mpd_t *result, const mpd_t *a, uint32_t *status)
1753 {
1754 if (result == a) return 1;
1755
1756 if (!mpd_qresize(result, a->len, status)) {
1757 return 0;
1758 }
1759
1760 mpd_copy_flags(result, a);
1761 result->exp = a->exp;
1762 result->digits = a->digits;
1763 result->len = a->len;
1764 memcpy(result->data, a->data, a->len * (sizeof *result->data));
1765
1766 return 1;
1767 }
1768
1769 /*
1770 * Copy to a decimal with a static buffer. The caller has to make sure that
1771 * the buffer is big enough. Cannot fail.
1772 */
1773 static void
1774 mpd_qcopy_static(mpd_t *result, const mpd_t *a)
1775 {
1776 if (result == a) return;
1777
1778 memcpy(result->data, a->data, a->len * (sizeof *result->data));
1779
1780 mpd_copy_flags(result, a);
1781 result->exp = a->exp;
1782 result->digits = a->digits;
1783 result->len = a->len;
1784 }
1785
1786 /*
1787 * Return a newly allocated copy of the operand. In case of an error,
1788 * status is set to MPD_Malloc_error and the return value is NULL.
1789 */
1790 mpd_t *
1791 mpd_qncopy(const mpd_t *a)
1792 {
1793 mpd_t *result;
1794
1795 if ((result = mpd_qnew_size(a->len)) == NULL) {
1796 return NULL;
1797 }
1798 memcpy(result->data, a->data, a->len * (sizeof *result->data));
1799 mpd_copy_flags(result, a);
1800 result->exp = a->exp;
1801 result->digits = a->digits;
1802 result->len = a->len;
1803
1804 return result;
1805 }
1806
1807 /*
1808 * Copy a decimal and set the sign to positive. In case of an error, the
1809 * status is set to MPD_Malloc_error.
1810 */
1811 int
1812 mpd_qcopy_abs(mpd_t *result, const mpd_t *a, uint32_t *status)
1813 {
1814 if (!mpd_qcopy(result, a, status)) {
1815 return 0;
1816 }
1817 mpd_set_positive(result);
1818 return 1;
1819 }
1820
1821 /*
1822 * Copy a decimal and negate the sign. In case of an error, the
1823 * status is set to MPD_Malloc_error.
1824 */
1825 int
1826 mpd_qcopy_negate(mpd_t *result, const mpd_t *a, uint32_t *status)
1827 {
1828 if (!mpd_qcopy(result, a, status)) {
1829 return 0;
1830 }
1831 _mpd_negate(result);
1832 return 1;
1833 }
1834
1835 /*
1836 * Copy a decimal, setting the sign of the first operand to the sign of the
1837 * second operand. In case of an error, the status is set to MPD_Malloc_error.
1838 */
1839 int
1840 mpd_qcopy_sign(mpd_t *result, const mpd_t *a, const mpd_t *b, uint32_t *status)
1841 {
1842 uint8_t sign_b = mpd_sign(b); /* result may equal b! */
1843
1844 if (!mpd_qcopy(result, a, status)) {
1845 return 0;
1846 }
1847 mpd_set_sign(result, sign_b);
1848 return 1;
1849 }
1850
1851
1852 /******************************************************************************/
1853 /* Comparisons */
1854 /******************************************************************************/
1855
1856 /*
1857 * For all functions that compare two operands and return an int the usual
1858 * convention applies to the return value:
1859 *
1860 * -1 if op1 < op2
1861 * 0 if op1 == op2
1862 * 1 if op1 > op2
1863 *
1864 * INT_MAX for error
1865 */
1866
1867
1868 /* Convenience macro. If a and b are not equal, return from the calling
1869 * function with the correct comparison value. */
1870 #define CMP_EQUAL_OR_RETURN(a, b) \
1871 if (a != b) { \
1872 if (a < b) { \
1873 return -1; \
1874 } \
1875 return 1; \
1876 }
1877
1878 /*
1879 * Compare the data of big and small. This function does the equivalent
1880 * of first shifting small to the left and then comparing the data of
1881 * big and small, except that no allocation for the left shift is needed.
1882 */
1883 static int
1884 _mpd_basecmp(mpd_uint_t *big, mpd_uint_t *small, mpd_size_t n, mpd_size_t m,
1885 mpd_size_t shift)
1886 {
1887 #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
1888 /* spurious uninitialized warnings */
1889 mpd_uint_t l=l, lprev=lprev, h=h;
1890 #else
1891 mpd_uint_t l, lprev, h;
1892 #endif
1893 mpd_uint_t q, r;
1894 mpd_uint_t ph, x;
1895
1896 assert(m > 0 && n >= m && shift > 0);
1897
1898 _mpd_div_word(&q, &r, (mpd_uint_t)shift, MPD_RDIGITS);
1899
1900 if (r != 0) {
1901
1902 ph = mpd_pow10[r];
1903
1904 --m; --n;
1905 _mpd_divmod_pow10(&h, &lprev, small[m--], MPD_RDIGITS-r);
1906 if (h != 0) {
1907 CMP_EQUAL_OR_RETURN(big[n], h)
1908 --n;
1909 }
1910 for (; m != MPD_SIZE_MAX; m--,n--) {
1911 _mpd_divmod_pow10(&h, &l, small[m], MPD_RDIGITS-r);
1912 x = ph * lprev + h;
1913 CMP_EQUAL_OR_RETURN(big[n], x)
1914 lprev = l;
1915 }
1916 x = ph * lprev;
1917 CMP_EQUAL_OR_RETURN(big[q], x)
1918 }
1919 else {
1920 while (--m != MPD_SIZE_MAX) {
1921 CMP_EQUAL_OR_RETURN(big[m+q], small[m])
1922 }
1923 }
1924
1925 return !_mpd_isallzero(big, q);
1926 }
1927
1928 /* Compare two decimals with the same adjusted exponent. */
1929 static int
1930 _mpd_cmp_same_adjexp(const mpd_t *a, const mpd_t *b)
1931 {
1932 mpd_ssize_t shift, i;
1933
1934 if (a->exp != b->exp) {
1935 /* Cannot wrap: a->exp + a->digits = b->exp + b->digits, so
1936 * a->exp - b->exp = b->digits - a->digits. */
1937 shift = a->exp - b->exp;
1938 if (shift > 0) {
1939 return -1 * _mpd_basecmp(b->data, a->data, b->len, a->le n, shift);
1940 }
1941 else {
1942 return _mpd_basecmp(a->data, b->data, a->len, b->len, -s hift);
1943 }
1944 }
1945
1946 /*
1947 * At this point adjexp(a) == adjexp(b) and a->exp == b->exp,
1948 * so a->digits == b->digits, therefore a->len == b->len.
1949 */
1950 for (i = a->len-1; i >= 0; --i) {
1951 CMP_EQUAL_OR_RETURN(a->data[i], b->data[i])
1952 }
1953
1954 return 0;
1955 }
1956
1957 /* Compare two numerical values. */
1958 static int
1959 _mpd_cmp(const mpd_t *a, const mpd_t *b)
1960 {
1961 mpd_ssize_t adjexp_a, adjexp_b;
1962
1963 /* equal pointers */
1964 if (a == b) {
1965 return 0;
1966 }
1967
1968 /* infinities */
1969 if (mpd_isinfinite(a)) {
1970 if (mpd_isinfinite(b)) {
1971 return mpd_isnegative(b) - mpd_isnegative(a);
1972 }
1973 return mpd_arith_sign(a);
1974 }
1975 if (mpd_isinfinite(b)) {
1976 return -mpd_arith_sign(b);
1977 }
1978
1979 /* zeros */
1980 if (mpd_iszerocoeff(a)) {
1981 if (mpd_iszerocoeff(b)) {
1982 return 0;
1983 }
1984 return -mpd_arith_sign(b);
1985 }
1986 if (mpd_iszerocoeff(b)) {
1987 return mpd_arith_sign(a);
1988 }
1989
1990 /* different signs */
1991 if (mpd_sign(a) != mpd_sign(b)) {
1992 return mpd_sign(b) - mpd_sign(a);
1993 }
1994
1995 /* different adjusted exponents */
1996 adjexp_a = mpd_adjexp(a);
1997 adjexp_b = mpd_adjexp(b);
1998 if (adjexp_a != adjexp_b) {
1999 if (adjexp_a < adjexp_b) {
2000 return -1 * mpd_arith_sign(a);
2001 }
2002 return mpd_arith_sign(a);
2003 }
2004
2005 /* same adjusted exponents */
2006 return _mpd_cmp_same_adjexp(a, b) * mpd_arith_sign(a);
2007 }
2008
2009 /* Compare the absolutes of two numerical values. */
2010 static int
2011 _mpd_cmp_abs(const mpd_t *a, const mpd_t *b)
2012 {
2013 mpd_ssize_t adjexp_a, adjexp_b;
2014
2015 /* equal pointers */
2016 if (a == b) {
2017 return 0;
2018 }
2019
2020 /* infinities */
2021 if (mpd_isinfinite(a)) {
2022 if (mpd_isinfinite(b)) {
2023 return 0;
2024 }
2025 return 1;
2026 }
2027 if (mpd_isinfinite(b)) {
2028 return -1;
2029 }
2030
2031 /* zeros */
2032 if (mpd_iszerocoeff(a)) {
2033 if (mpd_iszerocoeff(b)) {
2034 return 0;
2035 }
2036 return -1;
2037 }
2038 if (mpd_iszerocoeff(b)) {
2039 return 1;
2040 }
2041
2042 /* different adjusted exponents */
2043 adjexp_a = mpd_adjexp(a);
2044 adjexp_b = mpd_adjexp(b);
2045 if (adjexp_a != adjexp_b) {
2046 if (adjexp_a < adjexp_b) {
2047 return -1;
2048 }
2049 return 1;
2050 }
2051
2052 /* same adjusted exponents */
2053 return _mpd_cmp_same_adjexp(a, b);
2054 }
2055
2056 /* Compare two values and return an integer result. */
2057 int
2058 mpd_qcmp(const mpd_t *a, const mpd_t *b, uint32_t *status)
2059 {
2060 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2061 if (mpd_isnan(a) || mpd_isnan(b)) {
2062 *status |= MPD_Invalid_operation;
2063 return INT_MAX;
2064 }
2065 }
2066
2067 return _mpd_cmp(a, b);
2068 }
2069
2070 /*
2071 * Compare a and b, convert the the usual integer result to a decimal and
2072 * store it in 'result'. For convenience, the integer result of the comparison
2073 * is returned. Comparisons involving NaNs return NaN/INT_MAX.
2074 */
2075 int
2076 mpd_qcompare(mpd_t *result, const mpd_t *a, const mpd_t *b,
2077 const mpd_context_t *ctx, uint32_t *status)
2078 {
2079 int c;
2080
2081 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2082 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2083 return INT_MAX;
2084 }
2085 }
2086
2087 c = _mpd_cmp(a, b);
2088 _settriple(result, (c < 0), (c != 0), 0);
2089 return c;
2090 }
2091
2092 /* Same as mpd_compare(), but signal for all NaNs, i.e. also for quiet NaNs. */
2093 int
2094 mpd_qcompare_signal(mpd_t *result, const mpd_t *a, const mpd_t *b,
2095 const mpd_context_t *ctx, uint32_t *status)
2096 {
2097 int c;
2098
2099 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2100 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2101 *status |= MPD_Invalid_operation;
2102 return INT_MAX;
2103 }
2104 }
2105
2106 c = _mpd_cmp(a, b);
2107 _settriple(result, (c < 0), (c != 0), 0);
2108 return c;
2109 }
2110
2111 /* Compare the operands using a total order. */
2112 int
2113 mpd_cmp_total(const mpd_t *a, const mpd_t *b)
2114 {
2115 mpd_t aa, bb;
2116 int nan_a, nan_b;
2117 int c;
2118
2119 if (mpd_sign(a) != mpd_sign(b)) {
2120 return mpd_sign(b) - mpd_sign(a);
2121 }
2122
2123
2124 if (mpd_isnan(a)) {
2125 c = 1;
2126 if (mpd_isnan(b)) {
2127 nan_a = (mpd_isqnan(a)) ? 1 : 0;
2128 nan_b = (mpd_isqnan(b)) ? 1 : 0;
2129 if (nan_b == nan_a) {
2130 if (a->len > 0 && b->len > 0) {
2131 _mpd_copy_shared(&aa, a);
2132 _mpd_copy_shared(&bb, b);
2133 aa.exp = bb.exp = 0;
2134 /* compare payload */
2135 c = _mpd_cmp_abs(&aa, &bb);
2136 }
2137 else {
2138 c = (a->len > 0) - (b->len > 0);
2139 }
2140 }
2141 else {
2142 c = nan_a - nan_b;
2143 }
2144 }
2145 }
2146 else if (mpd_isnan(b)) {
2147 c = -1;
2148 }
2149 else {
2150 c = _mpd_cmp_abs(a, b);
2151 if (c == 0 && a->exp != b->exp) {
2152 c = (a->exp < b->exp) ? -1 : 1;
2153 }
2154 }
2155
2156 return c * mpd_arith_sign(a);
2157 }
2158
2159 /*
2160 * Compare a and b according to a total order, convert the usual integer result
2161 * to a decimal and store it in 'result'. For convenience, the integer result
2162 * of the comparison is returned.
2163 */
2164 int
2165 mpd_compare_total(mpd_t *result, const mpd_t *a, const mpd_t *b)
2166 {
2167 int c;
2168
2169 c = mpd_cmp_total(a, b);
2170 _settriple(result, (c < 0), (c != 0), 0);
2171 return c;
2172 }
2173
2174 /* Compare the magnitude of the operands using a total order. */
2175 int
2176 mpd_cmp_total_mag(const mpd_t *a, const mpd_t *b)
2177 {
2178 mpd_t aa, bb;
2179
2180 _mpd_copy_shared(&aa, a);
2181 _mpd_copy_shared(&bb, b);
2182
2183 mpd_set_positive(&aa);
2184 mpd_set_positive(&bb);
2185
2186 return mpd_cmp_total(&aa, &bb);
2187 }
2188
2189 /*
2190 * Compare the magnitude of a and b according to a total order, convert the
2191 * the usual integer result to a decimal and store it in 'result'.
2192 * For convenience, the integer result of the comparison is returned.
2193 */
2194 int
2195 mpd_compare_total_mag(mpd_t *result, const mpd_t *a, const mpd_t *b)
2196 {
2197 int c;
2198
2199 c = mpd_cmp_total_mag(a, b);
2200 _settriple(result, (c < 0), (c != 0), 0);
2201 return c;
2202 }
2203
2204 /* Determine an ordering for operands that are numerically equal. */
2205 static inline int
2206 _mpd_cmp_numequal(const mpd_t *a, const mpd_t *b)
2207 {
2208 int sign_a, sign_b;
2209 int c;
2210
2211 sign_a = mpd_sign(a);
2212 sign_b = mpd_sign(b);
2213 if (sign_a != sign_b) {
2214 c = sign_b - sign_a;
2215 }
2216 else {
2217 c = (a->exp < b->exp) ? -1 : 1;
2218 c *= mpd_arith_sign(a);
2219 }
2220
2221 return c;
2222 }
2223
2224
2225 /******************************************************************************/
2226 /* Shifting the coefficient */
2227 /******************************************************************************/
2228
2229 /*
2230 * Shift the coefficient of the operand to the left, no check for specials.
2231 * Both operands may be the same pointer. If the result length has to be
2232 * increased, mpd_qresize() might fail with MPD_Malloc_error.
2233 */
2234 int
2235 mpd_qshiftl(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status)
2236 {
2237 mpd_ssize_t size;
2238
2239 assert(n >= 0);
2240
2241 if (mpd_iszerocoeff(a) || n == 0) {
2242 return mpd_qcopy(result, a, status);
2243 }
2244
2245 size = mpd_digits_to_size(a->digits+n);
2246 if (!mpd_qresize(result, size, status)) {
2247 return 0; /* result is NaN */
2248 }
2249
2250 _mpd_baseshiftl(result->data, a->data, size, a->len, n);
2251
2252 mpd_copy_flags(result, a);
2253 result->len = size;
2254 result->exp = a->exp;
2255 result->digits = a->digits+n;
2256
2257 return 1;
2258 }
2259
2260 /* Determine the rounding indicator if all digits of the coefficient are shifted
2261 * out of the picture. */
2262 static mpd_uint_t
2263 _mpd_get_rnd(const mpd_uint_t *data, mpd_ssize_t len, int use_msd)
2264 {
2265 mpd_uint_t rnd = 0, rest = 0, word;
2266
2267 word = data[len-1];
2268 /* special treatment for the most significant digit if shift == digits * /
2269 if (use_msd) {
2270 _mpd_divmod_pow10(&rnd, &rest, word, mpd_word_digits(word)-1);
2271 if (len > 1 && rest == 0) {
2272 rest = !_mpd_isallzero(data, len-1);
2273 }
2274 }
2275 else {
2276 rest = !_mpd_isallzero(data, len);
2277 }
2278
2279 return (rnd == 0 || rnd == 5) ? rnd + !!rest : rnd;
2280 }
2281
2282 /*
2283 * Same as mpd_qshiftr(), but 'result' is a static array. It is the
2284 * caller's responsibility to make sure that the array is big enough.
2285 * The function cannot fail.
2286 */
2287 mpd_uint_t
2288 mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n)
2289 {
2290 mpd_uint_t rnd;
2291 mpd_ssize_t size;
2292
2293 assert(n >= 0);
2294
2295 if (mpd_iszerocoeff(a) || n == 0) {
2296 mpd_qcopy_static(result, a);
2297 return 0;
2298 }
2299
2300 if (n >= a->digits) {
2301 rnd = _mpd_get_rnd(a->data, a->len, (n==a->digits));
2302 mpd_zerocoeff(result);
2303 result->digits = 1;
2304 size = 1;
2305 }
2306 else {
2307 result->digits = a->digits-n;
2308 size = mpd_digits_to_size(result->digits);
2309 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2310 }
2311
2312 mpd_copy_flags(result, a);
2313 result->exp = a->exp;
2314 result->len = size;
2315
2316 return rnd;
2317 }
2318
2319 /*
2320 * Inplace shift of the coefficient to the right, no check for specials.
2321 * Returns the rounding indicator for mpd_rnd_incr().
2322 * The function cannot fail.
2323 */
2324 mpd_uint_t
2325 mpd_qshiftr_inplace(mpd_t *result, mpd_ssize_t n)
2326 {
2327 uint32_t dummy;
2328 mpd_uint_t rnd;
2329 mpd_ssize_t size;
2330
2331 assert(n >= 0);
2332
2333 if (mpd_iszerocoeff(result) || n == 0) {
2334 return 0;
2335 }
2336
2337 if (n >= result->digits) {
2338 rnd = _mpd_get_rnd(result->data, result->len, (n==result->digits ));
2339 mpd_zerocoeff(result);
2340 result->digits = 1;
2341 size = 1;
2342 }
2343 else {
2344 rnd = _mpd_baseshiftr(result->data, result->data, result->len, n );
2345 result->digits -= n;
2346 size = mpd_digits_to_size(result->digits);
2347 /* reducing the size cannot fail */
2348 mpd_qresize(result, size, &dummy);
2349 }
2350
2351 result->len = size;
2352
2353 return rnd;
2354 }
2355
2356 /*
2357 * Shift the coefficient of the operand to the right, no check for specials.
2358 * Both operands may be the same pointer. Returns the rounding indicator to
2359 * be used by mpd_rnd_incr(). If the result length has to be increased,
2360 * mpd_qcopy() or mpd_qresize() might fail with MPD_Malloc_error. In those
2361 * cases, MPD_UINT_MAX is returned.
2362 */
2363 mpd_uint_t
2364 mpd_qshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status)
2365 {
2366 mpd_uint_t rnd;
2367 mpd_ssize_t size;
2368
2369 assert(n >= 0);
2370
2371 if (mpd_iszerocoeff(a) || n == 0) {
2372 if (!mpd_qcopy(result, a, status)) {
2373 return MPD_UINT_MAX;
2374 }
2375 return 0;
2376 }
2377
2378 if (n >= a->digits) {
2379 rnd = _mpd_get_rnd(a->data, a->len, (n==a->digits));
2380 mpd_zerocoeff(result);
2381 result->digits = 1;
2382 size = 1;
2383 }
2384 else {
2385 result->digits = a->digits-n;
2386 size = mpd_digits_to_size(result->digits);
2387 if (result == a) {
2388 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2389 /* reducing the size cannot fail */
2390 mpd_qresize(result, size, status);
2391 }
2392 else {
2393 if (!mpd_qresize(result, size, status)) {
2394 return MPD_UINT_MAX;
2395 }
2396 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2397 }
2398 }
2399
2400 mpd_copy_flags(result, a);
2401 result->exp = a->exp;
2402 result->len = size;
2403
2404 return rnd;
2405 }
2406
2407
2408 /******************************************************************************/
2409 /* Miscellaneous operations */
2410 /******************************************************************************/
2411
2412 /* Logical And */
2413 void
2414 mpd_qand(mpd_t *result, const mpd_t *a, const mpd_t *b,
2415 const mpd_context_t *ctx, uint32_t *status)
2416 {
2417 const mpd_t *big = a, *small = b;
2418 mpd_uint_t x, y, z, xbit, ybit;
2419 int k, mswdigits;
2420 mpd_ssize_t i;
2421
2422 if (mpd_isspecial(a) || mpd_isspecial(b) ||
2423 mpd_isnegative(a) || mpd_isnegative(b) ||
2424 a->exp != 0 || b->exp != 0) {
2425 mpd_seterror(result, MPD_Invalid_operation, status);
2426 return;
2427 }
2428 if (b->digits > a->digits) {
2429 big = b;
2430 small = a;
2431 }
2432 if (!mpd_qresize(result, big->len, status)) {
2433 return;
2434 }
2435
2436
2437 /* full words */
2438 for (i = 0; i < small->len-1; i++) {
2439 x = small->data[i];
2440 y = big->data[i];
2441 z = 0;
2442 for (k = 0; k < MPD_RDIGITS; k++) {
2443 xbit = x % 10;
2444 x /= 10;
2445 ybit = y % 10;
2446 y /= 10;
2447 if (xbit > 1 || ybit > 1) {
2448 goto invalid_operation;
2449 }
2450 z += (xbit&ybit) ? mpd_pow10[k] : 0;
2451 }
2452 result->data[i] = z;
2453 }
2454 /* most significant word of small */
2455 x = small->data[i];
2456 y = big->data[i];
2457 z = 0;
2458 mswdigits = mpd_word_digits(x);
2459 for (k = 0; k < mswdigits; k++) {
2460 xbit = x % 10;
2461 x /= 10;
2462 ybit = y % 10;
2463 y /= 10;
2464 if (xbit > 1 || ybit > 1) {
2465 goto invalid_operation;
2466 }
2467 z += (xbit&ybit) ? mpd_pow10[k] : 0;
2468 }
2469 result->data[i++] = z;
2470
2471 /* scan the rest of y for digit > 1 */
2472 for (; k < MPD_RDIGITS; k++) {
2473 ybit = y % 10;
2474 y /= 10;
2475 if (ybit > 1) {
2476 goto invalid_operation;
2477 }
2478 }
2479 /* scan the rest of big for digit > 1 */
2480 for (; i < big->len; i++) {
2481 y = big->data[i];
2482 for (k = 0; k < MPD_RDIGITS; k++) {
2483 ybit = y % 10;
2484 y /= 10;
2485 if (ybit > 1) {
2486 goto invalid_operation;
2487 }
2488 }
2489 }
2490
2491 mpd_clear_flags(result);
2492 result->exp = 0;
2493 result->len = _mpd_real_size(result->data, small->len);
2494 mpd_qresize(result, result->len, status);
2495 mpd_setdigits(result);
2496 _mpd_cap(result, ctx);
2497 return;
2498
2499 invalid_operation:
2500 mpd_seterror(result, MPD_Invalid_operation, status);
2501 }
2502
2503 /* Class of an operand. Returns a pointer to the constant name. */
2504 const char *
2505 mpd_class(const mpd_t *a, const mpd_context_t *ctx)
2506 {
2507 if (mpd_isnan(a)) {
2508 if (mpd_isqnan(a))
2509 return "NaN";
2510 else
2511 return "sNaN";
2512 }
2513 else if (mpd_ispositive(a)) {
2514 if (mpd_isinfinite(a))
2515 return "+Infinity";
2516 else if (mpd_iszero(a))
2517 return "+Zero";
2518 else if (mpd_isnormal(a, ctx))
2519 return "+Normal";
2520 else
2521 return "+Subnormal";
2522 }
2523 else {
2524 if (mpd_isinfinite(a))
2525 return "-Infinity";
2526 else if (mpd_iszero(a))
2527 return "-Zero";
2528 else if (mpd_isnormal(a, ctx))
2529 return "-Normal";
2530 else
2531 return "-Subnormal";
2532 }
2533 }
2534
2535 /* Logical Xor */
2536 void
2537 mpd_qinvert(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
2538 uint32_t *status)
2539 {
2540 mpd_uint_t x, z, xbit;
2541 mpd_ssize_t i, digits, len;
2542 mpd_ssize_t q, r;
2543 int k;
2544
2545 if (mpd_isspecial(a) || mpd_isnegative(a) || a->exp != 0) {
2546 mpd_seterror(result, MPD_Invalid_operation, status);
2547 return;
2548 }
2549
2550 digits = (a->digits < ctx->prec) ? ctx->prec : a->digits;
2551 _mpd_idiv_word(&q, &r, digits, MPD_RDIGITS);
2552 len = (r == 0) ? q : q+1;
2553 if (!mpd_qresize(result, len, status)) {
2554 return;
2555 }
2556
2557 for (i = 0; i < len; i++) {
2558 x = (i < a->len) ? a->data[i] : 0;
2559 z = 0;
2560 for (k = 0; k < MPD_RDIGITS; k++) {
2561 xbit = x % 10;
2562 x /= 10;
2563 if (xbit > 1) {
2564 goto invalid_operation;
2565 }
2566 z += !xbit ? mpd_pow10[k] : 0;
2567 }
2568 result->data[i] = z;
2569 }
2570
2571 mpd_clear_flags(result);
2572 result->exp = 0;
2573 result->len = _mpd_real_size(result->data, len);
2574 mpd_qresize(result, result->len, status);
2575 mpd_setdigits(result);
2576 _mpd_cap(result, ctx);
2577 return;
2578
2579 invalid_operation:
2580 mpd_seterror(result, MPD_Invalid_operation, status);
2581 }
2582
2583 /* Exponent of the magnitude of the most significant digit of the operand. */
2584 void
2585 mpd_qlogb(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
2586 uint32_t *status)
2587 {
2588 if (mpd_isspecial(a)) {
2589 if (mpd_qcheck_nan(result, a, ctx, status)) {
2590 return;
2591 }
2592 mpd_setspecial(result, MPD_POS, MPD_INF);
2593 }
2594 else if (mpd_iszerocoeff(a)) {
2595 mpd_setspecial(result, MPD_NEG, MPD_INF);
2596 *status |= MPD_Division_by_zero;
2597 }
2598 else {
2599 mpd_qset_ssize(result, mpd_adjexp(a), ctx, status);
2600 }
2601 }
2602
2603 /* Logical Or */
2604 void
2605 mpd_qor(mpd_t *result, const mpd_t *a, const mpd_t *b,
2606 const mpd_context_t *ctx, uint32_t *status)
2607 {
2608 const mpd_t *big = a, *small = b;
2609 mpd_uint_t x, y, z, xbit, ybit;
2610 int k, mswdigits;
2611 mpd_ssize_t i;
2612
2613 if (mpd_isspecial(a) || mpd_isspecial(b) ||
2614 mpd_isnegative(a) || mpd_isnegative(b) ||
2615 a->exp != 0 || b->exp != 0) {
2616 mpd_seterror(result, MPD_Invalid_operation, status);
2617 return;
2618 }
2619 if (b->digits > a->digits) {
2620 big = b;
2621 small = a;
2622 }
2623 if (!mpd_qresize(result, big->len, status)) {
2624 return;
2625 }
2626
2627
2628 /* full words */
2629 for (i = 0; i < small->len-1; i++) {
2630 x = small->data[i];
2631 y = big->data[i];
2632 z = 0;
2633 for (k = 0; k < MPD_RDIGITS; k++) {
2634 xbit = x % 10;
2635 x /= 10;
2636 ybit = y % 10;
2637 y /= 10;
2638 if (xbit > 1 || ybit > 1) {
2639 goto invalid_operation;
2640 }
2641 z += (xbit|ybit) ? mpd_pow10[k] : 0;
2642 }
2643 result->data[i] = z;
2644 }
2645 /* most significant word of small */
2646 x = small->data[i];
2647 y = big->data[i];
2648 z = 0;
2649 mswdigits = mpd_word_digits(x);
2650 for (k = 0; k < mswdigits; k++) {
2651 xbit = x % 10;
2652 x /= 10;
2653 ybit = y % 10;
2654 y /= 10;
2655 if (xbit > 1 || ybit > 1) {
2656 goto invalid_operation;
2657 }
2658 z += (xbit|ybit) ? mpd_pow10[k] : 0;
2659 }
2660
2661 /* scan and copy the rest of y for digit > 1 */
2662 for (; k < MPD_RDIGITS; k++) {
2663 ybit = y % 10;
2664 y /= 10;
2665 if (ybit > 1) {
2666 goto invalid_operation;
2667 }
2668 z += ybit*mpd_pow10[k];
2669 }
2670 result->data[i++] = z;
2671 /* scan and copy the rest of big for digit > 1 */
2672 for (; i < big->len; i++) {
2673 y = big->data[i];
2674 for (k = 0; k < MPD_RDIGITS; k++) {
2675 ybit = y % 10;
2676 y /= 10;
2677 if (ybit > 1) {
2678 goto invalid_operation;
2679 }
2680 }
2681 result->data[i] = big->data[i];
2682 }
2683
2684 mpd_clear_flags(result);
2685 result->exp = 0;
2686 result->len = _mpd_real_size(result->data, big->len);
2687 mpd_qresize(result, result->len, status);
2688 mpd_setdigits(result);
2689 _mpd_cap(result, ctx);
2690 return;
2691
2692 invalid_operation:
2693 mpd_seterror(result, MPD_Invalid_operation, status);
2694 }
2695
2696 /*
2697 * Rotate the coefficient of a by b->data digits. b must be an integer with
2698 * exponent 0.
2699 */
2700 void
2701 mpd_qrotate(mpd_t *result, const mpd_t *a, const mpd_t *b,
2702 const mpd_context_t *ctx, uint32_t *status)
2703 {
2704 uint32_t workstatus = 0;
2705 MPD_NEW_STATIC(tmp,0,0,0,0);
2706 MPD_NEW_STATIC(big,0,0,0,0);
2707 MPD_NEW_STATIC(small,0,0,0,0);
2708 mpd_ssize_t n, lshift, rshift;
2709
2710 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2711 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2712 return;
2713 }
2714 }
2715 if (b->exp != 0 || mpd_isinfinite(b)) {
2716 mpd_seterror(result, MPD_Invalid_operation, status);
2717 return;
2718 }
2719
2720 n = mpd_qget_ssize(b, &workstatus);
2721 if (workstatus&MPD_Invalid_operation) {
2722 mpd_seterror(result, MPD_Invalid_operation, status);
2723 return;
2724 }
2725 if (n > ctx->prec || n < -ctx->prec) {
2726 mpd_seterror(result, MPD_Invalid_operation, status);
2727 return;
2728 }
2729 if (mpd_isinfinite(a)) {
2730 mpd_qcopy(result, a, status);
2731 return;
2732 }
2733
2734 if (n >= 0) {
2735 lshift = n;
2736 rshift = ctx->prec-n;
2737 }
2738 else {
2739 lshift = ctx->prec+n;
2740 rshift = -n;
2741 }
2742
2743 if (a->digits > ctx->prec) {
2744 if (!mpd_qcopy(&tmp, a, status)) {
2745 mpd_seterror(result, MPD_Malloc_error, status);
2746 goto finish;
2747 }
2748 _mpd_cap(&tmp, ctx);
2749 a = &tmp;
2750 }
2751
2752 if (!mpd_qshiftl(&big, a, lshift, status)) {
2753 mpd_seterror(result, MPD_Malloc_error, status);
2754 goto finish;
2755 }
2756 _mpd_cap(&big, ctx);
2757
2758 if (mpd_qshiftr(&small, a, rshift, status) == MPD_UINT_MAX) {
2759 mpd_seterror(result, MPD_Malloc_error, status);
2760 goto finish;
2761 }
2762 _mpd_qadd(result, &big, &small, ctx, status);
2763
2764
2765 finish:
2766 mpd_del(&tmp);
2767 mpd_del(&big);
2768 mpd_del(&small);
2769 }
2770
2771 /*
2772 * b must be an integer with exponent 0 and in the range +-2*(emax + prec).
2773 * XXX: In my opinion +-(2*emax + prec) would be more sensible.
2774 * The result is a with the value of b added to its exponent.
2775 */
2776 void
2777 mpd_qscaleb(mpd_t *result, const mpd_t *a, const mpd_t *b,
2778 const mpd_context_t *ctx, uint32_t *status)
2779 {
2780 uint32_t workstatus = 0;
2781 mpd_uint_t n, maxjump;
2782 #ifndef LEGACY_COMPILER
2783 int64_t exp;
2784 #else
2785 mpd_uint_t x;
2786 int x_sign, n_sign;
2787 mpd_ssize_t exp;
2788 #endif
2789
2790 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2791 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2792 return;
2793 }
2794 }
2795 if (b->exp != 0 || mpd_isinfinite(b)) {
2796 mpd_seterror(result, MPD_Invalid_operation, status);
2797 return;
2798 }
2799
2800 n = mpd_qabs_uint(b, &workstatus);
2801 /* the spec demands this */
2802 maxjump = 2 * (mpd_uint_t)(ctx->emax + ctx->prec);
2803
2804 if (n > maxjump || workstatus&MPD_Invalid_operation) {
2805 mpd_seterror(result, MPD_Invalid_operation, status);
2806 return;
2807 }
2808 if (mpd_isinfinite(a)) {
2809 mpd_qcopy(result, a, status);
2810 return;
2811 }
2812
2813 #ifndef LEGACY_COMPILER
2814 exp = a->exp + (int64_t)n * mpd_arith_sign(b);
2815 exp = (exp > MPD_EXP_INF) ? MPD_EXP_INF : exp;
2816 exp = (exp < MPD_EXP_CLAMP) ? MPD_EXP_CLAMP : exp;
2817 #else
2818 x = (a->exp < 0) ? -a->exp : a->exp;
2819 x_sign = (a->exp < 0) ? 1 : 0;
2820 n_sign = mpd_isnegative(b) ? 1 : 0;
2821
2822 if (x_sign == n_sign) {
2823 x = x + n;
2824 if (x < n) x = MPD_UINT_MAX;
2825 }
2826 else {
2827 x_sign = (x >= n) ? x_sign : n_sign;
2828 x = (x >= n) ? x - n : n - x;
2829 }
2830 if (!x_sign && x > MPD_EXP_INF) x = MPD_EXP_INF;
2831 if (x_sign && x > -MPD_EXP_CLAMP) x = -MPD_EXP_CLAMP;
2832 exp = x_sign ? -((mpd_ssize_t)x) : (mpd_ssize_t)x;
2833 #endif
2834
2835 mpd_qcopy(result, a, status);
2836 result->exp = (mpd_ssize_t)exp;
2837
2838 mpd_qfinalize(result, ctx, status);
2839 }
2840
2841 /*
2842 * Shift the coefficient by n digits, positive n is a left shift. In the case
2843 * of a left shift, the result is decapitated to fit the context precision. If
2844 * you don't want that, use mpd_shiftl().
2845 */
2846 void
2847 mpd_qshiftn(mpd_t *result, const mpd_t *a, mpd_ssize_t n, const mpd_context_t *c tx,
2848 uint32_t *status)
2849 {
2850 if (mpd_isspecial(a)) {
2851 if (mpd_qcheck_nan(result, a, ctx, status)) {
2852 return;
2853 }
2854 mpd_qcopy(result, a, status);
2855 return;
2856 }
2857
2858 if (n >= 0 && n <= ctx->prec) {
2859 mpd_qshiftl(result, a, n, status);
2860 _mpd_cap(result, ctx);
2861 }
2862 else if (n < 0 && n >= -ctx->prec) {
2863 if (!mpd_qcopy(result, a, status)) {
2864 return;
2865 }
2866 _mpd_cap(result, ctx);
2867 mpd_qshiftr_inplace(result, -n);
2868 }
2869 else {
2870 mpd_seterror(result, MPD_Invalid_operation, status);
2871 }
2872 }
2873
2874 /*
2875 * Same as mpd_shiftn(), but the shift is specified by the decimal b, which
2876 * must be an integer with a zero exponent. Infinities remain infinities.
2877 */
2878 void
2879 mpd_qshift(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *c tx,
2880 uint32_t *status)
2881 {
2882 uint32_t workstatus = 0;
2883 mpd_ssize_t n;
2884
2885 if (mpd_isspecial(a) || mpd_isspecial(b)) {
2886 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2887 return;
2888 }
2889 }
2890 if (b->exp != 0 || mpd_isinfinite(b)) {
2891 mpd_seterror(result, MPD_Invalid_operation, status);
2892 return;
2893 }
2894
2895 n = mpd_qget_ssize(b, &workstatus);
2896 if (workstatus&MPD_Invalid_operation) {
2897 mpd_seterror(result, MPD_Invalid_operation, status);
2898 return;
2899 }
2900 if (n > ctx->prec || n < -ctx->prec) {
2901 mpd_seterror(result, MPD_Invalid_operation, status);
2902 return;
2903 }
2904 if (mpd_isinfinite(a)) {
2905 mpd_qcopy(result, a, status);
2906 return;
2907 }
2908
2909 if (n >= 0) {
2910 mpd_qshiftl(result, a, n, status);
2911 _mpd_cap(result, ctx);
2912 }
2913 else {
2914 if (!mpd_qcopy(result, a, status)) {
2915 return;
2916 }
2917 _mpd_cap(result, ctx);
2918 mpd_qshiftr_inplace(result, -n);
2919 }
2920 }
2921
2922 /* Logical Xor */
2923 void
2924 mpd_qxor(mpd_t *result, const mpd_t *a, const mpd_t *b,
2925 const mpd_context_t *ctx, uint32_t *status)
2926 {
2927 const mpd_t *big = a, *small = b;
2928 mpd_uint_t x, y, z, xbit, ybit;
2929 int k, mswdigits;
2930 mpd_ssize_t i;
2931
2932 if (mpd_isspecial(a) || mpd_isspecial(b) ||
2933 mpd_isnegative(a) || mpd_isnegative(b) ||
2934 a->exp != 0 || b->exp != 0) {
2935 mpd_seterror(result, MPD_Invalid_operation, status);
2936 return;
2937 }
2938 if (b->digits > a->digits) {
2939 big = b;
2940 small = a;
2941 }
2942 if (!mpd_qresize(result, big->len, status)) {
2943 return;
2944 }
2945
2946
2947 /* full words */
2948 for (i = 0; i < small->len-1; i++) {
2949 x = small->data[i];
2950 y = big->data[i];
2951 z = 0;
2952 for (k = 0; k < MPD_RDIGITS; k++) {
2953 xbit = x % 10;
2954 x /= 10;
2955 ybit = y % 10;
2956 y /= 10;
2957 if (xbit > 1 || ybit > 1) {
2958 goto invalid_operation;
2959 }
2960 z += (xbit^ybit) ? mpd_pow10[k] : 0;
2961 }
2962 result->data[i] = z;
2963 }
2964 /* most significant word of small */
2965 x = small->data[i];
2966 y = big->data[i];
2967 z = 0;
2968 mswdigits = mpd_word_digits(x);
2969 for (k = 0; k < mswdigits; k++) {
2970 xbit = x % 10;
2971 x /= 10;
2972 ybit = y % 10;
2973 y /= 10;
2974 if (xbit > 1 || ybit > 1) {
2975 goto invalid_operation;
2976 }
2977 z += (xbit^ybit) ? mpd_pow10[k] : 0;
2978 }
2979
2980 /* scan and copy the rest of y for digit > 1 */
2981 for (; k < MPD_RDIGITS; k++) {
2982 ybit = y % 10;
2983 y /= 10;
2984 if (ybit > 1) {
2985 goto invalid_operation;
2986 }
2987 z += ybit*mpd_pow10[k];
2988 }
2989 result->data[i++] = z;
2990 /* scan and copy the rest of big for digit > 1 */
2991 for (; i < big->len; i++) {
2992 y = big->data[i];
2993 for (k = 0; k < MPD_RDIGITS; k++) {
2994 ybit = y % 10;
2995 y /= 10;
2996 if (ybit > 1) {
2997 goto invalid_operation;
2998 }
2999 }
3000 result->data[i] = big->data[i];
3001 }
3002
3003 mpd_clear_flags(result);
3004 result->exp = 0;
3005 result->len = _mpd_real_size(result->data, big->len);
3006 mpd_qresize(result, result->len, status);
3007 mpd_setdigits(result);
3008 _mpd_cap(result, ctx);
3009 return;
3010
3011 invalid_operation:
3012 mpd_seterror(result, MPD_Invalid_operation, status);
3013 }
3014
3015
3016 /******************************************************************************/
3017 /* Arithmetic operations */
3018 /******************************************************************************/
3019
3020 /*
3021 * The absolute value of a. If a is negative, the result is the same
3022 * as the result of the minus operation. Otherwise, the result is the
3023 * result of the plus operation.
3024 */
3025 void
3026 mpd_qabs(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
3027 uint32_t *status)
3028 {
3029 if (mpd_isspecial(a)) {
3030 if (mpd_qcheck_nan(result, a, ctx, status)) {
3031 return;
3032 }
3033 }
3034
3035 if (mpd_isnegative(a)) {
3036 mpd_qminus(result, a, ctx, status);
3037 }
3038 else {
3039 mpd_qplus(result, a, ctx, status);
3040 }
3041
3042 mpd_qfinalize(result, ctx, status);
3043 }
3044
3045 static inline void
3046 _mpd_ptrswap(mpd_t **a, mpd_t **b)
3047 {
3048 mpd_t *t = *a;
3049 *a = *b;
3050 *b = t;
3051 }
3052
3053 /* Add or subtract infinities. */
3054 static void
3055 _mpd_qaddsub_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
3056 uint32_t *status)
3057 {
3058 if (mpd_isinfinite(a)) {
3059 if (mpd_sign(a) != sign_b && mpd_isinfinite(b)) {
3060 mpd_seterror(result, MPD_Invalid_operation, status);
3061 }
3062 else {
3063 mpd_setspecial(result, mpd_sign(a), MPD_INF);
3064 }
3065 return;
3066 }
3067 assert(mpd_isinfinite(b));
3068 mpd_setspecial(result, sign_b, MPD_INF);
3069 }
3070
3071 /* Add or subtract non-special numbers. */
3072 static void
3073 _mpd_qaddsub(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
3074 const mpd_context_t *ctx, uint32_t *status)
3075 {
3076 mpd_t *big, *small;
3077 MPD_NEW_STATIC(big_aligned,0,0,0,0);
3078 MPD_NEW_CONST(tiny,0,0,0,1,1,1);
3079 mpd_uint_t carry;
3080 mpd_ssize_t newsize, shift;
3081 mpd_ssize_t exp, i;
3082 int swap = 0;
3083
3084
3085 /* compare exponents */
3086 big = (mpd_t *)a; small = (mpd_t *)b;
3087 if (big->exp != small->exp) {
3088 if (small->exp > big->exp) {
3089 _mpd_ptrswap(&big, &small);
3090 swap++;
3091 }
3092 if (!mpd_iszerocoeff(big)) {
3093 /* Test for adjexp(small) + big->digits < adjexp(big), i f big-digits > prec
3094 * Test for adjexp(small) + prec + 1 < adjexp(big), i f big-digits <= prec
3095 * If true, the magnitudes of the numbers are so far apa rt that one can as
3096 * well add or subtract 1*10**big->exp. */
3097 exp = big->exp - 1;
3098 exp += (big->digits > ctx->prec) ? 0 : big->digits-ctx-> prec-1;
3099 if (mpd_adjexp(small) < exp) {
3100 mpd_copy_flags(&tiny, small);
3101 tiny.exp = exp;
3102 tiny.digits = 1;
3103 tiny.len = 1;
3104 tiny.data[0] = mpd_iszerocoeff(small) ? 0 : 1;
3105 small = &tiny;
3106 }
3107 /* this cannot wrap: the difference is positive and <= m axprec+1 */
3108 shift = big->exp - small->exp;
3109 if (!mpd_qshiftl(&big_aligned, big, shift, status)) {
3110 mpd_seterror(result, MPD_Malloc_error, status);
3111 goto finish;
3112 }
3113 big = &big_aligned;
3114 }
3115 }
3116 result->exp = small->exp;
3117
3118
3119 /* compare length of coefficients */
3120 if (big->len < small->len) {
3121 _mpd_ptrswap(&big, &small);
3122 swap++;
3123 }
3124
3125 newsize = big->len;
3126 if (!mpd_qresize(result, newsize, status)) {
3127 goto finish;
3128 }
3129
3130 if (mpd_sign(a) == sign_b) {
3131
3132 carry = _mpd_baseadd(result->data, big->data, small->data,
3133 big->len, small->len);
3134
3135 if (carry) {
3136 newsize = big->len + 1;
3137 if (!mpd_qresize(result, newsize, status)) {
3138 goto finish;
3139 }
3140 result->data[newsize-1] = carry;
3141 }
3142
3143 result->len = newsize;
3144 mpd_set_flags(result, sign_b);
3145 }
3146 else {
3147 if (big->len == small->len) {
3148 for (i=big->len-1; i >= 0; --i) {
3149 if (big->data[i] != small->data[i]) {
3150 if (big->data[i] < small->data[i]) {
3151 _mpd_ptrswap(&big, &small);
3152 swap++;
3153 }
3154 break;
3155 }
3156 }
3157 }
3158
3159 _mpd_basesub(result->data, big->data, small->data,
3160 big->len, small->len);
3161 newsize = _mpd_real_size(result->data, big->len);
3162 /* resize to smaller cannot fail */
3163 (void)mpd_qresize(result, newsize, status);
3164
3165 result->len = newsize;
3166 sign_b = (swap & 1) ? sign_b : mpd_sign(a);
3167 mpd_set_flags(result, sign_b);
3168
3169 if (mpd_iszerocoeff(result)) {
3170 mpd_set_positive(result);
3171 if (ctx->round == MPD_ROUND_FLOOR) {
3172 mpd_set_negative(result);
3173 }
3174 }
3175 }
3176
3177 mpd_setdigits(result);
3178
3179 finish:
3180 mpd_del(&big_aligned);
3181 }
3182
3183 /* Add a and b. No specials, no finalizing. */
3184 static void
3185 _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
3186 const mpd_context_t *ctx, uint32_t *status)
3187 {
3188 _mpd_qaddsub(result, a, b, mpd_sign(b), ctx, status);
3189 }
3190
3191 /* Subtract b from a. No specials, no finalizing. */
3192 static void
3193 _mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b,
3194 const mpd_context_t *ctx, uint32_t *status)
3195 {
3196 _mpd_qaddsub(result, a, b, !mpd_sign(b), ctx, status);
3197 }
3198
3199 /* Add a and b. */
3200 void
3201 mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
3202 const mpd_context_t *ctx, uint32_t *status)
3203 {
3204 if (mpd_isspecial(a) || mpd_isspecial(b)) {
3205 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3206 return;
3207 }
3208 _mpd_qaddsub_inf(result, a, b, mpd_sign(b), status);
3209 return;
3210 }
3211
3212 _mpd_qaddsub(result, a, b, mpd_sign(b), ctx, status);
3213 mpd_qfinalize(result, ctx, status);
3214 }
3215
3216 /* Subtract b from a. */
3217 void
3218 mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b,
3219 const mpd_context_t *ctx, uint32_t *status)
3220 {
3221 if (mpd_isspecial(a) || mpd_isspecial(b)) {
3222 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3223 return;
3224 }
3225 _mpd_qaddsub_inf(result, a, b, !mpd_sign(b), status);
3226 return;
3227 }
3228
3229 _mpd_qaddsub(result, a, b, !mpd_sign(b), ctx, status);
3230 mpd_qfinalize(result, ctx, status);
3231 }
3232
3233 /* Add decimal and mpd_ssize_t. */
3234 void
3235 mpd_qadd_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
3236 const mpd_context_t *ctx, uint32_t *status)
3237 {
3238 mpd_context_t maxcontext;
3239 MPD_NEW_STATIC(bb,0,0,0,0);
3240
3241 mpd_maxcontext(&maxcontext);
3242 mpd_qsset_ssize(&bb, b, &maxcontext, status);
3243 mpd_qadd(result, a, &bb, ctx, status);
3244 mpd_del(&bb);
3245 }
3246
3247 /* Add decimal and mpd_uint_t. */
3248 void
3249 mpd_qadd_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
3250 const mpd_context_t *ctx, uint32_t *status)
3251 {
3252 mpd_context_t maxcontext;
3253 MPD_NEW_STATIC(bb,0,0,0,0);
3254
3255 mpd_maxcontext(&maxcontext);
3256 mpd_qsset_uint(&bb, b, &maxcontext, status);
3257 mpd_qadd(result, a, &bb, ctx, status);
3258 mpd_del(&bb);
3259 }
3260
3261 /* Subtract mpd_ssize_t from decimal. */
3262 void
3263 mpd_qsub_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
3264 const mpd_context_t *ctx, uint32_t *status)
3265 {
3266 mpd_context_t maxcontext;
3267 MPD_NEW_STATIC(bb,0,0,0,0);
3268
3269 mpd_maxcontext(&maxcontext);
3270 mpd_qsset_ssize(&bb, b, &maxcontext, status);
3271 mpd_qsub(result, a, &bb, ctx, status);
3272 mpd_del(&bb);
3273 }
3274
3275 /* Subtract mpd_uint_t from decimal. */
3276 void
3277 mpd_qsub_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
3278 const mpd_context_t *ctx, uint32_t *status)
3279 {
3280 mpd_context_t maxcontext;
3281 MPD_NEW_STATIC(bb,0,0,0,0);
3282
3283 mpd_maxcontext(&maxcontext);
3284 mpd_qsset_uint(&bb, b, &maxcontext, status);
3285 mpd_qsub(result, a, &bb, ctx, status);
3286 mpd_del(&bb);
3287 }
3288
3289 /* Add decimal and int32_t. */
3290 void
3291 mpd_qadd_i32(mpd_t *result, const mpd_t *a, int32_t b,
3292 const mpd_context_t *ctx, uint32_t *status)
3293 {
3294 mpd_qadd_ssize(result, a, b, ctx, status);
3295 }
3296
3297 /* Add decimal and uint32_t. */
3298 void
3299 mpd_qadd_u32(mpd_t *result, const mpd_t *a, uint32_t b,
3300 const mpd_context_t *ctx, uint32_t *status)
3301 {
3302 mpd_qadd_uint(result, a, b, ctx, status);
3303 }
3304
3305 #ifdef CONFIG_64
3306 /* Add decimal and int64_t. */
3307 void
3308 mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b,
3309 const mpd_context_t *ctx, uint32_t *status)
3310 {
3311 mpd_qadd_ssize(result, a, b, ctx, status);
3312 }
3313
3314 /* Add decimal and uint64_t. */
3315 void
3316 mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3317 const mpd_context_t *ctx, uint32_t *status)
3318 {
3319 mpd_qadd_uint(result, a, b, ctx, status);
3320 }
3321 #endif
3322
3323 /* Subtract int32_t from decimal. */
3324 void
3325 mpd_qsub_i32(mpd_t *result, const mpd_t *a, int32_t b,
3326 const mpd_context_t *ctx, uint32_t *status)
3327 {
3328 mpd_qsub_ssize(result, a, b, ctx, status);
3329 }
3330
3331 /* Subtract uint32_t from decimal. */
3332 void
3333 mpd_qsub_u32(mpd_t *result, const mpd_t *a, uint32_t b,
3334 const mpd_context_t *ctx, uint32_t *status)
3335 {
3336 mpd_qsub_uint(result, a, b, ctx, status);
3337 }
3338
3339 #ifdef CONFIG_64
3340 /* Subtract int64_t from decimal. */
3341 void
3342 mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b,
3343 const mpd_context_t *ctx, uint32_t *status)
3344 {
3345 mpd_qsub_ssize(result, a, b, ctx, status);
3346 }
3347
3348 /* Subtract uint64_t from decimal. */
3349 void
3350 mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3351 const mpd_context_t *ctx, uint32_t *status)
3352 {
3353 mpd_qsub_uint(result, a, b, ctx, status);
3354 }
3355 #endif
3356
3357
3358 /* Divide infinities. */
3359 static void
3360 _mpd_qdiv_inf(mpd_t *result, const mpd_t *a, const mpd_t *b,
3361 const mpd_context_t *ctx, uint32_t *status)
3362 {
3363 if (mpd_isinfinite(a)) {
3364 if (mpd_isinfinite(b)) {
3365 mpd_seterror(result, MPD_Invalid_operation, status);
3366 return;
3367 }
3368 mpd_setspecial(result, mpd_sign(a)^mpd_sign(b), MPD_INF);
3369 return;
3370 }
3371 assert(mpd_isinfinite(b));
3372 _settriple(result, mpd_sign(a)^mpd_sign(b), 0, mpd_etiny(ctx));
3373 *status |= MPD_Clamped;
3374 }
3375
3376 enum {NO_IDEAL_EXP, SET_IDEAL_EXP};
3377 /* Divide a by b. */
3378 static void
3379 _mpd_qdiv(int action, mpd_t *q, const mpd_t *a, const mpd_t *b,
3380 const mpd_context_t *ctx, uint32_t *status)
3381 {
3382 MPD_NEW_STATIC(aligned,0,0,0,0);
3383 mpd_uint_t ld;
3384 mpd_ssize_t shift, exp, tz;
3385 mpd_ssize_t newsize;
3386 mpd_ssize_t ideal_exp;
3387 mpd_uint_t rem;
3388 uint8_t sign_a = mpd_sign(a);
3389 uint8_t sign_b = mpd_sign(b);
3390
3391
3392 if (mpd_isspecial(a) || mpd_isspecial(b)) {
3393 if (mpd_qcheck_nans(q, a, b, ctx, status)) {
3394 return;
3395 }
3396 _mpd_qdiv_inf(q, a, b, ctx, status);
3397 return;
3398 }
3399 if (mpd_iszerocoeff(b)) {
3400 if (mpd_iszerocoeff(a)) {
3401 mpd_seterror(q, MPD_Division_undefined, status);
3402 }
3403 else {
3404 mpd_setspecial(q, sign_a^sign_b, MPD_INF);
3405 *status |= MPD_Division_by_zero;
3406 }
3407 return;
3408 }
3409 if (mpd_iszerocoeff(a)) {
3410 exp = a->exp - b->exp;
3411 _settriple(q, sign_a^sign_b, 0, exp);
3412 mpd_qfinalize(q, ctx, status);
3413 return;
3414 }
3415
3416 shift = (b->digits - a->digits) + ctx->prec + 1;
3417 ideal_exp = a->exp - b->exp;
3418 exp = ideal_exp - shift;
3419 if (shift > 0) {
3420 if (!mpd_qshiftl(&aligned, a, shift, status)) {
3421 mpd_seterror(q, MPD_Malloc_error, status);
3422 goto finish;
3423 }
3424 a = &aligned;
3425 }
3426 else if (shift < 0) {
3427 shift = -shift;
3428 if (!mpd_qshiftl(&aligned, b, shift, status)) {
3429 mpd_seterror(q, MPD_Malloc_error, status);
3430 goto finish;
3431 }
3432 b = &aligned;
3433 }
3434
3435
3436 newsize = a->len - b->len + 1;
3437 if ((q != b && q != a) || (q == b && newsize > b->len)) {
3438 if (!mpd_qresize(q, newsize, status)) {
3439 mpd_seterror(q, MPD_Malloc_error, status);
3440 goto finish;
3441 }
3442 }
3443
3444
3445 if (b->len == 1) {
3446 rem = _mpd_shortdiv(q->data, a->data, a->len, b->data[0]);
3447 }
3448 else if (a->len < 2*MPD_NEWTONDIV_CUTOFF &&
3449 b->len < MPD_NEWTONDIV_CUTOFF) {
3450 int ret = _mpd_basedivmod(q->data, NULL, a->data, b->data,
3451 a->len, b->len);
3452 if (ret < 0) {
3453 mpd_seterror(q, MPD_Malloc_error, status);
3454 goto finish;
3455 }
3456 rem = ret;
3457 }
3458 else {
3459 MPD_NEW_STATIC(r,0,0,0,0);
3460 _mpd_qbarrett_divmod(q, &r, a, b, status);
3461 if (mpd_isspecial(q) || mpd_isspecial(&r)) {
3462 mpd_del(&r);
3463 goto finish;
3464 }
3465 rem = !mpd_iszerocoeff(&r);
3466 mpd_del(&r);
3467 newsize = q->len;
3468 }
3469
3470 newsize = _mpd_real_size(q->data, newsize);
3471 /* resize to smaller cannot fail */
3472 mpd_qresize(q, newsize, status);
3473 q->len = newsize;
3474 mpd_setdigits(q);
3475
3476 shift = ideal_exp - exp;
3477 if (rem) {
3478 ld = mpd_lsd(q->data[0]);
3479 if (ld == 0 || ld == 5) {
3480 q->data[0] += 1;
3481 }
3482 }
3483 else if (action == SET_IDEAL_EXP && shift > 0) {
3484 tz = mpd_trail_zeros(q);
3485 shift = (tz > shift) ? shift : tz;
3486 mpd_qshiftr_inplace(q, shift);
3487 exp += shift;
3488 }
3489
3490 mpd_set_flags(q, sign_a^sign_b);
3491 q->exp = exp;
3492
3493
3494 finish:
3495 mpd_del(&aligned);
3496 mpd_qfinalize(q, ctx, status);
3497 }
3498
3499 /* Divide a by b. */
3500 void
3501 mpd_qdiv(mpd_t *q, const mpd_t *a, const mpd_t *b,
3502 const mpd_context_t *ctx, uint32_t *status)
3503 {
3504 _mpd_qdiv(SET_IDEAL_EXP, q, a, b, ctx, status);
3505 }
3506
3507 /* Internal function. */
3508 static void
3509 _mpd_qdivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b,
3510 const mpd_context_t *ctx, uint32_t *status)
3511 {
3512 MPD_NEW_STATIC(aligned,0,0,0,0);
3513 mpd_ssize_t qsize, rsize;
3514 mpd_ssize_t ideal_exp, expdiff, shift;
3515 uint8_t sign_a = mpd_sign(a);
3516 uint8_t sign_ab = mpd_sign(a)^mpd_sign(b);
3517
3518
3519 ideal_exp = (a->exp > b->exp) ? b->exp : a->exp;
3520 if (mpd_iszerocoeff(a)) {
3521 if (!mpd_qcopy(r, a, status)) {
3522 goto nanresult; /* GCOV_NOT_REACHED */
3523 }
3524 r->exp = ideal_exp;
3525 _settriple(q, sign_ab, 0, 0);
3526 return;
3527 }
3528
3529 expdiff = mpd_adjexp(a) - mpd_adjexp(b);
3530 if (expdiff < 0) {
3531 if (a->exp > b->exp) {
3532 /* positive and less than b->digits - a->digits */
3533 shift = a->exp - b->exp;
3534 if (!mpd_qshiftl(r, a, shift, status)) {
3535 goto nanresult;
3536 }
3537 r->exp = ideal_exp;
3538 }
3539 else {
3540 if (!mpd_qcopy(r, a, status)) {
3541 goto nanresult;
3542 }
3543 }
3544 _settriple(q, sign_ab, 0, 0);
3545 return;
3546 }
3547 if (expdiff > ctx->prec) {
3548 *status |= MPD_Division_impossible;
3549 goto nanresult;
3550 }
3551
3552
3553 /*
3554 * At this point we have:
3555 * (1) 0 <= a->exp + a->digits - b->exp - b->digits <= prec
3556 * (2) a->exp - b->exp >= b->digits - a->digits
3557 * (3) a->exp - b->exp <= prec + b->digits - a->digits
3558 */
3559 if (a->exp != b->exp) {
3560 shift = a->exp - b->exp;
3561 if (shift > 0) {
3562 /* by (3), after the shift a->digits <= prec + b->digits */
3563 if (!mpd_qshiftl(&aligned, a, shift, status)) {
3564 goto nanresult;
3565 }
3566 a = &aligned;
3567 }
3568 else {
3569 shift = -shift;
3570 /* by (2), after the shift b->digits <= a->digits */
3571 if (!mpd_qshiftl(&aligned, b, shift, status)) {
3572 goto nanresult;
3573 }
3574 b = &aligned;
3575 }
3576 }
3577
3578
3579 qsize = a->len - b->len + 1;
3580 if (!(q == a && qsize < a->len) && !(q == b && qsize < b->len)) {
3581 if (!mpd_qresize(q, qsize, status)) {
3582 goto nanresult;
3583 }
3584 }
3585
3586 rsize = b->len;
3587 if (!(r == a && rsize < a->len)) {
3588 if (!mpd_qresize(r, rsize, status)) {
3589 goto nanresult;
3590 }
3591 }
3592
3593 if (b->len == 1) {
3594 if (a->len == 1) {
3595 _mpd_div_word(&q->data[0], &r->data[0], a->data[0], b->d ata[0]);
3596 }
3597 else {
3598 r->data[0] = _mpd_shortdiv(q->data, a->data, a->len, b-> data[0]);
3599 }
3600 }
3601 else if (a->len < 2*MPD_NEWTONDIV_CUTOFF &&
3602 b->len < MPD_NEWTONDIV_CUTOFF) {
3603 int ret;
3604 ret = _mpd_basedivmod(q->data, r->data, a->data, b->data,
3605 a->len, b->len);
3606 if (ret == -1) {
3607 *status |= MPD_Malloc_error;
3608 goto nanresult;
3609 }
3610 }
3611 else {
3612 _mpd_qbarrett_divmod(q, r, a, b, status);
3613 if (mpd_isspecial(q) || mpd_isspecial(r)) {
3614 goto nanresult;
3615 }
3616 if (mpd_isinfinite(q) || q->digits > ctx->prec) {
3617 *status |= MPD_Division_impossible;
3618 goto nanresult;
3619 }
3620 qsize = q->len;
3621 rsize = r->len;
3622 }
3623
3624 qsize = _mpd_real_size(q->data, qsize);
3625 /* resize to smaller cannot fail */
3626 mpd_qresize(q, qsize, status);
3627 q->len = qsize;
3628 mpd_setdigits(q);
3629 mpd_set_flags(q, sign_ab);
3630 q->exp = 0;
3631 if (q->digits > ctx->prec) {
3632 *status |= MPD_Division_impossible;
3633 goto nanresult;
3634 }
3635
3636 rsize = _mpd_real_size(r->data, rsize);
3637 /* resize to smaller cannot fail */
3638 mpd_qresize(r, rsize, status);
3639 r->len = rsize;
3640 mpd_setdigits(r);
3641 mpd_set_flags(r, sign_a);
3642 r->exp = ideal_exp;
3643
3644 out:
3645 mpd_del(&aligned);
3646 return;
3647
3648 nanresult:
3649 mpd_setspecial(q, MPD_POS, MPD_NAN);
3650 mpd_setspecial(r, MPD_POS, MPD_NAN);
3651 goto out;
3652 }
3653
3654 /* Integer division with remainder. */
3655 void
3656 mpd_qdivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b,
3657 const mpd_context_t *ctx, uint32_t *status)
3658 {
3659 uint8_t sign = mpd_sign(a)^mpd_sign(b);
3660
3661 if (mpd_isspecial(a) || mpd_isspecial(b)) {
3662 if (mpd_qcheck_nans(q, a, b, ctx, status)) {
3663 mpd_qcopy(r, q, status);
3664 return;
3665 }
3666 if (mpd_isinfinite(a)) {
3667 if (mpd_isinfinite(b)) {
3668 mpd_setspecial(q, MPD_POS, MPD_NAN);
3669 }
3670 else {
3671 mpd_setspecial(q, sign, MPD_INF);
3672 }
3673 mpd_setspecial(r, MPD_POS, MPD_NAN);
3674 *status |= MPD_Invalid_operation;
3675 return;
3676 }
3677 if (mpd_isinfinite(b)) {
3678 if (!mpd_qcopy(r, a, status)) {
3679 mpd_seterror(q, MPD_Malloc_error, status);
3680 return;
3681 }
3682 mpd_qfinalize(r, ctx, status);
3683 _settriple(q, sign, 0, 0);
3684 return;
3685 }
3686 /* debug */
3687 abort(); /* GCOV_NOT_REACHED */
3688 }
3689 if (mpd_iszerocoeff(b)) {
3690 if (mpd_iszerocoeff(a)) {
3691 mpd_setspecial(q, MPD_POS, MPD_NAN);
3692 mpd_setspecial(r, MPD_POS, MPD_NAN);
3693 *status |= MPD_Division_undefined;
3694 }
3695 else {
3696 mpd_setspecial(q, sign, MPD_INF);
3697 mpd_setspecial(r, MPD_POS, MPD_NAN);
3698 *status |= (MPD_Division_by_zero|MPD_Invalid_operation);
3699 }
3700 return;
3701 }
3702
3703 _mpd_qdivmod(q, r, a, b, ctx, status);
3704 mpd_qfinalize(q, ctx, status);
3705 mpd_qfinalize(r, ctx, status);
3706 }
3707
3708 void
3709 mpd_qdivint(mpd_t *q, const mpd_t *a, const mpd_t *b,
3710 const mpd_context_t *ctx, uint32_t *status)
3711 {
3712 MPD_NEW_STATIC(r,0,0,0,0);
3713 uint8_t sign = mpd_sign(a)^mpd_sign(b);
3714
3715 if (mpd_isspecial(a) || mpd_isspecial(b)) {
3716 if (mpd_qcheck_nans(q, a, b, ctx, status)) {
3717 return;
3718 }
3719 if (mpd_isinfinite(a) && mpd_isinfinite(b)) {
3720 mpd_seterror(q, MPD_Invalid_operation, status);
3721 return;
3722 }
3723 if (mpd_isinfinite(a)) {
3724 mpd_setspecial(q, sign, MPD_INF);
3725 return;
3726 }
3727 if (mpd_isinfinite(b)) {
3728 _settriple(q, sign, 0, 0);
3729 return;
3730 }
3731 /* debug */
3732 abort(); /* GCOV_NOT_REACHED */
3733 }
3734 if (mpd_iszerocoeff(b)) {
3735 if (mpd_iszerocoeff(a)) {
3736 mpd_seterror(q, MPD_Division_undefined, status);
3737 }
3738 else {
3739 mpd_setspecial(q, sign, MPD_INF);
3740 *status |= MPD_Division_by_zero;
3741 }
3742 return;
3743 }
3744
3745
3746 _mpd_qdivmod(q, &r, a, b, ctx, status);
3747 mpd_del(&r);
3748 mpd_qfinalize(q, ctx, status);
3749 }
3750
3751 /* Divide decimal by mpd_ssize_t. */
3752 void
3753 mpd_qdiv_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
3754 const mpd_context_t *ctx, uint32_t *status)
3755 {
3756 mpd_context_t maxcontext;
3757 MPD_NEW_STATIC(bb,0,0,0,0);
3758
3759 mpd_maxcontext(&maxcontext);
3760 mpd_qsset_ssize(&bb, b, &maxcontext, status);
3761 mpd_qdiv(result, a, &bb, ctx, status);
3762 mpd_del(&bb);
3763 }
3764
3765 /* Divide decimal by mpd_uint_t. */
3766 void
3767 mpd_qdiv_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
3768 const mpd_context_t *ctx, uint32_t *status)
3769 {
3770 mpd_context_t maxcontext;
3771 MPD_NEW_STATIC(bb,0,0,0,0);
3772
3773 mpd_maxcontext(&maxcontext);
3774 mpd_qsset_uint(&bb, b, &maxcontext, status);
3775 mpd_qdiv(result, a, &bb, ctx, status);
3776 mpd_del(&bb);
3777 }
3778
3779 /* Divide decimal by int32_t. */
3780 void
3781 mpd_qdiv_i32(mpd_t *result, const mpd_t *a, int32_t b,
3782 const mpd_context_t *ctx, uint32_t *status)
3783 {
3784 mpd_qdiv_ssize(result, a, b, ctx, status);
3785 }
3786
3787 /* Divide decimal by uint32_t. */
3788 void
3789 mpd_qdiv_u32(mpd_t *result, const mpd_t *a, uint32_t b,
3790 const mpd_context_t *ctx, uint32_t *status)
3791 {
3792 mpd_qdiv_uint(result, a, b, ctx, status);
3793 }
3794
3795 #ifdef CONFIG_64
3796 /* Divide decimal by int64_t. */
3797 void
3798 mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b,
3799 const mpd_context_t *ctx, uint32_t *status)
3800 {
3801 mpd_qdiv_ssize(result, a, b, ctx, status);
3802 }
3803
3804 /* Divide decimal by uint64_t. */
3805 void
3806 mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3807 const mpd_context_t *ctx, uint32_t *status)
3808 {
3809 mpd_qdiv_uint(result, a, b, ctx, status);
3810 }
3811 #endif
3812
3813 #if defined(_MSC_VER)
3814 /* conversion from 'double' to 'mpd_ssize_t', possible loss of data */
3815 #pragma warning(disable:4244)
3816 #endif
3817 /*
3818 * Get the number of iterations for the Horner scheme in _mpd_qexp().
3819 */
3820 static inline mpd_ssize_t
3821 _mpd_get_exp_iterations(const mpd_t *a, mpd_ssize_t prec)
3822 {
3823 mpd_uint_t dummy;
3824 mpd_uint_t msdigits;
3825 double f;
3826
3827 /* 9 is MPD_RDIGITS for 32 bit platforms */
3828 _mpd_get_msdigits(&dummy, &msdigits, a, 9);
3829 f = ((double)msdigits + 1) / mpd_pow10[mpd_word_digits(msdigits)];
3830
3831 #ifdef CONFIG_64
3832 #ifdef USE_80BIT_LONG_DOUBLE
3833 return ceill((1.435*(long double)prec - 1.182)
3834 / log10l((long double)prec/f));
3835 #else
3836 /* prec > floor((1ULL<<53) / 1.435) */
3837 if (prec > 6276793905742851LL) {
3838 return MPD_SSIZE_MAX;
3839 }
3840 return ceil((1.435*(double)prec - 1.182) / log10((double)prec/f));
3841 #endif
3842 #else /* CONFIG_32 */
3843 return ceil((1.435*(double)prec - 1.182) / log10((double)prec/f));
3844 #if defined(_MSC_VER)
3845 #pragma warning(default:4244)
3846 #endif
3847 #endif
3848 }
3849
3850 /*
3851 * Internal function, specials have been dealt with.
3852 *
3853 * The algorithm is from Hull&Abrham, Variable Precision Exponential Function,
3854 * ACM Transactions on Mathematical Software, Vol. 12, No. 2, June 1986.
3855 *
3856 * Main differences:
3857 *
3858 * - The number of iterations for the Horner scheme is calculated using the
3859 * C log10() function.
3860 *
3861 * - The analysis for early abortion has been adapted for the mpd_t
3862 * ranges.
3863 */
3864 static void
3865 _mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
3866 uint32_t *status)
3867 {
3868 mpd_context_t workctx;
3869 MPD_NEW_STATIC(tmp,0,0,0,0);
3870 MPD_NEW_STATIC(sum,0,0,0,0);
3871 MPD_NEW_CONST(word,0,0,0,1,1,1);
3872 mpd_ssize_t j, n, t;
3873
3874 assert(!mpd_isspecial(a));
3875
3876 /*
3877 * We are calculating e^x = e^(r*10^t) = (e^r)^(10^t), where r < 1 and t >= 0.
3878 *
3879 * If t > 0, we have:
3880 *
3881 * (1) 0.1 <= r < 1, so e^r >= e^0.1. Overflow in the final power oper ation
3882 * will occur when (e^0.1)^(10^t) > 10^(emax+1). If we consider MA X_EMAX,
3883 * this will happen for t > 10 (32 bit) or (t > 19) (64 bit).
3884 *
3885 * (2) -1 < r <= -0.1, so e^r > e^-1. Underflow in the final power ope ration
3886 * will occur when (e^-1)^(10^t) < 10^(etiny-1). If we consider MI N_ETINY,
3887 * this will also happen for t > 10 (32 bit) or (t > 19) (64 bit).
3888 */
3889 #if defined(CONFIG_64)
3890 #define MPD_EXP_MAX_T 19
3891 #elif defined(CONFIG_32)
3892 #define MPD_EXP_MAX_T 10
3893 #endif
3894 t = a->digits + a->exp;
3895 t = (t > 0) ? t : 0;
3896 if (t > MPD_EXP_MAX_T) {
3897 if (mpd_ispositive(a)) {
3898 mpd_setspecial(result, MPD_POS, MPD_INF);
3899 *status |= MPD_Overflow|MPD_Inexact|MPD_Rounded;
3900 }
3901 else {
3902 _settriple(result, MPD_POS, 0, mpd_etiny(ctx));
3903 *status |= (MPD_Inexact|MPD_Rounded|MPD_Subnormal|
3904 MPD_Underflow|MPD_Clamped);
3905 }
3906 return;
3907 }
3908
3909 mpd_maxcontext(&workctx);
3910 workctx.prec = ctx->prec + t + 2;
3911 workctx.prec = (workctx.prec < 9) ? 9 : workctx.prec;
3912 workctx.round = MPD_ROUND_HALF_EVEN;
3913
3914 if ((n = _mpd_get_exp_iterations(a, workctx.prec)) == MPD_SSIZE_MAX) {
3915 mpd_seterror(result, MPD_Invalid_operation, status); /* GCOV_UNL IKELY */
3916 goto finish; /* GCOV_UNLIKELY */
3917 }
3918
3919 if (!mpd_qcopy(result, a, status)) {
3920 goto finish;
3921 }
3922 result->exp -= t;
3923
3924 _settriple(&sum, MPD_POS, 1, 0);
3925
3926 for (j = n-1; j >= 1; j--) {
3927 word.data[0] = j;
3928 mpd_setdigits(&word);
3929 mpd_qdiv(&tmp, result, &word, &workctx, &workctx.status);
3930 mpd_qmul(&sum, &sum, &tmp, &workctx, &workctx.status);
3931 mpd_qadd(&sum, &sum, &one, &workctx, &workctx.status);
3932 }
3933
3934 #ifdef CONFIG_64
3935 _mpd_qpow_uint(result, &sum, mpd_pow10[t], MPD_POS, &workctx, status);
3936 #else
3937 if (t <= MPD_MAX_POW10) {
3938 _mpd_qpow_uint(result, &sum, mpd_pow10[t], MPD_POS, &workctx, st atus);
3939 }
3940 else {
3941 t -= MPD_MAX_POW10;
3942 _mpd_qpow_uint(&tmp, &sum, mpd_pow10[MPD_MAX_POW10], MPD_POS,
3943 &workctx, status);
3944 _mpd_qpow_uint(result, &tmp, mpd_pow10[t], MPD_POS, &workctx, st atus);
3945 }
3946 #endif
3947
3948
3949 finish:
3950 mpd_del(&tmp);
3951 mpd_del(&sum);
3952 *status |= (workctx.status&MPD_Errors);
3953 *status |= (MPD_Inexact|MPD_Rounded);
3954 }
3955
3956 /* exp(a) */
3957 void
3958 mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
3959 uint32_t *status)
3960 {
3961 mpd_context_t workctx;
3962
3963 if (mpd_isspecial(a)) {
3964 if (mpd_qcheck_nan(result, a, ctx, status)) {
3965 return;
3966 }
3967 if (mpd_isnegative(a)) {
3968 _settriple(result, MPD_POS, 0, 0);
3969 }
3970 else {
3971 mpd_setspecial(result, MPD_POS, MPD_INF);
3972 }
3973 return;
3974 }
3975 if (mpd_iszerocoeff(a)) {
3976 _settriple(result, MPD_POS, 1, 0);
3977 return;
3978 }
3979
3980 workctx = *ctx;
3981 workctx.round = MPD_ROUND_HALF_EVEN;
3982
3983 if (ctx->allcr) {
3984 MPD_NEW_STATIC(t1, 0,0,0,0);
3985 MPD_NEW_STATIC(t2, 0,0,0,0);
3986 MPD_NEW_STATIC(ulp, 0,0,0,0);
3987 MPD_NEW_STATIC(aa, 0,0,0,0);
3988 mpd_ssize_t prec;
3989
3990 if (result == a) {
3991 if (!mpd_qcopy(&aa, a, status)) {
3992 mpd_seterror(result, MPD_Malloc_error, status);
3993 return;
3994 }
3995 a = &aa;
3996 }
3997
3998 workctx.clamp = 0;
3999 prec = ctx->prec + 3;
4000 while (1) {
4001 workctx.prec = prec;
4002 _mpd_qexp(result, a, &workctx, status);
4003 _ssettriple(&ulp, MPD_POS, 1,
4004 result->exp + result->digits-workctx.prec-1) ;
4005
4006 workctx.prec = ctx->prec;
4007 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
4008 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
4009 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
4010 mpd_qcmp(&t1, &t2, status) == 0) {
4011 workctx.clamp = ctx->clamp;
4012 mpd_check_underflow(result, &workctx, status);
4013 mpd_qfinalize(result, &workctx, status);
4014 break;
4015 }
4016 prec += MPD_RDIGITS;
4017 }
4018 mpd_del(&t1);
4019 mpd_del(&t2);
4020 mpd_del(&ulp);
4021 mpd_del(&aa);
4022 }
4023 else {
4024 _mpd_qexp(result, a, &workctx, status);
4025 mpd_check_underflow(result, &workctx, status);
4026 mpd_qfinalize(result, &workctx, status);
4027 }
4028 }
4029
4030 /* Fused multiply-add: (a * b) + c, with a single final rounding. */
4031 void
4032 mpd_qfma(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c,
4033 const mpd_context_t *ctx, uint32_t *status)
4034 {
4035 uint32_t workstatus = 0;
4036 mpd_t *cc = (mpd_t *)c;
4037
4038 if (result == c) {
4039 if ((cc = mpd_qncopy(c)) == NULL) {
4040 mpd_seterror(result, MPD_Malloc_error, status);
4041 return;
4042 }
4043 }
4044
4045 _mpd_qmul(result, a, b, ctx, &workstatus);
4046 if (!(workstatus&MPD_Invalid_operation)) {
4047 mpd_qadd(result, result, cc, ctx, &workstatus);
4048 }
4049
4050 if (cc != c) mpd_del(cc);
4051 *status |= workstatus;
4052 }
4053
4054 static inline int
4055 ln_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2], mpd_ssize_t maxprec,
4056 mpd_ssize_t initprec)
4057 {
4058 mpd_ssize_t k;
4059 int i;
4060
4061 assert(maxprec >= 2 && initprec >= 2);
4062 if (maxprec <= initprec) return -1;
4063
4064 i = 0; k = maxprec;
4065 do {
4066 k = (k+2) / 2;
4067 klist[i++] = k;
4068 } while (k > initprec);
4069
4070 return i-1;
4071 }
4072
4073 /* Two word initial approximations for ln(10) */
4074 #ifdef CONFIG_64
4075 #if MPD_RDIGITS != 19
4076 #error "mpdecimal.c: MPD_RDIGITS must be 19."
4077 #endif
4078 static mpd_uint_t mpd_ln10_data[MPD_MINALLOC_MAX] = {
4079 179914546843642076, 2302585092994045684
4080 };
4081 static mpd_uint_t mpd_ln10_init[2] = {
4082 179914546843642076, 2302585092994045684
4083 };
4084 #else
4085 #if MPD_RDIGITS != 9
4086 #error "mpdecimal.c: MPD_RDIGITS must be 9."
4087 #endif
4088 static mpd_uint_t mpd_ln10_data[MPD_MINALLOC_MAX] = {299404568, 230258509};
4089 static mpd_uint_t mpd_ln10_init[2] = {299404568, 230258509};
4090 #endif
4091 /* mpd_ln10 is cached in order to speed up computations */
4092 mpd_t mpd_ln10 = {MPD_STATIC|MPD_STATIC_DATA, -(2*MPD_RDIGITS-1),
4093 2*MPD_RDIGITS, 2, MPD_MINALLOC_MAX, mpd_ln10_data};
4094
4095 static void
4096 mpd_reset_ln10(void)
4097 {
4098 if (mpd_isdynamic_data(&mpd_ln10)) {
4099 mpd_free(mpd_ln10.data);
4100 }
4101 mpd_ln10.data = mpd_ln10_data;
4102 mpd_ln10_data[0] = mpd_ln10_init[0];
4103 mpd_ln10_data[1] = mpd_ln10_init[1];
4104 mpd_ln10.flags = MPD_STATIC|MPD_STATIC_DATA;
4105 mpd_ln10.exp = -(2*MPD_RDIGITS-1);
4106 mpd_ln10.digits = 2*MPD_RDIGITS;
4107 mpd_ln10.len = 2;
4108 mpd_ln10.alloc = MPD_MINALLOC_MAX;
4109 }
4110
4111 /*
4112 * Initializes or updates mpd_ln10. If mpd_ln10 is cached and has exactly the
4113 * requested precision, the function returns. If the cached precision is greater
4114 * than the requested precision, mpd_ln10 is shifted to the requested precision.
4115 *
4116 * The function can fail with MPD_Malloc_error.
4117 */
4118 void
4119 mpd_update_ln10(mpd_ssize_t maxprec, uint32_t *status)
4120 {
4121 mpd_context_t varcontext, maxcontext;
4122 MPD_NEW_STATIC(tmp, 0,0,0,0);
4123 MPD_NEW_CONST(static10, 0,0,2,1,1,10);
4124 mpd_ssize_t klist[MPD_MAX_PREC_LOG2];
4125 int i;
4126
4127 if (mpd_isspecial(&mpd_ln10)) {
4128 mpd_reset_ln10();
4129 }
4130
4131 if (mpd_ln10.digits > maxprec) {
4132 /* shift to smaller cannot fail */
4133 mpd_qshiftr_inplace(&mpd_ln10, mpd_ln10.digits-maxprec);
4134 mpd_ln10.exp = -(mpd_ln10.digits-1);
4135 return;
4136 }
4137 else if (mpd_ln10.digits == maxprec) {
4138 return;
4139 }
4140
4141 mpd_maxcontext(&maxcontext);
4142 mpd_maxcontext(&varcontext);
4143 varcontext.round = MPD_ROUND_TRUNC;
4144
4145 i = ln_schedule_prec(klist, maxprec+2, mpd_ln10.digits);
4146 for (; i >= 0; i--) {
4147 varcontext.prec = 2*klist[i]+3;
4148 mpd_ln10.flags ^= MPD_NEG;
4149 _mpd_qexp(&tmp, &mpd_ln10, &varcontext, status);
4150 mpd_ln10.flags ^= MPD_NEG;
4151 mpd_qmul(&tmp, &static10, &tmp, &varcontext, status);
4152 mpd_qsub(&tmp, &tmp, &one, &maxcontext, status);
4153 mpd_qadd(&mpd_ln10, &mpd_ln10, &tmp, &maxcontext, status);
4154 if (mpd_isspecial(&mpd_ln10)) {
4155 break;
4156 }
4157 }
4158
4159 mpd_del(&tmp);
4160 varcontext.prec = maxprec;
4161 varcontext.round = MPD_ROUND_HALF_EVEN;
4162 mpd_qfinalize(&mpd_ln10, &varcontext, status);
4163 }
4164
4165 /* Initial approximations for the ln() iteration */
4166 static const uint16_t lnapprox[900] = {
4167 /* index 0 - 400: log((i+100)/100) * 1000 */
4168 0, 10, 20, 30, 39, 49, 58, 68, 77, 86, 95, 104, 113, 122, 131, 140, 148, 157,
4169 166, 174, 182, 191, 199, 207, 215, 223, 231, 239, 247, 255, 262, 270, 278,
4170 285, 293, 300, 308, 315, 322, 329, 336, 344, 351, 358, 365, 372, 378, 385,
4171 392, 399, 406, 412, 419, 425, 432, 438, 445, 451, 457, 464, 470, 476, 482,
4172 489, 495, 501, 507, 513, 519, 525, 531, 536, 542, 548, 554, 560, 565, 571,
4173 577, 582, 588, 593, 599, 604, 610, 615, 621, 626, 631, 637, 642, 647, 652,
4174 658, 663, 668, 673, 678, 683, 688, 693, 698, 703, 708, 713, 718, 723, 728,
4175 732, 737, 742, 747, 751, 756, 761, 766, 770, 775, 779, 784, 788, 793, 798,
4176 802, 806, 811, 815, 820, 824, 829, 833, 837, 842, 846, 850, 854, 859, 863,
4177 867, 871, 876, 880, 884, 888, 892, 896, 900, 904, 908, 912, 916, 920, 924,
4178 928, 932, 936, 940, 944, 948, 952, 956, 959, 963, 967, 971, 975, 978, 982,
4179 986, 990, 993, 997, 1001, 1004, 1008, 1012, 1015, 1019, 1022, 1026, 1030,
4180 1033, 1037, 1040, 1044, 1047, 1051, 1054, 1058, 1061, 1065, 1068, 1072, 1075,
4181 1078, 1082, 1085, 1089, 1092, 1095, 1099, 1102, 1105, 1109, 1112, 1115, 1118,
4182 1122, 1125, 1128, 1131, 1135, 1138, 1141, 1144, 1147, 1151, 1154, 1157, 1160,
4183 1163, 1166, 1169, 1172, 1176, 1179, 1182, 1185, 1188, 1191, 1194, 1197, 1200,
4184 1203, 1206, 1209, 1212, 1215, 1218, 1221, 1224, 1227, 1230, 1233, 1235, 1238,
4185 1241, 1244, 1247, 1250, 1253, 1256, 1258, 1261, 1264, 1267, 1270, 1273, 1275,
4186 1278, 1281, 1284, 1286, 1289, 1292, 1295, 1297, 1300, 1303, 1306, 1308, 1311,
4187 1314, 1316, 1319, 1322, 1324, 1327, 1330, 1332, 1335, 1338, 1340, 1343, 1345,
4188 1348, 1351, 1353, 1356, 1358, 1361, 1364, 1366, 1369, 1371, 1374, 1376, 1379,
4189 1381, 1384, 1386, 1389, 1391, 1394, 1396, 1399, 1401, 1404, 1406, 1409, 1411,
4190 1413, 1416, 1418, 1421, 1423, 1426, 1428, 1430, 1433, 1435, 1437, 1440, 1442,
4191 1445, 1447, 1449, 1452, 1454, 1456, 1459, 1461, 1463, 1466, 1468, 1470, 1472,
4192 1475, 1477, 1479, 1482, 1484, 1486, 1488, 1491, 1493, 1495, 1497, 1500, 1502,
4193 1504, 1506, 1509, 1511, 1513, 1515, 1517, 1520, 1522, 1524, 1526, 1528, 1530,
4194 1533, 1535, 1537, 1539, 1541, 1543, 1545, 1548, 1550, 1552, 1554, 1556, 1558,
4195 1560, 1562, 1564, 1567, 1569, 1571, 1573, 1575, 1577, 1579, 1581, 1583, 1585,
4196 1587, 1589, 1591, 1593, 1595, 1597, 1599, 1601, 1603, 1605, 1607, 1609,
4197 /* index 401 - 899: -log((i+100)/1000) * 1000 */
4198 691, 689, 687, 685, 683, 681, 679, 677, 675, 673, 671, 669, 668, 666, 664,
4199 662, 660, 658, 656, 654, 652, 650, 648, 646, 644, 642, 641, 639, 637, 635,
4200 633, 631, 629, 627, 626, 624, 622, 620, 618, 616, 614, 612, 611, 609, 607,
4201 605, 603, 602, 600, 598, 596, 594, 592, 591, 589, 587, 585, 583, 582, 580,
4202 578, 576, 574, 573, 571, 569, 567, 566, 564, 562, 560, 559, 557, 555, 553,
4203 552, 550, 548, 546, 545, 543, 541, 540, 538, 536, 534, 533, 531, 529, 528,
4204 526, 524, 523, 521, 519, 518, 516, 514, 512, 511, 509, 508, 506, 504, 502,
4205 501, 499, 498, 496, 494, 493, 491, 489, 488, 486, 484, 483, 481, 480, 478,
4206 476, 475, 473, 472, 470, 468, 467, 465, 464, 462, 460, 459, 457, 456, 454,
4207 453, 451, 449, 448, 446, 445, 443, 442, 440, 438, 437, 435, 434, 432, 431,
4208 429, 428, 426, 425, 423, 422, 420, 419, 417, 416, 414, 412, 411, 410, 408,
4209 406, 405, 404, 402, 400, 399, 398, 396, 394, 393, 392, 390, 389, 387, 386,
4210 384, 383, 381, 380, 378, 377, 375, 374, 372, 371, 370, 368, 367, 365, 364,
4211 362, 361, 360, 358, 357, 355, 354, 352, 351, 350, 348, 347, 345, 344, 342,
4212 341, 340, 338, 337, 336, 334, 333, 331, 330, 328, 327, 326, 324, 323, 322,
4213 320, 319, 318, 316, 315, 313, 312, 311, 309, 308, 306, 305, 304, 302, 301,
4214 300, 298, 297, 296, 294, 293, 292, 290, 289, 288, 286, 285, 284, 282, 281,
4215 280, 278, 277, 276, 274, 273, 272, 270, 269, 268, 267, 265, 264, 263, 261,
4216 260, 259, 258, 256, 255, 254, 252, 251, 250, 248, 247, 246, 245, 243, 242,
4217 241, 240, 238, 237, 236, 234, 233, 232, 231, 229, 228, 227, 226, 224, 223,
4218 222, 221, 219, 218, 217, 216, 214, 213, 212, 211, 210, 208, 207, 206, 205,
4219 203, 202, 201, 200, 198, 197, 196, 195, 194, 192, 191, 190, 189, 188, 186,
4220 185, 184, 183, 182, 180, 179, 178, 177, 176, 174, 173, 172, 171, 170, 168,
4221 167, 166, 165, 164, 162, 161, 160, 159, 158, 157, 156, 154, 153, 152, 151,
4222 150, 148, 147, 146, 145, 144, 143, 142, 140, 139, 138, 137, 136, 135, 134,
4223 132, 131, 130, 129, 128, 127, 126, 124, 123, 122, 121, 120, 119, 118, 116,
4224 115, 114, 113, 112, 111, 110, 109, 108, 106, 105, 104, 103, 102, 101, 100,
4225 99, 98, 97, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, 84, 83, 82, 81, 80, 79,
4226 78, 77, 76, 75, 74, 73, 72, 70, 69, 68, 67, 66, 65, 64, 63, 62, 61, 60, 59,
4227 58, 57, 56, 54, 53, 52, 51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39,
4228 38, 37, 36, 35, 34, 33, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19,
4229 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1
4230 };
4231
4232 /* Internal ln() function that does not check for specials, zero or one. */
4233 static void
4234 _mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4235 uint32_t *status)
4236 {
4237 mpd_context_t varcontext, maxcontext;
4238 mpd_t *z = (mpd_t *) result;
4239 MPD_NEW_STATIC(v,0,0,0,0);
4240 MPD_NEW_STATIC(vtmp,0,0,0,0);
4241 MPD_NEW_STATIC(tmp,0,0,0,0);
4242 mpd_ssize_t klist[MPD_MAX_PREC_LOG2];
4243 mpd_ssize_t maxprec, shift, t;
4244 mpd_ssize_t a_digits, a_exp;
4245 mpd_uint_t dummy, x;
4246 int i;
4247
4248 assert(!mpd_isspecial(a) && !mpd_iszerocoeff(a));
4249
4250 /*
4251 * We are calculating ln(a) = ln(v * 10^t) = ln(v) + t*ln(10),
4252 * where 0.5 < v <= 5.
4253 */
4254 if (!mpd_qcopy(&v, a, status)) {
4255 mpd_seterror(result, MPD_Malloc_error, status);
4256 goto finish;
4257 }
4258
4259 /* Initial approximation: we have at least one non-zero digit */
4260 _mpd_get_msdigits(&dummy, &x, &v, 3);
4261 if (x < 10) x *= 10;
4262 if (x < 100) x *= 10;
4263 x -= 100;
4264
4265 /* a may equal z */
4266 a_digits = a->digits;
4267 a_exp = a->exp;
4268
4269 mpd_minalloc(z);
4270 mpd_clear_flags(z);
4271 z->data[0] = lnapprox[x];
4272 z->len = 1;
4273 z->exp = -3;
4274 mpd_setdigits(z);
4275
4276 if (x <= 400) {
4277 v.exp = -(a_digits - 1);
4278 t = a_exp + a_digits - 1;
4279 }
4280 else {
4281 v.exp = -a_digits;
4282 t = a_exp + a_digits;
4283 mpd_set_negative(z);
4284 }
4285
4286 mpd_maxcontext(&maxcontext);
4287 mpd_maxcontext(&varcontext);
4288 varcontext.round = MPD_ROUND_TRUNC;
4289
4290 maxprec = ctx->prec + 2;
4291 if (x <= 10 || x >= 805) {
4292 /* v is close to 1: Estimate the magnitude of the logarithm.
4293 * If v = 1 or ln(v) will underflow, skip the loop. Otherwise,
4294 * adjust the precision upwards in order to obtain a sufficient
4295 * number of significant digits.
4296 *
4297 * 1) x/(1+x) < ln(1+x) < x, for x > -1, x != 0
4298 *
4299 * 2) (v-1)/v < ln(v) < v-1
4300 */
4301 mpd_t *lower = &tmp;
4302 mpd_t *upper = &vtmp;
4303 int cmp = _mpd_cmp(&v, &one);
4304
4305 varcontext.round = MPD_ROUND_CEILING;
4306 varcontext.prec = maxprec;
4307 mpd_qsub(upper, &v, &one, &varcontext, &varcontext.status);
4308 varcontext.round = MPD_ROUND_FLOOR;
4309 mpd_qdiv(lower, upper, &v, &varcontext, &varcontext.status);
4310 varcontext.round = MPD_ROUND_TRUNC;
4311
4312 if (cmp < 0) {
4313 _mpd_ptrswap(&upper, &lower);
4314 }
4315 if (mpd_adjexp(upper) < mpd_etiny(ctx)) {
4316 _settriple(z, (cmp<0), 1, mpd_etiny(ctx)-1);
4317 goto postloop;
4318 }
4319 if (mpd_adjexp(lower) < 0) {
4320 maxprec = maxprec - mpd_adjexp(lower);
4321 }
4322 }
4323
4324 i = ln_schedule_prec(klist, maxprec, 2);
4325 for (; i >= 0; i--) {
4326 varcontext.prec = 2*klist[i]+3;
4327 z->flags ^= MPD_NEG;
4328 _mpd_qexp(&tmp, z, &varcontext, status);
4329 z->flags ^= MPD_NEG;
4330
4331 if (v.digits > varcontext.prec) {
4332 shift = v.digits - varcontext.prec;
4333 mpd_qshiftr(&vtmp, &v, shift, status);
4334 vtmp.exp += shift;
4335 mpd_qmul(&tmp, &vtmp, &tmp, &varcontext, status);
4336 }
4337 else {
4338 mpd_qmul(&tmp, &v, &tmp, &varcontext, status);
4339 }
4340
4341 mpd_qsub(&tmp, &tmp, &one, &maxcontext, status);
4342 mpd_qadd(z, z, &tmp, &maxcontext, status);
4343 if (mpd_isspecial(z)) {
4344 break;
4345 }
4346 }
4347
4348 postloop:
4349 mpd_update_ln10(maxprec+2, status);
4350 mpd_qmul_ssize(&tmp, &mpd_ln10, t, &maxcontext, status);
4351 varcontext.prec = maxprec+2;
4352 mpd_qadd(result, &tmp, z, &varcontext, status);
4353
4354
4355 finish:
4356 mpd_del(&v);
4357 mpd_del(&vtmp);
4358 mpd_del(&tmp);
4359 }
4360
4361 /* ln(a) */
4362 void
4363 mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4364 uint32_t *status)
4365 {
4366 mpd_context_t workctx;
4367 mpd_ssize_t adjexp, t;
4368
4369 if (mpd_isspecial(a)) {
4370 if (mpd_qcheck_nan(result, a, ctx, status)) {
4371 return;
4372 }
4373 if (mpd_isnegative(a)) {
4374 mpd_seterror(result, MPD_Invalid_operation, status);
4375 return;
4376 }
4377 mpd_setspecial(result, MPD_POS, MPD_INF);
4378 return;
4379 }
4380 if (mpd_iszerocoeff(a)) {
4381 mpd_setspecial(result, MPD_NEG, MPD_INF);
4382 return;
4383 }
4384 if (mpd_isnegative(a)) {
4385 mpd_seterror(result, MPD_Invalid_operation, status);
4386 return;
4387 }
4388 if (_mpd_cmp(a, &one) == 0) {
4389 _settriple(result, MPD_POS, 0, 0);
4390 return;
4391 }
4392 /* Check if the result will overflow.
4393 *
4394 * 1) adjexp(a) + 1 > log10(a) >= adjexp(a)
4395 *
4396 * 2) |log10(a)| >= adjexp(a), if adjexp(a) >= 0
4397 * |log10(a)| > -adjexp(a)-1, if adjexp(a) < 0
4398 *
4399 * 3) |log(a)| > 2*|log10(a)|
4400 */
4401 adjexp = mpd_adjexp(a);
4402 t = (adjexp < 0) ? -adjexp-1 : adjexp;
4403 t *= 2;
4404 if (mpd_exp_digits(t)-1 > ctx->emax) {
4405 *status |= MPD_Overflow|MPD_Inexact|MPD_Rounded;
4406 mpd_setspecial(result, (adjexp<0), MPD_INF);
4407 return;
4408 }
4409
4410 workctx = *ctx;
4411 workctx.round = MPD_ROUND_HALF_EVEN;
4412
4413 if (ctx->allcr) {
4414 MPD_NEW_STATIC(t1, 0,0,0,0);
4415 MPD_NEW_STATIC(t2, 0,0,0,0);
4416 MPD_NEW_STATIC(ulp, 0,0,0,0);
4417 MPD_NEW_STATIC(aa, 0,0,0,0);
4418 mpd_ssize_t prec;
4419
4420 if (result == a) {
4421 if (!mpd_qcopy(&aa, a, status)) {
4422 mpd_seterror(result, MPD_Malloc_error, status);
4423 return;
4424 }
4425 a = &aa;
4426 }
4427
4428 workctx.clamp = 0;
4429 prec = ctx->prec + 3;
4430 while (1) {
4431 workctx.prec = prec;
4432 _mpd_qln(result, a, &workctx, status);
4433 _ssettriple(&ulp, MPD_POS, 1,
4434 result->exp + result->digits-workctx.prec-1) ;
4435
4436 workctx.prec = ctx->prec;
4437 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
4438 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
4439 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
4440 mpd_qcmp(&t1, &t2, status) == 0) {
4441 workctx.clamp = ctx->clamp;
4442 mpd_check_underflow(result, &workctx, status);
4443 mpd_qfinalize(result, &workctx, status);
4444 break;
4445 }
4446 prec += MPD_RDIGITS;
4447 }
4448 mpd_del(&t1);
4449 mpd_del(&t2);
4450 mpd_del(&ulp);
4451 mpd_del(&aa);
4452 }
4453 else {
4454 _mpd_qln(result, a, &workctx, status);
4455 mpd_check_underflow(result, &workctx, status);
4456 mpd_qfinalize(result, &workctx, status);
4457 }
4458 }
4459
4460 /* Internal log10() function that does not check for specials, zero, ... */
4461 static void
4462 _mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4463 uint32_t *status)
4464 {
4465 mpd_context_t workctx;
4466
4467 mpd_maxcontext(&workctx);
4468 workctx.prec = ctx->prec + 3;
4469 _mpd_qln(result, a, &workctx, status);
4470 mpd_update_ln10(workctx.prec, status);
4471
4472 workctx = *ctx;
4473 workctx.round = MPD_ROUND_HALF_EVEN;
4474 _mpd_qdiv(NO_IDEAL_EXP, result, result, &mpd_ln10, &workctx, status);
4475 }
4476
4477 /* log10(a) */
4478 void
4479 mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4480 uint32_t *status)
4481 {
4482 mpd_context_t workctx;
4483 mpd_ssize_t adjexp, t;
4484
4485 workctx = *ctx;
4486 workctx.round = MPD_ROUND_HALF_EVEN;
4487
4488 if (mpd_isspecial(a)) {
4489 if (mpd_qcheck_nan(result, a, ctx, status)) {
4490 return;
4491 }
4492 if (mpd_isnegative(a)) {
4493 mpd_seterror(result, MPD_Invalid_operation, status);
4494 return;
4495 }
4496 mpd_setspecial(result, MPD_POS, MPD_INF);
4497 return;
4498 }
4499 if (mpd_iszerocoeff(a)) {
4500 mpd_setspecial(result, MPD_NEG, MPD_INF);
4501 return;
4502 }
4503 if (mpd_isnegative(a)) {
4504 mpd_seterror(result, MPD_Invalid_operation, status);
4505 return;
4506 }
4507 if (mpd_coeff_ispow10(a)) {
4508 uint8_t sign = 0;
4509 adjexp = mpd_adjexp(a);
4510 if (adjexp < 0) {
4511 sign = 1;
4512 adjexp = -adjexp;
4513 }
4514 _settriple(result, sign, adjexp, 0);
4515 mpd_qfinalize(result, &workctx, status);
4516 return;
4517 }
4518 /* Check if the result will overflow.
4519 *
4520 * 1) adjexp(a) + 1 > log10(a) >= adjexp(a)
4521 *
4522 * 2) |log10(a)| >= adjexp(a), if adjexp(a) >= 0
4523 * |log10(a)| > -adjexp(a)-1, if adjexp(a) < 0
4524 */
4525 adjexp = mpd_adjexp(a);
4526 t = (adjexp < 0) ? -adjexp-1 : adjexp;
4527 if (mpd_exp_digits(t)-1 > ctx->emax) {
4528 *status |= MPD_Overflow|MPD_Inexact|MPD_Rounded;
4529 mpd_setspecial(result, (adjexp<0), MPD_INF);
4530 return;
4531 }
4532
4533 if (ctx->allcr) {
4534 MPD_NEW_STATIC(t1, 0,0,0,0);
4535 MPD_NEW_STATIC(t2, 0,0,0,0);
4536 MPD_NEW_STATIC(ulp, 0,0,0,0);
4537 MPD_NEW_STATIC(aa, 0,0,0,0);
4538 mpd_ssize_t prec;
4539
4540 if (result == a) {
4541 if (!mpd_qcopy(&aa, a, status)) {
4542 mpd_seterror(result, MPD_Malloc_error, status);
4543 return;
4544 }
4545 a = &aa;
4546 }
4547
4548 workctx.clamp = 0;
4549 prec = ctx->prec + 3;
4550 while (1) {
4551 workctx.prec = prec;
4552 _mpd_qlog10(result, a, &workctx, status);
4553 _ssettriple(&ulp, MPD_POS, 1,
4554 result->exp + result->digits-workctx.prec-1) ;
4555
4556 workctx.prec = ctx->prec;
4557 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
4558 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
4559 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
4560 mpd_qcmp(&t1, &t2, status) == 0) {
4561 workctx.clamp = ctx->clamp;
4562 mpd_check_underflow(result, &workctx, status);
4563 mpd_qfinalize(result, &workctx, status);
4564 break;
4565 }
4566 prec += MPD_RDIGITS;
4567 }
4568 mpd_del(&t1);
4569 mpd_del(&t2);
4570 mpd_del(&ulp);
4571 mpd_del(&aa);
4572 }
4573 else {
4574 _mpd_qlog10(result, a, &workctx, status);
4575 mpd_check_underflow(result, &workctx, status);
4576 }
4577 }
4578
4579 /*
4580 * Maximum of the two operands. Attention: If one operand is a quiet NaN and the
4581 * other is numeric, the numeric operand is returned. This may not be what one
4582 * expects.
4583 */
4584 void
4585 mpd_qmax(mpd_t *result, const mpd_t *a, const mpd_t *b,
4586 const mpd_context_t *ctx, uint32_t *status)
4587 {
4588 int c;
4589
4590 if (mpd_isqnan(a) && !mpd_isnan(b)) {
4591 mpd_qcopy(result, b, status);
4592 }
4593 else if (mpd_isqnan(b) && !mpd_isnan(a)) {
4594 mpd_qcopy(result, a, status);
4595 }
4596 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
4597 return;
4598 }
4599 else {
4600 c = _mpd_cmp(a, b);
4601 if (c == 0) {
4602 c = _mpd_cmp_numequal(a, b);
4603 }
4604
4605 if (c < 0) {
4606 mpd_qcopy(result, b, status);
4607 }
4608 else {
4609 mpd_qcopy(result, a, status);
4610 }
4611 }
4612
4613 mpd_qfinalize(result, ctx, status);
4614 }
4615
4616 /*
4617 * Maximum magnitude: Same as mpd_max(), but compares the operands with their
4618 * sign ignored.
4619 */
4620 void
4621 mpd_qmax_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
4622 const mpd_context_t *ctx, uint32_t *status)
4623 {
4624 int c;
4625
4626 if (mpd_isqnan(a) && !mpd_isnan(b)) {
4627 mpd_qcopy(result, b, status);
4628 }
4629 else if (mpd_isqnan(b) && !mpd_isnan(a)) {
4630 mpd_qcopy(result, a, status);
4631 }
4632 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
4633 return;
4634 }
4635 else {
4636 c = _mpd_cmp_abs(a, b);
4637 if (c == 0) {
4638 c = _mpd_cmp_numequal(a, b);
4639 }
4640
4641 if (c < 0) {
4642 mpd_qcopy(result, b, status);
4643 }
4644 else {
4645 mpd_qcopy(result, a, status);
4646 }
4647 }
4648
4649 mpd_qfinalize(result, ctx, status);
4650 }
4651
4652 /*
4653 * Minimum of the two operands. Attention: If one operand is a quiet NaN and the
4654 * other is numeric, the numeric operand is returned. This may not be what one
4655 * expects.
4656 */
4657 void
4658 mpd_qmin(mpd_t *result, const mpd_t *a, const mpd_t *b,
4659 const mpd_context_t *ctx, uint32_t *status)
4660 {
4661 int c;
4662
4663 if (mpd_isqnan(a) && !mpd_isnan(b)) {
4664 mpd_qcopy(result, b, status);
4665 }
4666 else if (mpd_isqnan(b) && !mpd_isnan(a)) {
4667 mpd_qcopy(result, a, status);
4668 }
4669 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
4670 return;
4671 }
4672 else {
4673 c = _mpd_cmp(a, b);
4674 if (c == 0) {
4675 c = _mpd_cmp_numequal(a, b);
4676 }
4677
4678 if (c < 0) {
4679 mpd_qcopy(result, a, status);
4680 }
4681 else {
4682 mpd_qcopy(result, b, status);
4683 }
4684 }
4685
4686 mpd_qfinalize(result, ctx, status);
4687 }
4688
4689 /*
4690 * Minimum magnitude: Same as mpd_min(), but compares the operands with their
4691 * sign ignored.
4692 */
4693 void
4694 mpd_qmin_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
4695 const mpd_context_t *ctx, uint32_t *status)
4696 {
4697 int c;
4698
4699 if (mpd_isqnan(a) && !mpd_isnan(b)) {
4700 mpd_qcopy(result, b, status);
4701 }
4702 else if (mpd_isqnan(b) && !mpd_isnan(a)) {
4703 mpd_qcopy(result, a, status);
4704 }
4705 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
4706 return;
4707 }
4708 else {
4709 c = _mpd_cmp_abs(a, b);
4710 if (c == 0) {
4711 c = _mpd_cmp_numequal(a, b);
4712 }
4713
4714 if (c < 0) {
4715 mpd_qcopy(result, a, status);
4716 }
4717 else {
4718 mpd_qcopy(result, b, status);
4719 }
4720 }
4721
4722 mpd_qfinalize(result, ctx, status);
4723 }
4724
4725 /* Minimum space needed for the result array in _karatsuba_rec(). */
4726 static inline mpd_size_t
4727 _kmul_resultsize(mpd_size_t la, mpd_size_t lb)
4728 {
4729 mpd_size_t n, m;
4730
4731 n = add_size_t(la, lb);
4732 n = add_size_t(n, 1);
4733
4734 m = (la+1)/2 + 1;
4735 m = mul_size_t(m, 3);
4736
4737 return (m > n) ? m : n;
4738 }
4739
4740 /* Work space needed in _karatsuba_rec(). lim >= 4 */
4741 static inline mpd_size_t
4742 _kmul_worksize(mpd_size_t n, mpd_size_t lim)
4743 {
4744 mpd_size_t m;
4745
4746 if (n <= lim) {
4747 return 0;
4748 }
4749
4750 m = (n+1)/2 + 1;
4751
4752 return add_size_t(mul_size_t(m, 2), _kmul_worksize(m, lim));
4753 }
4754
4755
4756 #define MPD_KARATSUBA_BASECASE 16 /* must be >= 4 */
4757
4758 /*
4759 * Add the product of a and b to c.
4760 * c must be _kmul_resultsize(la, lb) in size.
4761 * w is used as a work array and must be _kmul_worksize(a, lim) in size.
4762 * Roman E. Maeder, Storage Allocation for the Karatsuba Integer Multiplication
4763 * Algorithm. In "Design and implementation of symbolic computation systems",
4764 * Springer, 1993, ISBN 354057235X, 9783540572350.
4765 */
4766 static void
4767 _karatsuba_rec(mpd_uint_t *c, const mpd_uint_t *a, const mpd_uint_t *b,
4768 mpd_uint_t *w, mpd_size_t la, mpd_size_t lb)
4769 {
4770 mpd_size_t m, lt;
4771
4772 assert (la >= lb && lb > 0);
4773
4774 if (la <= MPD_KARATSUBA_BASECASE) {
4775 _mpd_basemul(c, a, b, la, lb);
4776 return;
4777 }
4778
4779 m = (la+1)/2; // ceil(la/2)
4780
4781 /* lb <= m < la */
4782 if (lb <= m) {
4783
4784 /* lb can now be larger than la-m */
4785 if (lb > la-m) {
4786 lt = lb + lb + 1; // space needed for result array
4787 mpd_uint_zero(w, lt); // clear result array
4788 _karatsuba_rec(w, b, a+m, w+lt, lb, la-m); // b*ah
4789 }
4790 else {
4791 lt = (la-m) + (la-m) + 1; // space needed for result ar ray
4792 mpd_uint_zero(w, lt); // clear result array
4793 _karatsuba_rec(w, a+m, b, w+lt, la-m, lb); // ah*b
4794 }
4795 _mpd_baseaddto(c+m, w, (la-m)+lb); // add ah*b*B**m
4796
4797 lt = m + m + 1; // space needed for the result array
4798 mpd_uint_zero(w, lt); // clear result array
4799 _karatsuba_rec(w, a, b, w+lt, m, lb); // al*b
4800 _mpd_baseaddto(c, w, m+lb); // add al*b
4801
4802 return;
4803 }
4804
4805 /* la >= lb > m */
4806 memcpy(w, a, m * sizeof *w);
4807 w[m] = 0;
4808 _mpd_baseaddto(w, a+m, la-m);
4809
4810 memcpy(w+(m+1), b, m * sizeof *w);
4811 w[m+1+m] = 0;
4812 _mpd_baseaddto(w+(m+1), b+m, lb-m);
4813
4814 _karatsuba_rec(c+m, w, w+(m+1), w+2*(m+1), m+1, m+1);
4815
4816 lt = (la-m) + (la-m) + 1;
4817 mpd_uint_zero(w, lt);
4818
4819 _karatsuba_rec(w, a+m, b+m, w+lt, la-m, lb-m);
4820
4821 _mpd_baseaddto(c+2*m, w, (la-m) + (lb-m));
4822 _mpd_basesubfrom(c+m, w, (la-m) + (lb-m));
4823
4824 lt = m + m + 1;
4825 mpd_uint_zero(w, lt);
4826
4827 _karatsuba_rec(w, a, b, w+lt, m, m);
4828 _mpd_baseaddto(c, w, m+m);
4829 _mpd_basesubfrom(c+m, w, m+m);
4830
4831 return;
4832 }
4833
4834 /*
4835 * Multiply u and v, using Karatsuba multiplication. Returns a pointer
4836 * to the result or NULL in case of failure (malloc error).
4837 * Conditions: ulen >= vlen, ulen >= 4
4838 */
4839 mpd_uint_t *
4840 _mpd_kmul(const mpd_uint_t *u, const mpd_uint_t *v,
4841 mpd_size_t ulen, mpd_size_t vlen,
4842 mpd_size_t *rsize)
4843 {
4844 mpd_uint_t *result = NULL, *w = NULL;
4845 mpd_size_t m;
4846
4847 assert(ulen >= 4);
4848 assert(ulen >= vlen);
4849
4850 *rsize = _kmul_resultsize(ulen, vlen);
4851 if ((result = mpd_calloc(*rsize, sizeof *result)) == NULL) {
4852 return NULL;
4853 }
4854
4855 m = _kmul_worksize(ulen, MPD_KARATSUBA_BASECASE);
4856 if (m && ((w = mpd_calloc(m, sizeof *w)) == NULL)) {
4857 mpd_free(result);
4858 return NULL;
4859 }
4860
4861 _karatsuba_rec(result, u, v, w, ulen, vlen);
4862
4863
4864 if (w) mpd_free(w);
4865 return result;
4866 }
4867
4868
4869 /* Determine the minimum length for the number theoretic transform. */
4870 static inline mpd_size_t
4871 _mpd_get_transform_len(mpd_size_t rsize)
4872 {
4873 mpd_size_t log2rsize;
4874 mpd_size_t x, step;
4875
4876 assert(rsize >= 4);
4877 log2rsize = mpd_bsr(rsize);
4878
4879 if (rsize <= 1024) {
4880 x = ((mpd_size_t)1)<<log2rsize;
4881 return (rsize == x) ? x : x<<1;
4882 }
4883 else if (rsize <= MPD_MAXTRANSFORM_2N) {
4884 x = ((mpd_size_t)1)<<log2rsize;
4885 if (rsize == x) return x;
4886 step = x>>1;
4887 x += step;
4888 return (rsize <= x) ? x : x + step;
4889 }
4890 else if (rsize <= MPD_MAXTRANSFORM_2N+MPD_MAXTRANSFORM_2N/2) {
4891 return MPD_MAXTRANSFORM_2N+MPD_MAXTRANSFORM_2N/2;
4892 }
4893 else if (rsize <= 3*MPD_MAXTRANSFORM_2N) {
4894 return 3*MPD_MAXTRANSFORM_2N;
4895 }
4896 else {
4897 return MPD_SIZE_MAX;
4898 }
4899 }
4900
4901 #ifdef PPRO
4902 #ifndef _MSC_VER
4903 static inline unsigned short
4904 _mpd_get_control87(void)
4905 {
4906 unsigned short cw;
4907
4908 __asm__ __volatile__ ("fnstcw %0" : "=m" (cw));
4909 return cw;
4910 }
4911
4912 static inline void
4913 _mpd_set_control87(unsigned short cw)