Lines Matching defs:result
106 static void _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a,
112 static void _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
114 static inline void _mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b,
118 static inline void _mpd_qpow_uint(mpd_t *result, const mpd_t *base,
122 static mpd_uint_t mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n);
478 * Input invariant: MPD_MINALLOC <= result->alloc.
480 * Case nwords == result->alloc:
481 * 'result' is unchanged. Return 1.
483 * Case nwords > result->alloc:
485 * The value of 'result' does not change. Return 1.
487 * 'result' is NaN, status is updated with MPD_Malloc_error. Return 0.
489 * Case nwords < result->alloc:
491 * 'result' is unchanged. Return 1.
493 * The value of result is undefined (expected). Return 1.
499 mpd_qresize(mpd_t *result, mpd_ssize_t nwords, uint32_t *status)
501 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
502 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
503 assert(MPD_MINALLOC <= result->alloc);
506 if (nwords == result->alloc) {
509 if (mpd_isstatic_data(result)) {
510 if (nwords > result->alloc) {
511 return mpd_switch_to_dyn(result, nwords, status);
516 return mpd_realloc_dyn(result, nwords, status);
519 /* Same as mpd_qresize, but do not set the result no NaN on failure. */
521 mpd_qresize_cxx(mpd_t *result, mpd_ssize_t nwords)
523 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
524 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
525 assert(MPD_MINALLOC <= result->alloc);
528 if (nwords == result->alloc) {
531 if (mpd_isstatic_data(result)) {
532 if (nwords > result->alloc) {
533 return mpd_switch_to_dyn_cxx(result, nwords);
538 return mpd_realloc_dyn_cxx(result, nwords);
544 mpd_qresize_zero(mpd_t *result, mpd_ssize_t nwords, uint32_t *status)
546 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
547 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
548 assert(MPD_MINALLOC <= result->alloc);
551 if (nwords != result->alloc) {
552 if (mpd_isstatic_data(result)) {
553 if (nwords > result->alloc) {
554 return mpd_switch_to_dyn_zero(result, nwords, status);
557 else if (!mpd_realloc_dyn(result, nwords, status)) {
562 mpd_uint_zero(result->data, nwords);
573 mpd_minalloc(mpd_t *result)
575 assert(!mpd_isconst_data(result)); /* illegal operation for a const */
576 assert(!mpd_isshared_data(result)); /* illegal operation for a shared */
578 if (!mpd_isstatic_data(result) && result->alloc > MPD_MINALLOC) {
580 result->data = mpd_realloc(result->data, MPD_MINALLOC,
581 sizeof *result->data, &err);
583 result->alloc = MPD_MINALLOC;
589 mpd_resize(mpd_t *result, mpd_ssize_t nwords, mpd_context_t *ctx)
592 if (!mpd_qresize(result, nwords, &status)) {
600 mpd_resize_zero(mpd_t *result, mpd_ssize_t nwords, mpd_context_t *ctx)
603 if (!mpd_qresize_zero(result, nwords, &status)) {
615 /* Set digits. Assumption: result->len is initialized and > 0. */
617 mpd_setdigits(mpd_t *result)
619 mpd_ssize_t wdigits = mpd_word_digits(mpd_msword(result));
620 result->digits = wdigits + (result->len-1) * MPD_RDIGITS;
625 mpd_set_sign(mpd_t *result, uint8_t sign)
627 result->flags &= ~MPD_NEG;
628 result->flags |= sign;
633 mpd_signcpy(mpd_t *result, const mpd_t *a)
637 result->flags &= ~MPD_NEG;
638 result->flags |= sign;
643 mpd_set_infinity(mpd_t *result)
645 result->flags &= ~MPD_SPECIAL;
646 result->flags |= MPD_INF;
651 mpd_set_qnan(mpd_t *result)
653 result->flags &= ~MPD_SPECIAL;
654 result->flags |= MPD_NAN;
659 mpd_set_snan(mpd_t *result)
661 result->flags &= ~MPD_SPECIAL;
662 result->flags |= MPD_SNAN;
667 mpd_set_negative(mpd_t *result)
669 result->flags |= MPD_NEG;
674 mpd_set_positive(mpd_t *result)
676 result->flags &= ~MPD_NEG;
681 mpd_set_dynamic(mpd_t *result)
683 result->flags &= ~MPD_STATIC;
688 mpd_set_static(mpd_t *result)
690 result->flags |= MPD_STATIC;
695 mpd_set_dynamic_data(mpd_t *result)
697 result->flags &= ~MPD_DATAFLAGS;
702 mpd_set_static_data(mpd_t *result)
704 result->flags &= ~MPD_DATAFLAGS;
705 result->flags |= MPD_STATIC_DATA;
710 mpd_set_shared_data(mpd_t *result)
712 result->flags &= ~MPD_DATAFLAGS;
713 result->flags |= MPD_SHARED_DATA;
718 mpd_set_const_data(mpd_t *result)
720 result->flags &= ~MPD_DATAFLAGS;
721 result->flags |= MPD_CONST_DATA;
726 mpd_clear_flags(mpd_t *result)
728 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
733 mpd_set_flags(mpd_t *result, uint8_t flags)
735 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
736 result->flags |= flags;
739 /* Copy flags, preserving memory attributes of result. */
741 mpd_copy_flags(mpd_t *result, const mpd_t *a)
744 result->flags &= (MPD_STATIC|MPD_DATAFLAGS);
745 result->flags |= (aflags & ~(MPD_STATIC|MPD_DATAFLAGS));
777 mpd_zerocoeff(mpd_t *result)
779 mpd_minalloc(result);
780 result->digits = 1;
781 result->len = 1;
782 result->data[0] = 0;
787 mpd_qmaxcoeff(mpd_t *result, const mpd_context_t *ctx, uint32_t *status)
794 if (!mpd_qresize(result, len, status)) {
798 result->len = len;
799 result->digits = ctx->prec;
803 result->data[len--] = mpd_pow10[r]-1;
806 result->data[len] = MPD_RADIX-1;
815 _mpd_cap(mpd_t *result, const mpd_context_t *ctx)
820 if (result->len > 0 && result->digits > ctx->prec) {
825 result->data[len-1] %= mpd_pow10[r];
828 len = _mpd_real_size(result->data, len);
830 mpd_qresize(result, len, &dummy);
831 result->len = len;
832 mpd_setdigits(result);
834 if (mpd_iszero(result)) {
835 _settriple(result, mpd_sign(result), 0, result->exp);
844 _mpd_fix_nan(mpd_t *result, const mpd_context_t *ctx)
851 if (result->len > 0 && result->digits > prec) {
853 mpd_minalloc(result);
854 result->len = result->digits = 0;
861 result->data[len-1] %= mpd_pow10[r];
864 len = _mpd_real_size(result->data, len);
866 mpd_qresize(result, len, &dummy);
867 result->len = len;
868 mpd_setdigits(result);
869 if (mpd_iszerocoeff(result)) {
871 result->len = result->digits = 0;
882 * The result of the operation will be in lo. If the operation is impossible,
1054 _ssettriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp)
1056 mpd_set_flags(result, sign);
1057 result->exp = exp;
1058 _mpd_div_word(&result->data[1], &result->data[0], a, MPD_RADIX);
1059 result->len = (result->data[1] == 0) ? 1 : 2;
1060 mpd_setdigits(result);
1065 _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp)
1067 mpd_minalloc(result);
1068 mpd_set_flags(result, sign);
1069 result->exp = exp;
1070 _mpd_div_word(&result->data[1], &result->data[0], a, MPD_RADIX);
1071 result->len = (result->data[1] == 0) ? 1 : 2;
1072 mpd_setdigits(result);
1077 mpd_setspecial(mpd_t *result, uint8_t sign, uint8_t type)
1079 mpd_minalloc(result);
1080 result->flags &= ~(MPD_NEG|MPD_SPECIAL);
1081 result->flags |= (sign|type);
1082 result->exp = result->digits = result->len = 0;
1085 /* Set result of NaN with an error status */
1087 mpd_seterror(mpd_t *result, uint32_t flags, uint32_t *status)
1089 mpd_minalloc(result);
1090 mpd_set_qnan(result);
1091 mpd_set_positive(result);
1092 result->exp = result->digits = result->len = 0;
1098 mpd_qsset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx,
1117 _ssettriple(result, sign, u, 0);
1118 mpd_qfinalize(result, ctx, status);
1123 mpd_qsset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx,
1126 _ssettriple(result, MPD_POS, a, 0);
1127 mpd_qfinalize(result, ctx, status);
1132 mpd_qsset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx,
1135 mpd_qsset_ssize(result, a, ctx, status);
1140 mpd_qsset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx,
1143 mpd_qsset_uint(result, a, ctx, status);
1149 mpd_qsset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1152 mpd_qsset_ssize(result, a, ctx, status);
1157 mpd_qsset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1160 mpd_qsset_uint(result, a, ctx, status);
1166 mpd_qset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx,
1169 mpd_minalloc(result);
1170 mpd_qsset_ssize(result, a, ctx, status);
1175 mpd_qset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx,
1178 _settriple(result, MPD_POS, a, 0);
1179 mpd_qfinalize(result, ctx, status);
1184 mpd_qset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx,
1187 mpd_qset_ssize(result, a, ctx, status);
1192 mpd_qset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx,
1195 mpd_qset_uint(result, a, ctx, status);
1201 _c32setu64(mpd_t *result, uint64_t u, uint8_t sign, uint32_t *status)
1214 if (!mpd_qresize(result, len, status)) {
1218 result->data[i] = w[i];
1221 mpd_set_flags(result, sign);
1222 result->exp = 0;
1223 result->len = len;
1224 mpd_setdigits(result);
1228 _c32_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1231 _c32setu64(result, a, MPD_POS, status);
1232 mpd_qfinalize(result, ctx, status);
1237 _c32_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1255 _c32setu64(result, u, sign, status);
1256 mpd_qfinalize(result, ctx, status);
1263 mpd_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx,
1267 mpd_qset_ssize(result, a, ctx, status);
1269 _c32_qset_i64(result, a, ctx, status);
1275 mpd_qset_i64_exact(mpd_t *result, int64_t a, uint32_t *status)
1281 mpd_qset_ssize(result, a, &maxcontext, status);
1283 _c32_qset_i64(result, a, &maxcontext, status);
1288 mpd_seterror(result, MPD_Invalid_operation, status);
1295 mpd_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx,
1299 mpd_qset_uint(result, a, ctx, status);
1301 _c32_qset_u64(result, a, ctx, status);
1307 mpd_qset_u64_exact(mpd_t *result, uint64_t a, uint32_t *status)
1313 mpd_qset_uint(result, a, &maxcontext, status);
1315 _c32_qset_u64(result, a, &maxcontext, status);
1320 mpd_seterror(result, MPD_Invalid_operation, status);
1616 * Check if the operand is NaN, copy to result and return 1 if this is
1621 mpd_qcheck_nan(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
1626 mpd_qcopy(result, a, status);
1627 mpd_set_qnan(result);
1628 _mpd_fix_nan(result, ctx);
1635 * Check if either operand is NaN, copy to result and return 1 if this
1640 mpd_qcheck_nans(mpd_t *result, const mpd_t *a, const mpd_t *b,
1655 mpd_qcopy(result, choice, status);
1656 mpd_set_qnan(result);
1657 _mpd_fix_nan(result, ctx);
1664 * Check if one of the operands is NaN, copy to result and return 1 if this
1669 mpd_qcheck_3nans(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c,
1691 mpd_qcopy(result, choice, status);
1692 mpd_set_qnan(result);
1693 _mpd_fix_nan(result, ctx);
1742 * that case, the result is a power of ten with prec+1 digits.
1747 * check if the result has one digit too many.
1926 * rounding. If a result like 1.0000000000e-101 is finalized, there
1964 mpd_qfinalize(mpd_t *result, const mpd_context_t *ctx, uint32_t *status)
1966 if (mpd_isspecial(result)) {
1967 if (mpd_isnan(result)) {
1968 _mpd_fix_nan(result, ctx);
1973 _mpd_check_exp(result, ctx, status);
1974 _mpd_check_round(result, ctx, status);
2000 mpd_qcopy(mpd_t *result, const mpd_t *a, uint32_t *status)
2002 if (result == a) return 1;
2004 if (!mpd_qresize(result, a->len, status)) {
2008 mpd_copy_flags(result, a);
2009 result->exp = a->exp;
2010 result->digits = a->digits;
2011 result->len = a->len;
2012 memcpy(result->data, a->data, a->len * (sizeof *result->data));
2017 /* Same as mpd_qcopy, but do not set the result to NaN on failure. */
2019 mpd_qcopy_cxx(mpd_t *result, const mpd_t *a)
2021 if (result == a) return 1;
2023 if (!mpd_qresize_cxx(result, a->len)) {
2027 mpd_copy_flags(result, a);
2028 result->exp = a->exp;
2029 result->digits = a->digits;
2030 result->len = a->len;
2031 memcpy(result->data, a->data, a->len * (sizeof *result->data));
2041 mpd_qcopy_static(mpd_t *result, const mpd_t *a)
2043 if (result == a) return;
2045 memcpy(result->data, a->data, a->len * (sizeof *result->data));
2047 mpd_copy_flags(result, a);
2048 result->exp = a->exp;
2049 result->digits = a->digits;
2050 result->len = a->len;
2060 mpd_t *result;
2062 if ((result = mpd_qnew_size(a->len)) == NULL) {
2065 memcpy(result->data, a->data, a->len * (sizeof *result->data));
2066 mpd_copy_flags(result, a);
2067 result->exp = a->exp;
2068 result->digits = a->digits;
2069 result->len = a->len;
2071 return result;
2079 mpd_qcopy_abs(mpd_t *result, const mpd_t *a, uint32_t *status)
2081 if (!mpd_qcopy(result, a, status)) {
2084 mpd_set_positive(result);
2093 mpd_qcopy_negate(mpd_t *result, const mpd_t *a, uint32_t *status)
2095 if (!mpd_qcopy(result, a, status)) {
2098 _mpd_negate(result);
2107 mpd_qcopy_sign(mpd_t *result, const mpd_t *a, const mpd_t *b, uint32_t *status)
2109 uint8_t sign_b = mpd_sign(b); /* result may equal b! */
2111 if (!mpd_qcopy(result, a, status)) {
2114 mpd_set_sign(result, sign_b);
2323 /* Compare two values and return an integer result. */
2338 * Compare a and b, convert the usual integer result to a decimal and
2339 * store it in 'result'. For convenience, the integer result of the comparison
2343 mpd_qcompare(mpd_t *result, const mpd_t *a, const mpd_t *b,
2349 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2355 _settriple(result, (c < 0), (c != 0), 0);
2361 mpd_qcompare_signal(mpd_t *result, const mpd_t *a, const mpd_t *b,
2367 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2374 _settriple(result, (c < 0), (c != 0), 0);
2427 * Compare a and b according to a total order, convert the usual integer result
2428 * to a decimal and store it in 'result'. For convenience, the integer result
2432 mpd_compare_total(mpd_t *result, const mpd_t *a, const mpd_t *b)
2437 _settriple(result, (c < 0), (c != 0), 0);
2458 * the usual integer result to a decimal and store it in 'result'.
2459 * For convenience, the integer result of the comparison is returned.
2462 mpd_compare_total_mag(mpd_t *result, const mpd_t *a, const mpd_t *b)
2467 _settriple(result, (c < 0), (c != 0), 0);
2498 * Both operands may be the same pointer. If the result length has to be
2502 mpd_qshiftl(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status)
2510 return mpd_qcopy(result, a, status);
2514 if (!mpd_qresize(result, size, status)) {
2515 return 0; /* result is NaN */
2518 _mpd_baseshiftl(result->data, a->data, size, a->len, n);
2520 mpd_copy_flags(result, a);
2521 result->exp = a->exp;
2522 result->digits = a->digits+n;
2523 result->len = size;
2551 * Same as mpd_qshiftr(), but 'result' is an mpd_t with a static coefficient.
2556 mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n)
2565 mpd_qcopy_static(result, a);
2571 mpd_zerocoeff(result);
2574 result->digits = a->digits-n;
2575 size = mpd_digits_to_size(result->digits);
2576 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2577 result->len = size;
2580 mpd_copy_flags(result, a);
2581 result->exp = a->exp;
2592 mpd_qshiftr_inplace(mpd_t *result, mpd_ssize_t n)
2598 assert(!mpd_isspecial(result));
2601 if (mpd_iszerocoeff(result) || n == 0) {
2605 if (n >= result->digits) {
2606 rnd = _mpd_get_rnd(result->data, result->len, (n==result->digits));
2607 mpd_zerocoeff(result);
2610 rnd = _mpd_baseshiftr(result->data, result->data, result->len, n);
2611 result->digits -= n;
2612 size = mpd_digits_to_size(result->digits);
2614 mpd_qresize(result, size, &dummy);
2615 result->len = size;
2624 * be used by mpd_rnd_incr(). If the result length has to be increased,
2629 mpd_qshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status)
2638 if (!mpd_qcopy(result, a, status)) {
2646 mpd_zerocoeff(result);
2649 result->digits = a->digits-n;
2650 size = mpd_digits_to_size(result->digits);
2651 if (result == a) {
2652 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2654 mpd_qresize(result, size, status);
2657 if (!mpd_qresize(result, size, status)) {
2660 rnd = _mpd_baseshiftr(result->data, a->data, a->len, n);
2662 result->len = size;
2665 mpd_copy_flags(result, a);
2666 result->exp = a->exp;
2678 mpd_qand(mpd_t *result, const mpd_t *a, const mpd_t *b,
2689 mpd_seterror(result, MPD_Invalid_operation, status);
2696 if (!mpd_qresize(result, big->len, status)) {
2716 result->data[i] = z;
2733 result->data[i++] = z;
2755 mpd_clear_flags(result);
2756 result->exp = 0;
2757 result->len = _mpd_real_size(result->data, small->len);
2758 mpd_qresize(result, result->len, status);
2759 mpd_setdigits(result);
2760 _mpd_cap(result, ctx);
2764 mpd_seterror(result, MPD_Invalid_operation, status);
2801 mpd_qinvert(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
2810 mpd_seterror(result, MPD_Invalid_operation, status);
2817 if (!mpd_qresize(result, len, status)) {
2832 result->data[i] = z;
2835 mpd_clear_flags(result);
2836 result->exp = 0;
2837 result->len = _mpd_real_size(result->data, len);
2838 mpd_qresize(result, result->len, status);
2839 mpd_setdigits(result);
2840 _mpd_cap(result, ctx);
2844 mpd_seterror(result, MPD_Invalid_operation, status);
2849 mpd_qlogb(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
2853 if (mpd_qcheck_nan(result, a, ctx, status)) {
2856 mpd_setspecial(result, MPD_POS, MPD_INF);
2859 mpd_setspecial(result, MPD_NEG, MPD_INF);
2863 mpd_qset_ssize(result, mpd_adjexp(a), ctx, status);
2869 mpd_qor(mpd_t *result, const mpd_t *a, const mpd_t *b,
2880 mpd_seterror(result, MPD_Invalid_operation, status);
2887 if (!mpd_qresize(result, big->len, status)) {
2907 result->data[i] = z;
2934 result->data[i++] = z;
2945 result->data[i] = big->data[i];
2948 mpd_clear_flags(result);
2949 result->exp = 0;
2950 result->len = _mpd_real_size(result->data, big->len);
2951 mpd_qresize(result, result->len, status);
2952 mpd_setdigits(result);
2953 _mpd_cap(result, ctx);
2957 mpd_seterror(result, MPD_Invalid_operation, status);
2965 mpd_qrotate(mpd_t *result, const mpd_t *a, const mpd_t *b,
2975 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
2980 mpd_seterror(result, MPD_Invalid_operation, status);
2986 mpd_seterror(result, MPD_Invalid_operation, status);
2990 mpd_seterror(result, MPD_Invalid_operation, status);
2994 mpd_qcopy(result, a, status);
3009 mpd_seterror(result, MPD_Malloc_error, status);
3017 mpd_seterror(result, MPD_Malloc_error, status);
3023 mpd_seterror(result, MPD_Malloc_error, status);
3026 _mpd_qadd(result, &big, &small, ctx, status);
3038 * The result is a with the value of b added to its exponent.
3041 mpd_qscaleb(mpd_t *result, const mpd_t *a, const mpd_t *b,
3055 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3060 mpd_seterror(result, MPD_Invalid_operation, status);
3069 mpd_seterror(result, MPD_Invalid_operation, status);
3073 mpd_qcopy(result, a, status);
3099 mpd_qcopy(result, a, status);
3100 result->exp = (mpd_ssize_t)exp;
3102 mpd_qfinalize(result, ctx, status);
3107 * of a left shift, the result is decapitated to fit the context precision. If
3111 mpd_qshiftn(mpd_t *result, const mpd_t *a, mpd_ssize_t n, const mpd_context_t *ctx,
3115 if (mpd_qcheck_nan(result, a, ctx, status)) {
3118 mpd_qcopy(result, a, status);
3123 mpd_qshiftl(result, a, n, status);
3124 _mpd_cap(result, ctx);
3127 if (!mpd_qcopy(result, a, status)) {
3130 _mpd_cap(result, ctx);
3131 mpd_qshiftr_inplace(result, -n);
3134 mpd_seterror(result, MPD_Invalid_operation, status);
3143 mpd_qshift(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx,
3150 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3155 mpd_seterror(result, MPD_Invalid_operation, status);
3161 mpd_seterror(result, MPD_Invalid_operation, status);
3165 mpd_seterror(result, MPD_Invalid_operation, status);
3169 mpd_qcopy(result, a, status);
3174 mpd_qshiftl(result, a, n, status);
3175 _mpd_cap(result, ctx);
3178 if (!mpd_qcopy(result, a, status)) {
3181 _mpd_cap(result, ctx);
3182 mpd_qshiftr_inplace(result, -n);
3188 mpd_qxor(mpd_t *result, const mpd_t *a, const mpd_t *b,
3199 mpd_seterror(result, MPD_Invalid_operation, status);
3206 if (!mpd_qresize(result, big->len, status)) {
3226 result->data[i] = z;
3253 result->data[i++] = z;
3264 result->data[i] = big->data[i];
3267 mpd_clear_flags(result);
3268 result->exp = 0;
3269 result->len = _mpd_real_size(result->data, big->len);
3270 mpd_qresize(result, result->len, status);
3271 mpd_setdigits(result);
3272 _mpd_cap(result, ctx);
3276 mpd_seterror(result, MPD_Invalid_operation, status);
3285 * The absolute value of a. If a is negative, the result is the same
3286 * as the result of the minus operation. Otherwise, the result is the
3287 * result of the plus operation.
3290 mpd_qabs(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
3294 if (mpd_qcheck_nan(result, a, ctx, status)) {
3300 mpd_qminus(result, a, ctx, status);
3303 mpd_qplus(result, a, ctx, status);
3317 _mpd_qaddsub_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
3322 mpd_seterror(result, MPD_Invalid_operation, status);
3325 mpd_setspecial(result, mpd_sign(a), MPD_INF);
3330 mpd_setspecial(result, sign_b, MPD_INF);
3335 _mpd_qaddsub(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
3401 mpd_seterror(result, MPD_Malloc_error, status);
3407 result->exp = small->exp;
3417 if (!mpd_qresize(result, newsize, status)) {
3423 carry = _mpd_baseadd(result->data, big->data, small->data,
3428 if (!mpd_qresize(result, newsize, status)) {
3431 result->data[newsize-1] = carry;
3434 result->len = newsize;
3435 mpd_set_flags(result, sign_b);
3450 _mpd_basesub(result->data, big->data, small->data,
3452 newsize = _mpd_real_size(result->data, big->len);
3454 (void)mpd_qresize(result, newsize, status);
3456 result->len = newsize;
3458 mpd_set_flags(result, sign_b);
3460 if (mpd_iszerocoeff(result)) {
3461 mpd_set_positive(result);
3463 mpd_set_negative(result);
3468 mpd_setdigits(result);
3476 _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
3479 _mpd_qaddsub(result, a, b, mpd_sign(b), ctx, status);
3484 _mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b,
3487 _mpd_qaddsub(result, a, b, !mpd_sign(b), ctx, status);
3492 mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
3496 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3499 _mpd_qaddsub_inf(result, a, b, mpd_sign(b), status);
3503 _mpd_qaddsub(result, a, b, mpd_sign(b), ctx, status);
3504 mpd_qfinalize(result, ctx, status);
3507 /* Add a and b. Set NaN/Invalid_operation if the result is inexact. */
3509 _mpd_qadd_exact(mpd_t *result, const mpd_t *a, const mpd_t *b,
3514 mpd_qadd(result, a, b, ctx, &workstatus);
3517 mpd_seterror(result, MPD_Invalid_operation, status);
3523 mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b,
3527 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
3530 _mpd_qaddsub_inf(result, a, b, !mpd_sign(b), status);
3534 _mpd_qaddsub(result, a, b, !mpd_sign(b), ctx, status);
3535 mpd_qfinalize(result, ctx, status);
3538 /* Subtract b from a. Set NaN/Invalid_operation if the result is inexact. */
3540 _mpd_qsub_exact(mpd_t *result, const mpd_t *a, const mpd_t *b,
3545 mpd_qsub(result, a, b, ctx, &workstatus);
3548 mpd_seterror(result, MPD_Invalid_operation, status);
3554 mpd_qadd_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
3562 mpd_qadd(result, a, &bb, ctx, status);
3568 mpd_qadd_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
3576 mpd_qadd(result, a, &bb, ctx, status);
3582 mpd_qsub_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
3590 mpd_qsub(result, a, &bb, ctx, status);
3596 mpd_qsub_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
3604 mpd_qsub(result, a, &bb, ctx, status);
3610 mpd_qadd_i32(mpd_t *result, const mpd_t *a, int32_t b,
3613 mpd_qadd_ssize(result, a, b, ctx, status);
3618 mpd_qadd_u32(mpd_t *result, const mpd_t *a, uint32_t b,
3621 mpd_qadd_uint(result, a, b, ctx, status);
3627 mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b,
3630 mpd_qadd_ssize(result, a, b, ctx, status);
3635 mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3638 mpd_qadd_uint(result, a, b, ctx, status);
3643 mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b,
3651 mpd_qadd(result, a, &bb, ctx, status);
3657 mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3665 mpd_qadd(result, a, &bb, ctx, status);
3672 mpd_qsub_i32(mpd_t *result, const mpd_t *a, int32_t b,
3675 mpd_qsub_ssize(result, a, b, ctx, status);
3680 mpd_qsub_u32(mpd_t *result, const mpd_t *a, uint32_t b,
3683 mpd_qsub_uint(result, a, b, ctx, status);
3689 mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b,
3692 mpd_qsub_ssize(result, a, b, ctx, status);
3697 mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3700 mpd_qsub_uint(result, a, b, ctx, status);
3705 mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b,
3713 mpd_qsub(result, a, &bb, ctx, status);
3719 mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b,
3727 mpd_qsub(result, a, &bb, ctx, status);
3735 _mpd_qdiv_inf(mpd_t *result, const mpd_t *a, const mpd_t *b,
3740 mpd_seterror(result, MPD_Invalid_operation, status);
3743 mpd_setspecial(result, mpd_sign(a)^mpd_sign(b), MPD_INF);
3747 _settriple(result, mpd_sign(a)^mpd_sign(b), 0, mpd_etiny(ctx));
3904 * the operation with a lower precision in case the result is exact.
3907 * when the result is exact. If a_coeff' * 1 / b_coeff' is in lowest
4189 mpd_qdiv_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
4197 mpd_qdiv(result, a, &bb, ctx, status);
4203 mpd_qdiv_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
4211 mpd_qdiv(result, a, &bb, ctx, status);
4217 mpd_qdiv_i32(mpd_t *result, const mpd_t *a, int32_t b,
4220 mpd_qdiv_ssize(result, a, b, ctx, status);
4225 mpd_qdiv_u32(mpd_t *result, const mpd_t *a, uint32_t b,
4228 mpd_qdiv_uint(result, a, b, ctx, status);
4234 mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b,
4237 mpd_qdiv_ssize(result, a, b, ctx, status);
4242 mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b,
4245 mpd_qdiv_uint(result, a, b, ctx, status);
4250 mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b,
4258 mpd_qdiv(result, a, &bb, ctx, status);
4264 mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b,
4272 mpd_qdiv(result, a, &bb, ctx, status);
4277 /* Pad the result with trailing zeros if it has fewer digits than prec. */
4279 _mpd_zeropad(mpd_t *result, const mpd_context_t *ctx, uint32_t *status)
4281 if (!mpd_isspecial(result) && !mpd_iszero(result) &&
4282 result->digits < ctx->prec) {
4283 mpd_ssize_t shift = ctx->prec - result->digits;
4284 mpd_qshiftl(result, result, shift, status);
4285 result->exp -= shift;
4289 /* Check if the result is guaranteed to be one. */
4291 _mpd_qexp_check_one(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4301 _settriple(result, 0, 1, 0);
4353 * and Underflow, two cases must be considered for the error of the result:
4355 * 1) abs(a) <= 9 * 10**(-prec-1) ==> result == 1
4362 * Relative error: abs(result - e**x) < 0.5 * 10**(-prec) * e**x
4382 _mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4394 _settriple(result, MPD_POS, 1, 0);
4420 mpd_setspecial(result, MPD_POS, MPD_INF);
4424 _settriple(result, MPD_POS, 0, mpd_etiny(ctx));
4432 if (_mpd_qexp_check_one(result, a, ctx, status)) {
4441 if (!mpd_qcopy(result, a, status)) {
4444 result->exp -= t;
4453 n = _mpd_get_exp_iterations(result, workctx.prec);
4455 mpd_seterror(result, MPD_Invalid_operation, status); /* GCOV_UNLIKELY */
4464 mpd_qdiv(&tmp, result, &word, &workctx, &workctx.status);
4469 _mpd_qpow_uint(result, &sum, mpd_pow10[t], MPD_POS, &workctx, status);
4472 _mpd_qpow_uint(result, &sum, mpd_pow10[t], MPD_POS, &workctx, status);
4478 _mpd_qpow_uint(result, &tmp, mpd_pow10[t], MPD_POS, &workctx, status);
4490 mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4496 if (mpd_qcheck_nan(result, a, ctx, status)) {
4500 _settriple(result, MPD_POS, 0, 0);
4503 mpd_setspecial(result, MPD_POS, MPD_INF);
4508 _settriple(result, MPD_POS, 1, 0);
4524 if (result == a) {
4526 mpd_seterror(result, MPD_Malloc_error, status);
4538 _mpd_qexp(result, a, &workctx, &workstatus);
4541 ulpexp = result->exp + result->digits - workctx.prec;
4543 /* The effective work precision is result->digits. */
4544 ulpexp = result->exp;
4550 * 1) abs(result - e**x) < 0.5 * 10**(-prec) * e**x
4551 * 2) result - ulp < e**x < result + ulp
4552 * 3) result - ulp < result < result + ulp
4554 * If round(result-ulp)==round(result+ulp), then
4555 * round(result)==round(e**x). Therefore the result
4562 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
4563 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
4564 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
4567 _mpd_zeropad(result, &workctx, status);
4568 mpd_check_underflow(result, &workctx, status);
4569 mpd_qfinalize(result, &workctx, status);
4580 _mpd_qexp(result, a, &workctx, status);
4581 _mpd_zeropad(result, &workctx, status);
4582 mpd_check_underflow(result, &workctx, status);
4583 mpd_qfinalize(result, &workctx, status);
4589 mpd_qfma(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c,
4595 if (result == c) {
4597 mpd_seterror(result, MPD_Malloc_error, status);
4603 _mpd_qmul(result, a, b, ctx, &workstatus);
4605 mpd_qadd(result, result, c, ctx, &workstatus);
4622 * 2) abs(log(v) - result) < 10**(-2*k_n-1 + 1) <= 10**-maxprec.
4699 * Set 'result' to log(10).
4700 * Ulp error: abs(result - log(10)) < ulp(log(10))
4701 * Relative error: abs(result - log(10)) < 5 * 10**-prec * log(10)
4707 mpd_qln10(mpd_t *result, mpd_ssize_t prec, uint32_t *status)
4722 rnd = mpd_qshiftr(result, &_mpd_ln10, shift, status);
4724 mpd_seterror(result, MPD_Malloc_error, status);
4727 result->exp = -(result->digits-1);
4732 _mpd_apply_round_excess(result, rnd, &maxcontext, status);
4740 i = ln_schedule_prec(klist, prec+2, -result->exp);
4743 result->flags ^= MPD_NEG;
4744 _mpd_qexp(&tmp, result, &varcontext, status);
4745 result->flags ^= MPD_NEG;
4748 mpd_qadd(result, result, &tmp, &maxcontext, status);
4749 if (mpd_isspecial(result)) {
4756 mpd_qfinalize(result, &maxcontext, status);
4839 * Relative error: abs(result - log(a)) < 0.1 * 10**-prec * abs(log(a))
4842 _mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
4846 mpd_t *z = result;
4863 mpd_seterror(result, MPD_Malloc_error, status);
4921 mpd_seterror(result, MPD_Malloc_error, status);
4931 _settriple(result, (cmp<0), 1, mpd_etiny(ctx)-1);
4939 * the result is (using 10**adjexp(x) <= abs(x)):
4974 * t * log(10) == 0, the result does not change and the analysis
4998 mpd_qadd(result, &tmp, z, &maxcontext, status);
5010 mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
5017 if (mpd_qcheck_nan(result, a, ctx, status)) {
5021 mpd_seterror(result, MPD_Invalid_operation, status);
5024 mpd_setspecial(result, MPD_POS, MPD_INF);
5028 mpd_setspecial(result, MPD_NEG, MPD_INF);
5032 mpd_seterror(result, MPD_Invalid_operation, status);
5036 _settriple(result, MPD_POS, 0, 0);
5040 * Check if the result will overflow (0 < x, x != 1):
5065 mpd_setspecial(result, (adjexp<0), MPD_INF);
5079 if (result == a) {
5081 mpd_seterror(result, MPD_Malloc_error, status);
5091 _mpd_qln(result, a, &workctx, status);
5093 result->exp + result->digits-workctx.prec);
5096 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
5097 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
5098 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
5101 mpd_check_underflow(result, &workctx, status);
5102 mpd_qfinalize(result, &workctx, status);
5113 _mpd_qln(result, a, &workctx, status);
5114 mpd_check_underflow(result, &workctx, status);
5115 mpd_qfinalize(result, &workctx, status);
5122 * Relative error: abs(result - log10(a)) < 0.1 * 10**-prec * abs(log10(a))
5124 * Ulp error: abs(result - log10(a)) < ulp(log10(a))
5128 _mpd_qlog10(int action, mpd_t *result, const mpd_t *a,
5137 * in _mpd_qln() does not change the final result. */
5138 _mpd_qln(result, a, &workctx, status);
5147 _mpd_qdiv(NO_IDEAL_EXP, result, result, &ln10, &workctx, status);
5154 mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
5164 if (mpd_qcheck_nan(result, a, ctx, status)) {
5168 mpd_seterror(result, MPD_Invalid_operation, status);
5171 mpd_setspecial(result, MPD_POS, MPD_INF);
5175 mpd_setspecial(result, MPD_NEG, MPD_INF);
5179 mpd_seterror(result, MPD_Invalid_operation, status);
5189 _settriple(result, sign, adjexp, 0);
5190 mpd_qfinalize(result, &workctx, status);
5194 * Check if the result will overflow (0 < x, x != 1):
5217 mpd_setspecial(result, (adjexp<0), MPD_INF);
5228 if (result == a) {
5230 mpd_seterror(result, MPD_Malloc_error, status);
5240 _mpd_qlog10(SKIP_FINALIZE, result, a, &workctx, status);
5242 result->exp + result->digits-workctx.prec);
5245 mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
5246 mpd_qsub(&t2, result, &ulp, &workctx, &workctx.status);
5247 if (mpd_isspecial(result) || mpd_iszerocoeff(result) ||
5250 mpd_check_underflow(result, &workctx, status);
5251 mpd_qfinalize(result, &workctx, status);
5262 _mpd_qlog10(DO_FINALIZE, result, a, &workctx, status);
5263 mpd_check_underflow(result, &workctx, status);
5273 mpd_qmax(mpd_t *result, const mpd_t *a, const mpd_t *b,
5279 mpd_qcopy(result, b, status);
5282 mpd_qcopy(result, a, status);
5284 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
5294 mpd_qcopy(result, b, status);
5297 mpd_qcopy(result, a, status);
5301 mpd_qfinalize(result, ctx, status);
5309 mpd_qmax_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
5315 mpd_qcopy(result, b, status);
5318 mpd_qcopy(result, a, status);
5320 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
5330 mpd_qcopy(result, b, status);
5333 mpd_qcopy(result, a, status);
5337 mpd_qfinalize(result, ctx, status);
5346 mpd_qmin(mpd_t *result, const mpd_t *a, const mpd_t *b,
5352 mpd_qcopy(result, b, status);
5355 mpd_qcopy(result, a, status);
5357 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
5367 mpd_qcopy(result, a, status);
5370 mpd_qcopy(result, b, status);
5374 mpd_qfinalize(result, ctx, status);
5382 mpd_qmin_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
5388 mpd_qcopy(result, b, status);
5391 mpd_qcopy(result, a, status);
5393 else if (mpd_qcheck_nans(result, a, b, ctx, status)) {
5403 mpd_qcopy(result, a, status);
5406 mpd_qcopy(result, b, status);
5410 mpd_qfinalize(result, ctx, status);
5413 /* Minimum space needed for the result array in _karatsuba_rec(). */
5475 lt = lb + lb + 1; /* space needed for result array */
5476 mpd_uint_zero(w, lt); /* clear result array */
5480 lt = (la-m) + (la-m) + 1; /* space needed for result array */
5481 mpd_uint_zero(w, lt); /* clear result array */
5486 lt = m + m + 1; /* space needed for the result array */
5487 mpd_uint_zero(w, lt); /* clear result array */
5525 * to the result or NULL in case of failure (malloc error).
5533 mpd_uint_t *result = NULL, *w = NULL;
5540 if ((result = mpd_calloc(*rsize, sizeof *result)) == NULL) {
5546 mpd_free(result);
5550 _karatsuba_rec(result, u, v, w, ulen, vlen);
5554 return result;
5649 * a pointer to the result or NULL in case of failure (malloc error).
5753 mpd_uint_t *result;
5756 if ((result = _mpd_fntmul(a, b, la, lb, &dummy)) == NULL) {
5759 memcpy(c, result, (la+lb) * (sizeof *result));
5760 mpd_free(result);
5772 lt = lb + lb + 1; /* space needed for result array */
5773 mpd_uint_zero(w, lt); /* clear result array */
5779 lt = (la-m) + (la-m) + 1; /* space needed for result array */
5780 mpd_uint_zero(w, lt); /* clear result array */
5787 lt = m + m + 1; /* space needed for the result array */
5788 mpd_uint_zero(w, lt); /* clear result array */
5834 * base case. Returns a pointer to the result or NULL in case of failure
5842 mpd_uint_t *result = NULL, *w = NULL;
5849 if ((result = mpd_calloc(*rsize, sizeof *result)) == NULL) {
5855 mpd_free(result); /* GCOV_UNLIKELY */
5859 if (!_karatsuba_rec_fnt(result, u, v, w, ulen, vlen)) {
5860 mpd_free(result);
5861 result = NULL;
5866 return result;
5872 _mpd_qmul_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, uint32_t *status)
5876 mpd_seterror(result, MPD_Invalid_operation, status);
5879 mpd_setspecial(result, mpd_sign(a)^mpd_sign(b), MPD_INF);
5885 mpd_seterror(result, MPD_Invalid_operation, status);
5888 mpd_setspecial(result, mpd_sign(a)^mpd_sign(b), MPD_INF);
5894 * does NOT finalize the result. This is for use in mpd_fma().
5897 _mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b,
5907 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
5910 _mpd_qmul_inf(result, a, b, status);
5921 _mpd_singlemul(result->data, big->data[0], small->data[0]);
5937 if (!mpd_qresize(result, rsize, status)) {
5941 result->data[i] = rbuf[i];
5969 mpd_seterror(result, MPD_Malloc_error, status);
5973 if (mpd_isdynamic_data(result)) {
5974 mpd_free(result->data);
5976 result->data = rdata;
5977 result->alloc = rsize;
5978 mpd_set_dynamic_data(result);
5982 mpd_set_flags(result, mpd_sign(a)^mpd_sign(b));
5983 result->exp = big->exp + small->exp;
5984 result->len = _mpd_real_size(result->data, rsize);
5986 mpd_qresize(result, result->len, status);
5987 mpd_setdigits(result);
5992 mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b,
5995 _mpd_qmul(result, a, b, ctx, status);
5996 mpd_qfinalize(result, ctx, status);
5999 /* Multiply a and b. Set NaN/Invalid_operation if the result is inexact. */
6001 _mpd_qmul_exact(mpd_t *result, const mpd_t *a, const mpd_t *b,
6006 mpd_qmul(result, a, b, ctx, &workstatus);
6009 mpd_seterror(result, MPD_Invalid_operation, status);
6015 mpd_qmul_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
6023 mpd_qmul(result, a, &bb, ctx, status);
6029 mpd_qmul_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
6037 mpd_qmul(result, a, &bb, ctx, status);
6042 mpd_qmul_i32(mpd_t *result, const mpd_t *a, int32_t b,
6045 mpd_qmul_ssize(result, a, b, ctx, status);
6049 mpd_qmul_u32(mpd_t *result, const mpd_t *a, uint32_t b,
6052 mpd_qmul_uint(result, a, b, ctx, status);
6057 mpd_qmul_i64(mpd_t *result, const mpd_t *a, int64_t b,
6060 mpd_qmul_ssize(result, a, b, ctx, status);
6064 mpd_qmul_u64(mpd_t *result, const mpd_t *a, uint64_t b,
6067 mpd_qmul_uint(result, a, b, ctx, status);
6072 mpd_qmul_i64(mpd_t *result, const mpd_t *a, int64_t b,
6080 mpd_qmul(result, a, &bb, ctx, status);
6086 mpd_qmul_u64(mpd_t *result, const mpd_t *a, uint64_t b,
6094 mpd_qmul(result, a, &bb, ctx, status);
6101 mpd_qminus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
6105 if (mpd_qcheck_nan(result, a, ctx, status)) {
6111 mpd_qcopy_abs(result, a, status);
6114 mpd_qcopy_negate(result, a, status);
6117 mpd_qfinalize(result, ctx, status);
6122 mpd_qplus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
6126 if (mpd_qcheck_nan(result, a, ctx, status)) {
6132 mpd_qcopy_abs(result, a, status);
6135 mpd_qcopy(result, a, status);
6138 mpd_qfinalize(result, ctx, status);
6143 mpd_qnext_minus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
6150 if (mpd_qcheck_nan(result, a, ctx, status)) {
6156 mpd_qcopy(result, a, status);
6160 mpd_clear_flags(result);
6161 mpd_qmaxcoeff(result, ctx, status);
6162 if (mpd_isnan(result)) {
6165 result->exp = mpd_etop(ctx);
6173 if (!mpd_qcopy(result, a, status)) {
6177 mpd_qfinalize(result, &workctx, &workctx.status);
6184 mpd_qsub(result, a, &tiny, &workctx, &workctx.status);
6190 mpd_qnext_plus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
6197 if (mpd_qcheck_nan(result, a, ctx, status)) {
6203 mpd_qcopy(result, a, status);
6206 mpd_clear_flags(result);
6207 mpd_qmaxcoeff(result, ctx, status);
6208 if (mpd_isnan(result)) {
6211 mpd_set_flags(result, MPD_NEG);
6212 result->exp = mpd_etop(ctx);
6220 if (!mpd_qcopy(result, a, status)) {
6224 mpd_qfinalize(result, &workctx, &workctx.status);
6231 mpd_qadd(result, a, &tiny, &workctx, &workctx.status);
6240 mpd_qnext_toward(mpd_t *result, const mpd_t *a, const mpd_t *b,
6245 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
6251 mpd_qcopy_sign(result, a, b, status);
6256 mpd_qnext_plus(result, a, ctx, status);
6259 mpd_qnext_minus(result, a, ctx, status);
6262 if (mpd_isinfinite(result)) {
6265 else if (mpd_adjexp(result) < ctx->emin) {
6267 if (mpd_iszero(result)) {
6281 * result = x**k * (1 + err)**(k-1)
6284 _mpd_qpow_uint(mpd_t *result, const mpd_t *base, mpd_uint_t exp,
6291 _settriple(result, resultsign, 1, 0); /* GCOV_NOT_REACHED */
6295 if (!mpd_qcopy(result, base, status)) {
6301 mpd_qmul(result, result, result, ctx, &workstatus);
6303 mpd_qmul(result, result, base, ctx, &workstatus);
6305 if (mpd_isspecial(result) ||
6306 (mpd_iszerocoeff(result) && (workstatus & MPD_Clamped))) {
6312 mpd_set_sign(result, resultsign);
6323 * result = x**k * (1 + err)**k
6326 _mpd_qpow_mpd(mpd_t *result, mpd_t *tbase, mpd_t *texp, uint8_t resultsign,
6337 mpd_qcopy(result, &one, status);
6341 mpd_qmul(result, result, tbase, ctx, &workstatus);
6343 if (mpd_isspecial(result) ||
6344 (mpd_iszerocoeff(result) && (workstatus & MPD_Clamped))) {
6351 mpd_seterror(result, workstatus&MPD_Errors, status);
6355 mpd_set_sign(result, resultsign);
6361 * abs(result - base**exp) < 0.1 * 10**-prec * abs(base**exp)
6364 _mpd_qpow_int(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6384 mpd_setspecial(result, MPD_POS, MPD_NAN);
6390 mpd_setspecial(result, MPD_POS, MPD_NAN);
6398 mpd_setspecial(result, MPD_POS, MPD_NAN); /* GCOV_UNLIKELY */
6401 _mpd_qpow_mpd(result, &tbase, &texp, resultsign, &workctx, status);
6404 _mpd_qpow_uint(result, &tbase, n, resultsign, &workctx, status);
6407 if (mpd_isinfinite(result)) {
6409 _settriple(result, resultsign, 1, MPD_EXP_INF);
6415 mpd_qfinalize(result, ctx, status);
6419 * If the exponent is infinite and base equals one, the result is one
6420 * with a coefficient of length prec. Otherwise, result is undefined.
6424 _qcheck_pow_one_inf(mpd_t *result, const mpd_t *base, uint8_t resultsign,
6432 mpd_qshiftl(result, &one, shift, status);
6433 result->exp = -shift;
6434 mpd_set_flags(result, resultsign);
6442 * If abs(base) equals one, calculate the correct power of one result.
6443 * Otherwise, result is undefined. Return the value of the comparison
6449 _qcheck_pow_one(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6460 _settriple(result, resultsign, 1, 0);
6464 mpd_qmul_ssize(result, exp, -base->exp, ctx, &workstatus);
6470 shift = mpd_qget_ssize(result, &workstatus);
6471 /* shift is MPD_SSIZE_MAX if result is too large */
6484 if (!mpd_qshiftl(result, &one, shift, status)) {
6487 result->exp = -shift;
6488 mpd_set_flags(result, resultsign);
6525 * is strictly increasing, the end result is a lower bound.
6597 _qcheck_pow_bounds(mpd_t *result, const mpd_t *x, const mpd_t *y,
6610 mpd_seterror(result, MPD_Malloc_error, status);
6619 _settriple(result, resultsign, 1, MPD_EXP_INF);
6620 mpd_qfinalize(result, ctx, status);
6628 _settriple(result, resultsign, 1, mpd_etiny(ctx)-1);
6629 mpd_qfinalize(result, ctx, status);
6644 _mpd_qpow_exact(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6653 * Relative error: abs(result - e**y) < e**y * 1/5 * 10**(-prec - 1)
6656 _mpd_qpow_real(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6663 mpd_seterror(result, MPD_Malloc_error, status);
6695 * 7) abs(result - e**y) < e**y * 1/5 * 10**(-prec - 1)
6697 mpd_qln(result, base, &workctx, &workctx.status);
6698 mpd_qmul(result, result, &texp, &workctx, &workctx.status);
6699 mpd_qexp(result, result, &workctx, status);
6708 mpd_qpow(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6716 if (mpd_qcheck_nans(result, base, exp, ctx, status)) {
6727 mpd_seterror(result, MPD_Invalid_operation, status);
6730 mpd_setspecial(result, resultsign, MPD_INF);
6733 _settriple(result, resultsign, 0, 0);
6739 mpd_seterror(result, MPD_Invalid_operation, status);
6745 cmp = _qcheck_pow_one_inf(result, base, resultsign, ctx, status);
6752 _settriple(result, resultsign, 0, 0);
6755 mpd_setspecial(result, resultsign, MPD_INF);
6762 _settriple(result, resultsign, 1, 0);
6765 _settriple(result, resultsign, 0, 0);
6768 mpd_setspecial(result, resultsign, MPD_INF);
6773 _settriple(result, resultsign, 1, 0);
6776 if (_qcheck_pow_one(result, base, exp, resultsign, ctx, status) == 0) {
6779 if (_qcheck_pow_bounds(result, base, exp, resultsign, ctx, status)) {
6784 _mpd_qpow_int(result, base, exp, resultsign, ctx, status);
6787 _mpd_qpow_real(result, base, exp, ctx, status);
6788 if (!mpd_isspecial(result) && _mpd_cmp(result, &one) == 0) {
6790 mpd_qshiftl(result, &one, shift, status);
6791 result->exp = -shift;
6793 if (mpd_isinfinite(result)) {
6795 _settriple(result, MPD_POS, 1, MPD_EXP_INF);
6797 mpd_qfinalize(result, ctx, status);
6806 _mpd_qpowmod_uint(mpd_t *result, mpd_t *base, mpd_uint_t exp,
6814 mpd_qcopy(result, &one, status);
6818 _mpd_qmul_exact(result, result, base, &maxcontext, status);
6819 mpd_qrem(result, result, mod, &maxcontext, status);
6829 mpd_qpowmod(mpd_t *result, const mpd_t *base, const mpd_t *exp,
6847 if (mpd_qcheck_3nans(result, base, exp, mod, ctx, status)) {
6850 mpd_seterror(result, MPD_Invalid_operation, status);
6856 mpd_seterror(result, MPD_Invalid_operation, status);
6860 mpd_seterror(result, MPD_Invalid_operation, status);
6864 mpd_seterror(result, MPD_Invalid_operation, status);
6871 mpd_seterror(result, MPD_Invalid_operation, status);
6875 _settriple(result, sign, r, 0);
6879 mpd_seterror(result, MPD_Invalid_operation, status);
6883 _settriple(result, sign, 0, 0);
6891 mpd_seterror(result, maxcontext.status&MPD_Errors, status);
6908 mpd_qshiftl(result, &one, tbase_exp, status);
6909 mpd_qrem(result, result, &tmod, &maxcontext, status);
6910 _mpd_qmul_exact(&tbase, &tbase, result, &maxcontext, status);
6929 mpd_qcopy(result, &one, status);
6932 _mpd_qmul_exact(result, result, &tbase, &maxcontext, status);
6933 mpd_qrem(result, result, &tmod, &maxcontext, status);
6940 mpd_isspecial(&tmod) || mpd_isspecial(result)) {
6945 mpd_set_sign(result, sign);
6956 mpd_setspecial(result, MPD_POS, MPD_NAN);
6961 mpd_qquantize(mpd_t *result, const mpd_t *a, const mpd_t *b,
6970 if (mpd_qcheck_nans(result, a, b, ctx, status)) {
6974 mpd_qcopy(result, a, status);
6977 mpd_seterror(result, MPD_Invalid_operation, status);
6982 mpd_seterror(result, MPD_Invalid_operation, status);
6987 _settriple(result, mpd_sign(a), 0, b->exp);
6988 mpd_qfinalize(result, ctx, status);
6995 mpd_seterror(result, MPD_Invalid_operation, status);
7001 if (!mpd_qshiftl(result, a, shift, status)) {
7004 result->exp = b_exp;
7010 rnd = mpd_qshiftr(result, a, shift, status);
7014 result->exp = b_exp;
7015 if (!_mpd_apply_round_fit(result, rnd, ctx, status)) {
7024 if (mpd_adjexp(result) > ctx->emax ||
7025 mpd_adjexp(result) < mpd_etiny(ctx)) {
7026 mpd_seterror(result, MPD_Invalid_operation, status);
7031 mpd_qfinalize(result, ctx, status);
7035 mpd_qreduce(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7042 if (mpd_qcheck_nan(result, a, ctx, status)) {
7045 mpd_qcopy(result, a, status);
7049 if (!mpd_qcopy(result, a, status)) {
7052 mpd_qfinalize(result, ctx, status);
7053 if (mpd_isspecial(result)) {
7056 if (mpd_iszero(result)) {
7057 _settriple(result, sign_a, 0, 0);
7061 shift = mpd_trail_zeros(result);
7063 /* After the finalizing above result->exp <= maxexp. */
7064 maxshift = maxexp - result->exp;
7067 mpd_qshiftr_inplace(result, shift);
7068 result->exp += shift;
7204 _mpd_qrescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp,
7211 mpd_qcopy(result, a, status);
7216 _settriple(result, mpd_sign(a), 0, exp);
7224 mpd_seterror(result, MPD_Invalid_operation, status);
7227 if (!mpd_qshiftl(result, a, shift, status)) {
7230 result->exp = exp;
7234 rnd = mpd_qshiftr(result, a, shift, status);
7238 result->exp = exp;
7239 _mpd_apply_round_excess(result, rnd, ctx, status);
7246 if (mpd_issubnormal(result, ctx)) {
7257 * result->digits <= MPD_MAX_PREC+1
7260 mpd_qrescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp,
7264 mpd_seterror(result, MPD_Invalid_operation, status);
7268 _mpd_qrescale(result, a, exp, ctx, status);
7272 * Same as mpd_qrescale, but with relaxed restrictions. The result of this
7277 * result->digits <= MPD_MAX_PREC+1
7280 mpd_qrescale_fmt(mpd_t *result, const mpd_t *a, mpd_ssize_t exp,
7284 mpd_seterror(result, MPD_Invalid_operation, status);
7288 _mpd_qrescale(result, a, exp, ctx, status);
7294 _mpd_qround_to_integral(int action, mpd_t *result, const mpd_t *a,
7300 if (mpd_qcheck_nan(result, a, ctx, status)) {
7303 mpd_qcopy(result, a, status);
7307 mpd_qcopy(result, a, status);
7311 _settriple(result, mpd_sign(a), 0, 0);
7315 rnd = mpd_qshiftr(result, a, -a->exp, status);
7319 result->exp = 0;
7322 _mpd_apply_round_excess(result, rnd, ctx, status);
7333 mpd_qround_to_intx(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7336 (void)_mpd_qround_to_integral(TO_INT_EXACT, result, a, ctx, status);
7340 mpd_qround_to_int(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7343 (void)_mpd_qround_to_integral(TO_INT_SILENT, result, a, ctx, status);
7347 mpd_qtrunc(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7351 mpd_seterror(result, MPD_Invalid_operation, status);
7355 (void)_mpd_qround_to_integral(TO_INT_TRUNC, result, a, ctx, status);
7359 mpd_qfloor(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7365 mpd_seterror(result, MPD_Invalid_operation, status);
7370 (void)_mpd_qround_to_integral(TO_INT_SILENT, result, a,
7375 mpd_qceil(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7381 mpd_seterror(result, MPD_Invalid_operation, status);
7386 (void)_mpd_qround_to_integral(TO_INT_SILENT, result, a,
7452 * Reciprocal, calculated with Newton's Method. Assumption: result != a.
7458 _mpd_qreciprocal(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7462 mpd_t *z = result; /* current approximation */
7473 assert(result != a);
7497 /* Let s := z**2, exact result */
7511 /* Let s := 2*z, exact result */
7535 * inexact result is fixed by a small loop, using at most one iteration.
7610 /* Fix the result. At this point -b < r < 2*b, so the correction loop
7674 * 2) abs(sqrt(v) - result) < 10**(-2*k_n-1 + 2) <= 10**-maxprec.
7748 * Set 'result' to 1/sqrt(a).
7749 * Relative error: abs(result - 1/sqrt(a)) < 10**-prec * 1/sqrt(a)
7752 _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7757 mpd_t *z = result; /* current approximation */
7775 if (result == a) {
7777 mpd_seterror(result, MPD_Malloc_error, status);
7834 tz = mpd_trail_zeros(result);
7835 shift = ideal_exp - result->exp;
7838 mpd_qshiftr_inplace(result, shift);
7839 result->exp += shift;
7851 mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7857 if (mpd_qcheck_nan(result, a, ctx, status)) {
7861 mpd_seterror(result, MPD_Invalid_operation, status);
7865 _settriple(result, MPD_POS, 0, mpd_etiny(ctx));
7870 mpd_setspecial(result, mpd_sign(a), MPD_INF);
7875 mpd_seterror(result, MPD_Invalid_operation, status);
7882 _mpd_qinvroot(result, a, &workctx, status);
7883 mpd_qfinalize(result, ctx, status);
7889 _mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
7905 if (mpd_qcheck_nan(result, a, ctx, status)) {
7909 mpd_seterror(result, MPD_Invalid_operation, status);
7912 mpd_setspecial(result, MPD_POS, MPD_INF);
7916 _settriple(result, mpd_sign(a), 0, ideal_exp);
7917 mpd_qfinalize(result, ctx, status);
7921 mpd_seterror(result, MPD_Invalid_operation, status);
7956 /* find result = floor(sqrt(c)) using Newton's method */
7957 if (!mpd_qshiftl(result, &one, prec, status)) {
7962 _mpd_qdivmod(&q, &r, &c, result, &maxcontext, &maxcontext.status);
7963 if (mpd_isspecial(result) || mpd_isspecial(&q)) {
7964 mpd_seterror(result, maxcontext.status&MPD_Errors, status);
7967 if (_mpd_cmp(result, &q) <= 0) {
7970 _mpd_qadd_exact(result, result, &q, &maxcontext, &maxcontext.status);
7971 if (mpd_isspecial(result)) {
7972 mpd_seterror(result, maxcontext.status&MPD_Errors, status);
7975 _mpd_qdivmod(result, &r, result, &two, &maxcontext, &maxcontext.status);
7979 _mpd_qmul_exact(&r, result, result, &maxcontext, &maxcontext.status);
7981 mpd_seterror(result, maxcontext.status&MPD_Errors, status);
7989 mpd_qshiftr_inplace(result, shift);
7992 if (!mpd_qshiftl(result, result, -shift, status)) {
7999 int lsd = (int)mpd_lsd(result->data[0]);
8001 result->data[0] += 1;
8005 result->exp = ideal_exp;
8014 mpd_qfinalize(result, &maxcontext, status);
8018 mpd_seterror(result, MPD_Malloc_error, status);
8023 mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
8029 if (result == a) {
8031 mpd_seterror(result, MPD_Malloc_error, status);
8037 _mpd_qsqrt(result, a, ctx, &xstatus);
8042 * a lower context precision in case the result is exact.
8044 * If the result is exact, an upper bound for the number of digits
8058 _mpd_qsqrt(result, a, &workctx, &ystatus);
8061 mpd_seterror(result, ystatus, status);
8362 * The least significant word of the result is (*rdata)[0].
8441 * The least significant word of the result is (*rdata)[0].
8538 mpd_qimport_u16(mpd_t *result,
8544 mpd_ssize_t rlen; /* length of the result */
8552 mpd_seterror(result, MPD_Invalid_operation, status);
8558 mpd_seterror(result, MPD_Malloc_error, status);
8565 if (!mpd_qresize(result, rlen, status)) {
8569 n = _coeff_from_u16(result, rlen, usrc, srclen, srcbase, status);
8574 mpd_set_flags(result, srcsign);
8575 result->exp = 0;
8576 result->len = n;
8577 mpd_setdigits(result);
8579 mpd_qresize(result, result->len, status);
8580 mpd_qfinalize(result, ctx, status);
8592 mpd_qimport_u32(mpd_t *result,
8597 mpd_ssize_t rlen; /* length of the result */
8604 mpd_seterror(result, MPD_Invalid_operation, status);
8608 if (!mpd_qresize(result, rlen, status)) {
8613 n = _coeff_from_smaller_base(result, rlen, MPD_RADIX,
8618 if (!mpd_qresize(result, srclen, status)) {
8621 memcpy(result->data, srcdata, srclen * (sizeof *srcdata));
8625 n = _coeff_from_smaller_base(result, rlen, MPD_RADIX,
8632 mpd_seterror(result, MPD_Malloc_error, status);
8639 n = _coeff_from_larger_base(result, rlen, MPD_RADIX,
8650 mpd_set_flags(result, srcsign);
8651 result->exp = 0;
8652 result->len = n;
8653 mpd_setdigits(result);
8655 mpd_qresize(result, result->len, status);
8656 mpd_qfinalize(result, ctx, status);
8742 _set_uint128_coeff_exp(mpd_t *result, uint64_t hi, uint64_t lo, mpd_ssize_t exp)
8754 if (!mpd_qresize(result, len, &status)) {
8759 result->data[i] = data[i];
8762 result->exp = exp;
8763 result->len = len;
8764 mpd_setdigits(result);
8770 mpd_from_uint128_triple(mpd_t *result, const mpd_uint128_triple_t *triple, uint32_t *status)
8803 mpd_setspecial(result, sign, flags);
8809 if (_set_uint128_coeff_exp(result, hi, lo, exp) < 0) {
8821 mpd_setspecial(result, sign, MPD_INF);
8832 mpd_set_flags(result, flags);
8841 if (_set_uint128_coeff_exp(result, hi, lo, exp) < 0) {
8846 mpd_qfinalize(result, &maxcontext, &workstatus);
8859 mpd_seterror(result, MPD_Conversion_syntax, status);
8863 mpd_seterror(result, MPD_Malloc_error, status);