/* * Copyright (c) 2008-2020 Stefan Krah. All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. */ #include "mpdecimal.h" #include <assert.h> #include <limits.h> #include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include "basearith.h" #include "bits.h" #include "constants.h" #include "convolute.h" #include "crt.h" #include "mpalloc.h" #include "typearith.h" #ifdef PPRO #if defined(_MSC_VER) #include <float.h> #pragma float_control(precise, on) #pragma fenv_access(on) #elif !defined(__OpenBSD__) && !defined(__NetBSD__) /* C99 */ #include <fenv.h> #pragma STDC FENV_ACCESS ON #endif #endif /* Disable warning that is part of -Wextra since gcc 7.0. */ #if defined(__GNUC__) && !defined(__INTEL_COMPILER) && __GNUC__ >= 7 #pragma GCC diagnostic ignored "-Wimplicit-fallthrough" #endif #if defined(_MSC_VER) #define ALWAYS_INLINE … #elif defined (__IBMC__) || defined(LEGACY_COMPILER) #define ALWAYS_INLINE #undef inline #define inline #else #ifdef TEST_COVERAGE #define ALWAYS_INLINE #else #define ALWAYS_INLINE … #endif #endif /* ClangCL claims to support 128-bit int, but doesn't */ #if defined(__SIZEOF_INT128__) && defined(__clang__) && defined(_MSC_VER) #undef __SIZEOF_INT128__ #endif #define MPD_NEWTONDIV_CUTOFF … #define MPD_NEW_STATIC(name, flags, exp, digits, len) … #define MPD_NEW_CONST(name, flags, exp, digits, len, alloc, initval) … #define MPD_NEW_SHARED(name, a) … static mpd_uint_t data_one[1] = …; static mpd_uint_t data_zero[1] = …; static const mpd_t one = …; static const mpd_t minus_one = …; static const mpd_t zero = …; static inline void _mpd_check_exp(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status); static void _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp); static inline mpd_ssize_t _mpd_real_size(mpd_uint_t *data, mpd_ssize_t size); static int _mpd_cmp_abs(const mpd_t *a, const mpd_t *b); static void _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status); static inline void _mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status); static void _mpd_base_ndivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b, uint32_t *status); static inline void _mpd_qpow_uint(mpd_t *result, const mpd_t *base, mpd_uint_t exp, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status); static mpd_uint_t mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n); /******************************************************************************/ /* Version */ /******************************************************************************/ const char * mpd_version(void) { … } /******************************************************************************/ /* Performance critical inline functions */ /******************************************************************************/ #ifdef CONFIG_64 /* Digits in a word, primarily useful for the most significant word. */ ALWAYS_INLINE int mpd_word_digits(mpd_uint_t word) { … } #else ALWAYS_INLINE int mpd_word_digits(mpd_uint_t word) { if (word < mpd_pow10[4]) { if (word < mpd_pow10[2]) { return (word < mpd_pow10[1]) ? 1 : 2; } return (word < mpd_pow10[3]) ? 3 : 4; } if (word < mpd_pow10[6]) { return (word < mpd_pow10[5]) ? 5 : 6; } if (word < mpd_pow10[8]) { return (word < mpd_pow10[7]) ? 7 : 8; } return (word < mpd_pow10[9]) ? 9 : 10; } #endif /* Adjusted exponent */ ALWAYS_INLINE mpd_ssize_t mpd_adjexp(const mpd_t *dec) { … } /* Etiny */ ALWAYS_INLINE mpd_ssize_t mpd_etiny(const mpd_context_t *ctx) { … } /* Etop: used for folding down in IEEE clamping */ ALWAYS_INLINE mpd_ssize_t mpd_etop(const mpd_context_t *ctx) { … } /* Most significant word */ ALWAYS_INLINE mpd_uint_t mpd_msword(const mpd_t *dec) { … } /* Most significant digit of a word */ inline mpd_uint_t mpd_msd(mpd_uint_t word) { … } /* Least significant digit of a word */ ALWAYS_INLINE mpd_uint_t mpd_lsd(mpd_uint_t word) { … } /* Coefficient size needed to store 'digits' */ mpd_ssize_t mpd_digits_to_size(mpd_ssize_t digits) { … } /* Number of digits in the exponent. Not defined for MPD_SSIZE_MIN. */ inline int mpd_exp_digits(mpd_ssize_t exp) { … } /* Canonical */ ALWAYS_INLINE int mpd_iscanonical(const mpd_t *dec) { … } /* Finite */ ALWAYS_INLINE int mpd_isfinite(const mpd_t *dec) { … } /* Infinite */ ALWAYS_INLINE int mpd_isinfinite(const mpd_t *dec) { … } /* NaN */ ALWAYS_INLINE int mpd_isnan(const mpd_t *dec) { … } /* Negative */ ALWAYS_INLINE int mpd_isnegative(const mpd_t *dec) { … } /* Positive */ ALWAYS_INLINE int mpd_ispositive(const mpd_t *dec) { … } /* qNaN */ ALWAYS_INLINE int mpd_isqnan(const mpd_t *dec) { … } /* Signed */ ALWAYS_INLINE int mpd_issigned(const mpd_t *dec) { … } /* sNaN */ ALWAYS_INLINE int mpd_issnan(const mpd_t *dec) { … } /* Special */ ALWAYS_INLINE int mpd_isspecial(const mpd_t *dec) { … } /* Zero */ ALWAYS_INLINE int mpd_iszero(const mpd_t *dec) { … } /* Test for zero when specials have been ruled out already */ ALWAYS_INLINE int mpd_iszerocoeff(const mpd_t *dec) { … } /* Normal */ inline int mpd_isnormal(const mpd_t *dec, const mpd_context_t *ctx) { … } /* Subnormal */ inline int mpd_issubnormal(const mpd_t *dec, const mpd_context_t *ctx) { … } /* Odd word */ ALWAYS_INLINE int mpd_isoddword(mpd_uint_t word) { … } /* Odd coefficient */ ALWAYS_INLINE int mpd_isoddcoeff(const mpd_t *dec) { … } /* 0 if dec is positive, 1 if dec is negative */ ALWAYS_INLINE uint8_t mpd_sign(const mpd_t *dec) { … } /* 1 if dec is positive, -1 if dec is negative */ ALWAYS_INLINE int mpd_arith_sign(const mpd_t *dec) { … } /* Radix */ ALWAYS_INLINE long mpd_radix(void) { … } /* Dynamic decimal */ ALWAYS_INLINE int mpd_isdynamic(const mpd_t *dec) { … } /* Static decimal */ ALWAYS_INLINE int mpd_isstatic(const mpd_t *dec) { … } /* Data of decimal is dynamic */ ALWAYS_INLINE int mpd_isdynamic_data(const mpd_t *dec) { … } /* Data of decimal is static */ ALWAYS_INLINE int mpd_isstatic_data(const mpd_t *dec) { … } /* Data of decimal is shared */ ALWAYS_INLINE int mpd_isshared_data(const mpd_t *dec) { … } /* Data of decimal is const */ ALWAYS_INLINE int mpd_isconst_data(const mpd_t *dec) { … } /******************************************************************************/ /* Inline memory handling */ /******************************************************************************/ /* Fill destination with zeros */ ALWAYS_INLINE void mpd_uint_zero(mpd_uint_t *dest, mpd_size_t len) { … } /* Free a decimal */ ALWAYS_INLINE void mpd_del(mpd_t *dec) { … } /* * Resize the coefficient. Existing data up to 'nwords' is left untouched. * Return 1 on success, 0 otherwise. * * Input invariant: MPD_MINALLOC <= result->alloc. * * Case nwords == result->alloc: * 'result' is unchanged. Return 1. * * Case nwords > result->alloc: * Case realloc success: * The value of 'result' does not change. Return 1. * Case realloc failure: * 'result' is NaN, status is updated with MPD_Malloc_error. Return 0. * * Case nwords < result->alloc: * Case is_static_data or realloc failure [1]: * 'result' is unchanged. Return 1. * Case realloc success: * The value of result is undefined (expected). Return 1. * * * [1] In that case the old (now oversized) area is still valid. */ ALWAYS_INLINE int mpd_qresize(mpd_t *result, mpd_ssize_t nwords, uint32_t *status) { … } /* Same as mpd_qresize, but do not set the result no NaN on failure. */ static ALWAYS_INLINE int mpd_qresize_cxx(mpd_t *result, mpd_ssize_t nwords) { … } /* Same as mpd_qresize, but the complete coefficient (including the old * memory area!) is initialized to zero. */ ALWAYS_INLINE int mpd_qresize_zero(mpd_t *result, mpd_ssize_t nwords, uint32_t *status) { … } /* * Reduce memory size for the coefficient to MPD_MINALLOC. In theory, * realloc may fail even when reducing the memory size. But in that case * the old memory area is always big enough, so checking for MPD_Malloc_error * is not imperative. */ ALWAYS_INLINE void mpd_minalloc(mpd_t *result) { … } int mpd_resize(mpd_t *result, mpd_ssize_t nwords, mpd_context_t *ctx) { … } int mpd_resize_zero(mpd_t *result, mpd_ssize_t nwords, mpd_context_t *ctx) { … } /******************************************************************************/ /* Set attributes of a decimal */ /******************************************************************************/ /* Set digits. Assumption: result->len is initialized and > 0. */ inline void mpd_setdigits(mpd_t *result) { … } /* Set sign */ ALWAYS_INLINE void mpd_set_sign(mpd_t *result, uint8_t sign) { … } /* Copy sign from another decimal */ ALWAYS_INLINE void mpd_signcpy(mpd_t *result, const mpd_t *a) { … } /* Set infinity */ ALWAYS_INLINE void mpd_set_infinity(mpd_t *result) { … } /* Set qNaN */ ALWAYS_INLINE void mpd_set_qnan(mpd_t *result) { … } /* Set sNaN */ ALWAYS_INLINE void mpd_set_snan(mpd_t *result) { … } /* Set to negative */ ALWAYS_INLINE void mpd_set_negative(mpd_t *result) { … } /* Set to positive */ ALWAYS_INLINE void mpd_set_positive(mpd_t *result) { … } /* Set to dynamic */ ALWAYS_INLINE void mpd_set_dynamic(mpd_t *result) { … } /* Set to static */ ALWAYS_INLINE void mpd_set_static(mpd_t *result) { … } /* Set data to dynamic */ ALWAYS_INLINE void mpd_set_dynamic_data(mpd_t *result) { … } /* Set data to static */ ALWAYS_INLINE void mpd_set_static_data(mpd_t *result) { … } /* Set data to shared */ ALWAYS_INLINE void mpd_set_shared_data(mpd_t *result) { … } /* Set data to const */ ALWAYS_INLINE void mpd_set_const_data(mpd_t *result) { … } /* Clear flags, preserving memory attributes. */ ALWAYS_INLINE void mpd_clear_flags(mpd_t *result) { … } /* Set flags, preserving memory attributes. */ ALWAYS_INLINE void mpd_set_flags(mpd_t *result, uint8_t flags) { … } /* Copy flags, preserving memory attributes of result. */ ALWAYS_INLINE void mpd_copy_flags(mpd_t *result, const mpd_t *a) { … } /* Initialize a workcontext from ctx. Set traps, flags and newtrap to 0. */ static inline void mpd_workcontext(mpd_context_t *workctx, const mpd_context_t *ctx) { … } /******************************************************************************/ /* Getting and setting parts of decimals */ /******************************************************************************/ /* Flip the sign of a decimal */ static inline void _mpd_negate(mpd_t *dec) { … } /* Set coefficient to zero */ void mpd_zerocoeff(mpd_t *result) { … } /* Set the coefficient to all nines. */ void mpd_qmaxcoeff(mpd_t *result, const mpd_context_t *ctx, uint32_t *status) { … } /* * Cut off the most significant digits so that the rest fits in ctx->prec. * Cannot fail. */ static void _mpd_cap(mpd_t *result, const mpd_context_t *ctx) { … } /* * Cut off the most significant digits of a NaN payload so that the rest * fits in ctx->prec - ctx->clamp. Cannot fail. */ static void _mpd_fix_nan(mpd_t *result, const mpd_context_t *ctx) { … } /* * Get n most significant digits from a decimal, where 0 < n <= MPD_UINT_DIGITS. * Assumes MPD_UINT_DIGITS == MPD_RDIGITS+1, which is true for 32 and 64 bit * machines. * * The result of the operation will be in lo. If the operation is impossible, * hi will be nonzero. This is used to indicate an error. */ static inline void _mpd_get_msdigits(mpd_uint_t *hi, mpd_uint_t *lo, const mpd_t *dec, unsigned int n) { … } /******************************************************************************/ /* Gathering information about a decimal */ /******************************************************************************/ /* The real size of the coefficient without leading zero words. */ static inline mpd_ssize_t _mpd_real_size(mpd_uint_t *data, mpd_ssize_t size) { … } /* Return number of trailing zeros. No errors are possible. */ mpd_ssize_t mpd_trail_zeros(const mpd_t *dec) { … } /* Integer: Undefined for specials */ static int _mpd_isint(const mpd_t *dec) { … } /* Integer */ int mpd_isinteger(const mpd_t *dec) { … } /* Word is a power of 10 */ static int mpd_word_ispow10(mpd_uint_t word) { … } /* Coefficient is a power of 10 */ static int mpd_coeff_ispow10(const mpd_t *dec) { … } /* All digits of a word are nines */ static int mpd_word_isallnine(mpd_uint_t word) { … } /* All digits of the coefficient are nines */ static int mpd_coeff_isallnine(const mpd_t *dec) { … } /* Odd decimal: Undefined for non-integers! */ int mpd_isodd(const mpd_t *dec) { … } /* Even: Undefined for non-integers! */ int mpd_iseven(const mpd_t *dec) { … } /******************************************************************************/ /* Getting and setting decimals */ /******************************************************************************/ /* Internal function: Set a static decimal from a triple, no error checking. */ static void _ssettriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp) { … } /* Internal function: Set a decimal from a triple, no error checking. */ static void _settriple(mpd_t *result, uint8_t sign, mpd_uint_t a, mpd_ssize_t exp) { … } /* Set a special number from a triple */ void mpd_setspecial(mpd_t *result, uint8_t sign, uint8_t type) { … } /* Set result of NaN with an error status */ void mpd_seterror(mpd_t *result, uint32_t flags, uint32_t *status) { … } /* quietly set a static decimal from an mpd_ssize_t */ void mpd_qsset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a static decimal from an mpd_uint_t */ void mpd_qsset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a static decimal from an int32_t */ void mpd_qsset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a static decimal from a uint32_t */ void mpd_qsset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx, uint32_t *status) { … } #ifdef CONFIG_64 /* quietly set a static decimal from an int64_t */ void mpd_qsset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a static decimal from a uint64_t */ void mpd_qsset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx, uint32_t *status) { … } #endif /* quietly set a decimal from an mpd_ssize_t */ void mpd_qset_ssize(mpd_t *result, mpd_ssize_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a decimal from an mpd_uint_t */ void mpd_qset_uint(mpd_t *result, mpd_uint_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a decimal from an int32_t */ void mpd_qset_i32(mpd_t *result, int32_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a decimal from a uint32_t */ void mpd_qset_u32(mpd_t *result, uint32_t a, const mpd_context_t *ctx, uint32_t *status) { … } #if defined(CONFIG_32) && !defined(LEGACY_COMPILER) /* set a decimal from a uint64_t */ static void _c32setu64(mpd_t *result, uint64_t u, uint8_t sign, uint32_t *status) { mpd_uint_t w[3]; uint64_t q; int i, len; len = 0; do { q = u / MPD_RADIX; w[len] = (mpd_uint_t)(u - q * MPD_RADIX); u = q; len++; } while (u != 0); if (!mpd_qresize(result, len, status)) { return; } for (i = 0; i < len; i++) { result->data[i] = w[i]; } mpd_set_flags(result, sign); result->exp = 0; result->len = len; mpd_setdigits(result); } static void _c32_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx, uint32_t *status) { _c32setu64(result, a, MPD_POS, status); mpd_qfinalize(result, ctx, status); } /* set a decimal from an int64_t */ static void _c32_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx, uint32_t *status) { uint64_t u; uint8_t sign = MPD_POS; if (a < 0) { if (a == INT64_MIN) { u = (uint64_t)INT64_MAX + (-(INT64_MIN+INT64_MAX)); } else { u = -a; } sign = MPD_NEG; } else { u = a; } _c32setu64(result, u, sign, status); mpd_qfinalize(result, ctx, status); } #endif /* CONFIG_32 && !LEGACY_COMPILER */ #ifndef LEGACY_COMPILER /* quietly set a decimal from an int64_t */ void mpd_qset_i64(mpd_t *result, int64_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a decimal from an int64_t, use a maxcontext for conversion */ void mpd_qset_i64_exact(mpd_t *result, int64_t a, uint32_t *status) { … } /* quietly set a decimal from a uint64_t */ void mpd_qset_u64(mpd_t *result, uint64_t a, const mpd_context_t *ctx, uint32_t *status) { … } /* quietly set a decimal from a uint64_t, use a maxcontext for conversion */ void mpd_qset_u64_exact(mpd_t *result, uint64_t a, uint32_t *status) { … } #endif /* !LEGACY_COMPILER */ /* * Quietly get an mpd_uint_t from a decimal. Assumes * MPD_UINT_DIGITS == MPD_RDIGITS+1, which is true for * 32 and 64 bit machines. * * If the operation is impossible, MPD_Invalid_operation is set. */ static mpd_uint_t _mpd_qget_uint(int use_sign, const mpd_t *a, uint32_t *status) { … } /* * Sets Invalid_operation for: * - specials * - negative numbers (except negative zero) * - non-integers * - overflow */ mpd_uint_t mpd_qget_uint(const mpd_t *a, uint32_t *status) { … } /* Same as above, but gets the absolute value, i.e. the sign is ignored. */ mpd_uint_t mpd_qabs_uint(const mpd_t *a, uint32_t *status) { … } /* quietly get an mpd_ssize_t from a decimal */ mpd_ssize_t mpd_qget_ssize(const mpd_t *a, uint32_t *status) { … } #if defined(CONFIG_32) && !defined(LEGACY_COMPILER) /* * Quietly get a uint64_t from a decimal. If the operation is impossible, * MPD_Invalid_operation is set. */ static uint64_t _c32_qget_u64(int use_sign, const mpd_t *a, uint32_t *status) { MPD_NEW_STATIC(tmp,0,0,20,3); mpd_context_t maxcontext; uint64_t ret; tmp_data[0] = 709551615; tmp_data[1] = 446744073; tmp_data[2] = 18; if (mpd_isspecial(a)) { *status |= MPD_Invalid_operation; return UINT64_MAX; } if (mpd_iszero(a)) { return 0; } if (use_sign && mpd_isnegative(a)) { *status |= MPD_Invalid_operation; return UINT64_MAX; } if (!_mpd_isint(a)) { *status |= MPD_Invalid_operation; return UINT64_MAX; } if (_mpd_cmp_abs(a, &tmp) > 0) { *status |= MPD_Invalid_operation; return UINT64_MAX; } mpd_maxcontext(&maxcontext); mpd_qrescale(&tmp, a, 0, &maxcontext, &maxcontext.status); maxcontext.status &= ~MPD_Rounded; if (maxcontext.status != 0) { *status |= (maxcontext.status|MPD_Invalid_operation); /* GCOV_NOT_REACHED */ return UINT64_MAX; /* GCOV_NOT_REACHED */ } ret = 0; switch (tmp.len) { case 3: ret += (uint64_t)tmp_data[2] * 1000000000000000000ULL; case 2: ret += (uint64_t)tmp_data[1] * 1000000000ULL; case 1: ret += tmp_data[0]; break; default: abort(); /* GCOV_NOT_REACHED */ } return ret; } static int64_t _c32_qget_i64(const mpd_t *a, uint32_t *status) { uint64_t u; int isneg; u = _c32_qget_u64(0, a, status); if (*status&MPD_Invalid_operation) { return INT64_MAX; } isneg = mpd_isnegative(a); if (u <= INT64_MAX) { return isneg ? -((int64_t)u) : (int64_t)u; } else if (isneg && u+(INT64_MIN+INT64_MAX) == INT64_MAX) { return INT64_MIN; } *status |= MPD_Invalid_operation; return INT64_MAX; } #endif /* CONFIG_32 && !LEGACY_COMPILER */ #ifdef CONFIG_64 /* quietly get a uint64_t from a decimal */ uint64_t mpd_qget_u64(const mpd_t *a, uint32_t *status) { … } /* quietly get an int64_t from a decimal */ int64_t mpd_qget_i64(const mpd_t *a, uint32_t *status) { … } /* quietly get a uint32_t from a decimal */ uint32_t mpd_qget_u32(const mpd_t *a, uint32_t *status) { … } /* quietly get an int32_t from a decimal */ int32_t mpd_qget_i32(const mpd_t *a, uint32_t *status) { … } #else #ifndef LEGACY_COMPILER /* quietly get a uint64_t from a decimal */ uint64_t mpd_qget_u64(const mpd_t *a, uint32_t *status) { uint32_t workstatus = 0; uint64_t x = _c32_qget_u64(1, a, &workstatus); *status |= workstatus; return x; } /* quietly get an int64_t from a decimal */ int64_t mpd_qget_i64(const mpd_t *a, uint32_t *status) { uint32_t workstatus = 0; int64_t x = _c32_qget_i64(a, &workstatus); *status |= workstatus; return x; } #endif /* quietly get a uint32_t from a decimal */ uint32_t mpd_qget_u32(const mpd_t *a, uint32_t *status) { return mpd_qget_uint(a, status); } /* quietly get an int32_t from a decimal */ int32_t mpd_qget_i32(const mpd_t *a, uint32_t *status) { return mpd_qget_ssize(a, status); } #endif /******************************************************************************/ /* Filtering input of functions, finalizing output of functions */ /******************************************************************************/ /* * Check if the operand is NaN, copy to result and return 1 if this is * the case. Copying can fail since NaNs are allowed to have a payload that * does not fit in MPD_MINALLOC. */ int mpd_qcheck_nan(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * Check if either operand is NaN, copy to result and return 1 if this * is the case. Copying can fail since NaNs are allowed to have a payload * that does not fit in MPD_MINALLOC. */ int mpd_qcheck_nans(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Check if one of the operands is NaN, copy to result and return 1 if this * is the case. Copying can fail since NaNs are allowed to have a payload * that does not fit in MPD_MINALLOC. */ static int mpd_qcheck_3nans(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c, const mpd_context_t *ctx, uint32_t *status) { … } /* Check if rounding digit 'rnd' leads to an increment. */ static inline int _mpd_rnd_incr(const mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx) { … } /* * Apply rounding to a decimal that has been right-shifted into a full * precision decimal. If an increment leads to an overflow of the precision, * adjust the coefficient and the exponent and check the new exponent for * overflow. */ static inline void _mpd_apply_round(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx, uint32_t *status) { … } /* * Apply rounding to a decimal. Allow overflow of the precision. */ static inline void _mpd_apply_round_excess(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx, uint32_t *status) { … } /* * Apply rounding to a decimal that has been right-shifted into a decimal * with full precision or less. Return failure if an increment would * overflow the precision. */ static inline int _mpd_apply_round_fit(mpd_t *dec, mpd_uint_t rnd, const mpd_context_t *ctx, uint32_t *status) { … } /* Check a normal number for overflow, underflow, clamping. If the operand is modified, it will be zero, special or (sub)normal with a coefficient that fits into the current context precision. */ static inline void _mpd_check_exp(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status) { … } /* Transcendental functions do not always set Underflow reliably, * since they only use as much precision as is necessary for correct * rounding. If a result like 1.0000000000e-101 is finalized, there * is no rounding digit that would trigger Underflow. But we can * assume Inexact, so a short check suffices. */ static inline void mpd_check_underflow(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status) { … } /* Check if a normal number must be rounded after the exponent has been checked. */ static inline void _mpd_check_round(mpd_t *dec, const mpd_context_t *ctx, uint32_t *status) { … } /* Finalize all operations. */ void mpd_qfinalize(mpd_t *result, const mpd_context_t *ctx, uint32_t *status) { … } /******************************************************************************/ /* Copying */ /******************************************************************************/ /* Internal function: Copy a decimal, share data with src: USE WITH CARE! */ static inline void _mpd_copy_shared(mpd_t *dest, const mpd_t *src) { … } /* * Copy a decimal. In case of an error, status is set to MPD_Malloc_error. */ int mpd_qcopy(mpd_t *result, const mpd_t *a, uint32_t *status) { … } /* Same as mpd_qcopy, but do not set the result to NaN on failure. */ int mpd_qcopy_cxx(mpd_t *result, const mpd_t *a) { … } /* * Copy to a decimal with a static buffer. The caller has to make sure that * the buffer is big enough. Cannot fail. */ static void mpd_qcopy_static(mpd_t *result, const mpd_t *a) { … } /* * Return a newly allocated copy of the operand. In case of an error, * status is set to MPD_Malloc_error and the return value is NULL. */ mpd_t * mpd_qncopy(const mpd_t *a) { … } /* * Copy a decimal and set the sign to positive. In case of an error, the * status is set to MPD_Malloc_error. */ int mpd_qcopy_abs(mpd_t *result, const mpd_t *a, uint32_t *status) { … } /* * Copy a decimal and negate the sign. In case of an error, the * status is set to MPD_Malloc_error. */ int mpd_qcopy_negate(mpd_t *result, const mpd_t *a, uint32_t *status) { … } /* * Copy a decimal, setting the sign of the first operand to the sign of the * second operand. In case of an error, the status is set to MPD_Malloc_error. */ int mpd_qcopy_sign(mpd_t *result, const mpd_t *a, const mpd_t *b, uint32_t *status) { … } /******************************************************************************/ /* Comparisons */ /******************************************************************************/ /* * For all functions that compare two operands and return an int the usual * convention applies to the return value: * * -1 if op1 < op2 * 0 if op1 == op2 * 1 if op1 > op2 * * INT_MAX for error */ /* Convenience macro. If a and b are not equal, return from the calling * function with the correct comparison value. */ #define CMP_EQUAL_OR_RETURN(a, b) … /* * Compare the data of big and small. This function does the equivalent * of first shifting small to the left and then comparing the data of * big and small, except that no allocation for the left shift is needed. */ static int _mpd_basecmp(mpd_uint_t *big, mpd_uint_t *small, mpd_size_t n, mpd_size_t m, mpd_size_t shift) { … } /* Compare two decimals with the same adjusted exponent. */ static int _mpd_cmp_same_adjexp(const mpd_t *a, const mpd_t *b) { … } /* Compare two numerical values. */ static int _mpd_cmp(const mpd_t *a, const mpd_t *b) { … } /* Compare the absolutes of two numerical values. */ static int _mpd_cmp_abs(const mpd_t *a, const mpd_t *b) { … } /* Compare two values and return an integer result. */ int mpd_qcmp(const mpd_t *a, const mpd_t *b, uint32_t *status) { … } /* * Compare a and b, convert the usual integer result to a decimal and * store it in 'result'. For convenience, the integer result of the comparison * is returned. Comparisons involving NaNs return NaN/INT_MAX. */ int mpd_qcompare(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Same as mpd_compare(), but signal for all NaNs, i.e. also for quiet NaNs. */ int mpd_qcompare_signal(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Compare the operands using a total order. */ int mpd_cmp_total(const mpd_t *a, const mpd_t *b) { … } /* * Compare a and b according to a total order, convert the usual integer result * to a decimal and store it in 'result'. For convenience, the integer result * of the comparison is returned. */ int mpd_compare_total(mpd_t *result, const mpd_t *a, const mpd_t *b) { … } /* Compare the magnitude of the operands using a total order. */ int mpd_cmp_total_mag(const mpd_t *a, const mpd_t *b) { … } /* * Compare the magnitude of a and b according to a total order, convert the * the usual integer result to a decimal and store it in 'result'. * For convenience, the integer result of the comparison is returned. */ int mpd_compare_total_mag(mpd_t *result, const mpd_t *a, const mpd_t *b) { … } /* Determine an ordering for operands that are numerically equal. */ static inline int _mpd_cmp_numequal(const mpd_t *a, const mpd_t *b) { … } /******************************************************************************/ /* Shifting the coefficient */ /******************************************************************************/ /* * Shift the coefficient of the operand to the left, no check for specials. * Both operands may be the same pointer. If the result length has to be * increased, mpd_qresize() might fail with MPD_Malloc_error. */ int mpd_qshiftl(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status) { … } /* Determine the rounding indicator if all digits of the coefficient are shifted * out of the picture. */ static mpd_uint_t _mpd_get_rnd(const mpd_uint_t *data, mpd_ssize_t len, int use_msd) { … } /* * Same as mpd_qshiftr(), but 'result' is an mpd_t with a static coefficient. * It is the caller's responsibility to ensure that the coefficient is big * enough. The function cannot fail. */ static mpd_uint_t mpd_qsshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n) { … } /* * Inplace shift of the coefficient to the right, no check for specials. * Returns the rounding indicator for mpd_rnd_incr(). * The function cannot fail. */ mpd_uint_t mpd_qshiftr_inplace(mpd_t *result, mpd_ssize_t n) { … } /* * Shift the coefficient of the operand to the right, no check for specials. * Both operands may be the same pointer. Returns the rounding indicator to * be used by mpd_rnd_incr(). If the result length has to be increased, * mpd_qcopy() or mpd_qresize() might fail with MPD_Malloc_error. In those * cases, MPD_UINT_MAX is returned. */ mpd_uint_t mpd_qshiftr(mpd_t *result, const mpd_t *a, mpd_ssize_t n, uint32_t *status) { … } /******************************************************************************/ /* Miscellaneous operations */ /******************************************************************************/ /* Logical And */ void mpd_qand(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Class of an operand. Returns a pointer to the constant name. */ const char * mpd_class(const mpd_t *a, const mpd_context_t *ctx) { … } /* Logical Xor */ void mpd_qinvert(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* Exponent of the magnitude of the most significant digit of the operand. */ void mpd_qlogb(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* Logical Or */ void mpd_qor(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Rotate the coefficient of 'a' by 'b' digits. 'b' must be an integer with * exponent 0. */ void mpd_qrotate(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * b must be an integer with exponent 0 and in the range +-2*(emax + prec). * XXX: In my opinion +-(2*emax + prec) would be more sensible. * The result is a with the value of b added to its exponent. */ void mpd_qscaleb(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Shift the coefficient by n digits, positive n is a left shift. In the case * of a left shift, the result is decapitated to fit the context precision. If * you don't want that, use mpd_shiftl(). */ void mpd_qshiftn(mpd_t *result, const mpd_t *a, mpd_ssize_t n, const mpd_context_t *ctx, uint32_t *status) { … } /* * Same as mpd_shiftn(), but the shift is specified by the decimal b, which * must be an integer with a zero exponent. Infinities remain infinities. */ void mpd_qshift(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Logical Xor */ void mpd_qxor(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /******************************************************************************/ /* Arithmetic operations */ /******************************************************************************/ /* * The absolute value of a. If a is negative, the result is the same * as the result of the minus operation. Otherwise, the result is the * result of the plus operation. */ void mpd_qabs(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } static inline void _mpd_ptrswap(const mpd_t **a, const mpd_t **b) { … } /* Add or subtract infinities. */ static void _mpd_qaddsub_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b, uint32_t *status) { … } /* Add or subtract non-special numbers. */ static void _mpd_qaddsub(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add a and b. No specials, no finalizing. */ static void _mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract b from a. No specials, no finalizing. */ static void _mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add a and b. */ void mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add a and b. Set NaN/Invalid_operation if the result is inexact. */ static void _mpd_qadd_exact(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract b from a. */ void mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract b from a. Set NaN/Invalid_operation if the result is inexact. */ static void _mpd_qsub_exact(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add decimal and mpd_ssize_t. */ void mpd_qadd_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add decimal and mpd_uint_t. */ void mpd_qadd_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract mpd_ssize_t from decimal. */ void mpd_qsub_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract mpd_uint_t from decimal. */ void mpd_qsub_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add decimal and int32_t. */ void mpd_qadd_i32(mpd_t *result, const mpd_t *a, int32_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add decimal and uint32_t. */ void mpd_qadd_u32(mpd_t *result, const mpd_t *a, uint32_t b, const mpd_context_t *ctx, uint32_t *status) { … } #ifdef CONFIG_64 /* Add decimal and int64_t. */ void mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Add decimal and uint64_t. */ void mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { … } #elif !defined(LEGACY_COMPILER) /* Add decimal and int64_t. */ void mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_i64(&bb, b, &maxcontext, status); mpd_qadd(result, a, &bb, ctx, status); mpd_del(&bb); } /* Add decimal and uint64_t. */ void mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_u64(&bb, b, &maxcontext, status); mpd_qadd(result, a, &bb, ctx, status); mpd_del(&bb); } #endif /* Subtract int32_t from decimal. */ void mpd_qsub_i32(mpd_t *result, const mpd_t *a, int32_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract uint32_t from decimal. */ void mpd_qsub_u32(mpd_t *result, const mpd_t *a, uint32_t b, const mpd_context_t *ctx, uint32_t *status) { … } #ifdef CONFIG_64 /* Subtract int64_t from decimal. */ void mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Subtract uint64_t from decimal. */ void mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { … } #elif !defined(LEGACY_COMPILER) /* Subtract int64_t from decimal. */ void mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_i64(&bb, b, &maxcontext, status); mpd_qsub(result, a, &bb, ctx, status); mpd_del(&bb); } /* Subtract uint64_t from decimal. */ void mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_u64(&bb, b, &maxcontext, status); mpd_qsub(result, a, &bb, ctx, status); mpd_del(&bb); } #endif /* Divide infinities. */ static void _mpd_qdiv_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } enum { … }; /* Divide a by b. */ static void _mpd_qdiv(int action, mpd_t *q, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide a by b. */ void mpd_qdiv(mpd_t *q, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Internal function. */ static void _mpd_qdivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Integer division with remainder. */ void mpd_qdivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qdivint(mpd_t *q, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide decimal by mpd_ssize_t. */ void mpd_qdiv_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide decimal by mpd_uint_t. */ void mpd_qdiv_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide decimal by int32_t. */ void mpd_qdiv_i32(mpd_t *result, const mpd_t *a, int32_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide decimal by uint32_t. */ void mpd_qdiv_u32(mpd_t *result, const mpd_t *a, uint32_t b, const mpd_context_t *ctx, uint32_t *status) { … } #ifdef CONFIG_64 /* Divide decimal by int64_t. */ void mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Divide decimal by uint64_t. */ void mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { … } #elif !defined(LEGACY_COMPILER) /* Divide decimal by int64_t. */ void mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_i64(&bb, b, &maxcontext, status); mpd_qdiv(result, a, &bb, ctx, status); mpd_del(&bb); } /* Divide decimal by uint64_t. */ void mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_u64(&bb, b, &maxcontext, status); mpd_qdiv(result, a, &bb, ctx, status); mpd_del(&bb); } #endif /* Pad the result with trailing zeros if it has fewer digits than prec. */ static void _mpd_zeropad(mpd_t *result, const mpd_context_t *ctx, uint32_t *status) { … } /* Check if the result is guaranteed to be one. */ static int _mpd_qexp_check_one(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * Get the number of iterations for the Horner scheme in _mpd_qexp(). */ static inline mpd_ssize_t _mpd_get_exp_iterations(const mpd_t *r, mpd_ssize_t p) { … } /* * Internal function, specials have been dealt with. Apart from Overflow * and Underflow, two cases must be considered for the error of the result: * * 1) abs(a) <= 9 * 10**(-prec-1) ==> result == 1 * * Absolute error: abs(1 - e**x) < 10**(-prec) * ------------------------------------------- * * 2) abs(a) > 9 * 10**(-prec-1) * * Relative error: abs(result - e**x) < 0.5 * 10**(-prec) * e**x * ------------------------------------------------------------- * * The algorithm is from Hull&Abrham, Variable Precision Exponential Function, * ACM Transactions on Mathematical Software, Vol. 12, No. 2, June 1986. * * Main differences: * * - The number of iterations for the Horner scheme is calculated using * 53-bit floating point arithmetic. * * - In the error analysis for ER (relative error accumulated in the * evaluation of the truncated series) the reduced operand r may * have any number of digits. * ACL2 proof: exponent-relative-error * * - The analysis for early abortion has been adapted for the mpd_t * ranges. */ static void _mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* exp(a) */ void mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* Fused multiply-add: (a * b) + c, with a single final rounding. */ void mpd_qfma(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_t *c, const mpd_context_t *ctx, uint32_t *status) { … } /* * Schedule the optimal precision increase for the Newton iteration. * v := input operand * z_0 := initial approximation * initprec := natural number such that abs(log(v) - z_0) < 10**-initprec * maxprec := target precision * * For convenience the output klist contains the elements in reverse order: * klist := [k_n-1, ..., k_0], where * 1) k_0 <= initprec and * 2) abs(log(v) - result) < 10**(-2*k_n-1 + 1) <= 10**-maxprec. */ static inline int ln_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2], mpd_ssize_t maxprec, mpd_ssize_t initprec) { … } /* The constants have been verified with both decimal.py and mpfr. */ #ifdef CONFIG_64 #if MPD_RDIGITS != 19 #error "mpdecimal.c: MPD_RDIGITS must be 19." #endif static const mpd_uint_t mpd_ln10_data[MPD_MINALLOC_MAX] = …; #else #if MPD_RDIGITS != 9 #error "mpdecimal.c: MPD_RDIGITS must be 9." #endif static const mpd_uint_t mpd_ln10_data[MPD_MINALLOC_MAX] = { 401682692UL, 708474699UL, 720754403UL, 30896345UL, 602301057UL, 765871416UL, 192920333UL, 763113569UL, 589402567UL, 956890167UL, 82413146UL, 589257242UL, 245544057UL, 811364292UL, 734206705UL, 868569356UL, 167465505UL, 775026849UL, 706480002UL, 18064450UL, 636167921UL, 569476834UL, 734507478UL, 156591213UL, 148046637UL, 283552201UL, 677432162UL, 470806855UL, 880840126UL, 417480036UL, 210510171UL, 940440022UL, 939147961UL, 893431493UL, 436515504UL, 440424327UL, 654366747UL, 821988674UL, 622228769UL, 884616336UL, 537773262UL, 350530896UL, 319852839UL, 989482623UL, 468084379UL, 720832555UL, 168948290UL, 736909878UL, 675666628UL, 546508280UL, 863340952UL, 404228624UL, 834196778UL, 508959829UL, 23599720UL, 967735248UL, 96757260UL, 603332790UL, 862877297UL, 760110148UL, 468436420UL, 401799145UL, 299404568UL, 230258509UL }; #endif /* _mpd_ln10 is used directly for precisions smaller than MINALLOC_MAX*RDIGITS. Otherwise, it serves as the initial approximation for calculating ln(10). */ static const mpd_t _mpd_ln10 = …; /* * Set 'result' to log(10). * Ulp error: abs(result - log(10)) < ulp(log(10)) * Relative error: abs(result - log(10)) < 5 * 10**-prec * log(10) * * NOTE: The relative error is not derived from the ulp error, but * calculated separately using the fact that 23/10 < log(10) < 24/10. */ void mpd_qln10(mpd_t *result, mpd_ssize_t prec, uint32_t *status) { … } /* * Initial approximations for the ln() iteration. The values have the * following properties (established with both decimal.py and mpfr): * * Index 0 - 400, logarithms of x in [1.00, 5.00]: * abs(lnapprox[i] * 10**-3 - log((i+100)/100)) < 10**-2 * abs(lnapprox[i] * 10**-3 - log((i+1+100)/100)) < 10**-2 * * Index 401 - 899, logarithms of x in (0.500, 0.999]: * abs(-lnapprox[i] * 10**-3 - log((i+100)/1000)) < 10**-2 * abs(-lnapprox[i] * 10**-3 - log((i+1+100)/1000)) < 10**-2 */ static const uint16_t lnapprox[900] = …; /* * Internal ln() function that does not check for specials, zero or one. * Relative error: abs(result - log(a)) < 0.1 * 10**-prec * abs(log(a)) */ static void _mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* ln(a) */ void mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * Internal log10() function that does not check for specials, zero or one. * Case SKIP_FINALIZE: * Relative error: abs(result - log10(a)) < 0.1 * 10**-prec * abs(log10(a)) * Case DO_FINALIZE: * Ulp error: abs(result - log10(a)) < ulp(log10(a)) */ enum { … }; static void _mpd_qlog10(int action, mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* log10(a) */ void mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * Maximum of the two operands. Attention: If one operand is a quiet NaN and the * other is numeric, the numeric operand is returned. This may not be what one * expects. */ void mpd_qmax(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Maximum magnitude: Same as mpd_max(), but compares the operands with their * sign ignored. */ void mpd_qmax_mag(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Minimum of the two operands. Attention: If one operand is a quiet NaN and the * other is numeric, the numeric operand is returned. This may not be what one * expects. */ void mpd_qmin(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Minimum magnitude: Same as mpd_min(), but compares the operands with their * sign ignored. */ void mpd_qmin_mag(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Minimum space needed for the result array in _karatsuba_rec(). */ static inline mpd_size_t _kmul_resultsize(mpd_size_t la, mpd_size_t lb) { … } /* Work space needed in _karatsuba_rec(). lim >= 4 */ static inline mpd_size_t _kmul_worksize(mpd_size_t n, mpd_size_t lim) { … } #define MPD_KARATSUBA_BASECASE … /* * Add the product of a and b to c. * c must be _kmul_resultsize(la, lb) in size. * w is used as a work array and must be _kmul_worksize(a, lim) in size. * Roman E. Maeder, Storage Allocation for the Karatsuba Integer Multiplication * Algorithm. In "Design and implementation of symbolic computation systems", * Springer, 1993, ISBN 354057235X, 9783540572350. */ static void _karatsuba_rec(mpd_uint_t *c, const mpd_uint_t *a, const mpd_uint_t *b, mpd_uint_t *w, mpd_size_t la, mpd_size_t lb) { … } /* * Multiply u and v, using Karatsuba multiplication. Returns a pointer * to the result or NULL in case of failure (malloc error). * Conditions: ulen >= vlen, ulen >= 4 */ static mpd_uint_t * _mpd_kmul(const mpd_uint_t *u, const mpd_uint_t *v, mpd_size_t ulen, mpd_size_t vlen, mpd_size_t *rsize) { … } /* * Determine the minimum length for the number theoretic transform. Valid * transform lengths are 2**n or 3*2**n, where 2**n <= MPD_MAXTRANSFORM_2N. * The function finds the shortest length m such that rsize <= m. */ static inline mpd_size_t _mpd_get_transform_len(mpd_size_t rsize) { … } #ifdef PPRO #ifndef _MSC_VER static inline unsigned short _mpd_get_control87(void) { unsigned short cw; __asm__ __volatile__ ("fnstcw %0" : "=m" (cw)); return cw; } static inline void _mpd_set_control87(unsigned short cw) { __asm__ __volatile__ ("fldcw %0" : : "m" (cw)); } #endif static unsigned int mpd_set_fenv(void) { unsigned int cw; #ifdef _MSC_VER unsigned int flags = _EM_INVALID|_EM_DENORMAL|_EM_ZERODIVIDE|_EM_OVERFLOW| _EM_UNDERFLOW|_EM_INEXACT|_RC_CHOP|_PC_64; unsigned int mask = _MCW_EM|_MCW_RC|_MCW_PC; unsigned int dummy; __control87_2(0, 0, &cw, NULL); __control87_2(flags, mask, &dummy, NULL); #else cw = _mpd_get_control87(); _mpd_set_control87(cw|0xF3F); #endif return cw; } static void mpd_restore_fenv(unsigned int cw) { #ifdef _MSC_VER unsigned int mask = _MCW_EM|_MCW_RC|_MCW_PC; unsigned int dummy; __control87_2(cw, mask, &dummy, NULL); #else _mpd_set_control87((unsigned short)cw); #endif } #endif /* PPRO */ /* * Multiply u and v, using the fast number theoretic transform. Returns * a pointer to the result or NULL in case of failure (malloc error). */ static mpd_uint_t * _mpd_fntmul(const mpd_uint_t *u, const mpd_uint_t *v, mpd_size_t ulen, mpd_size_t vlen, mpd_size_t *rsize) { … } /* * Karatsuba multiplication with FNT/basemul as the base case. */ static int _karatsuba_rec_fnt(mpd_uint_t *c, const mpd_uint_t *a, const mpd_uint_t *b, mpd_uint_t *w, mpd_size_t la, mpd_size_t lb) { … } /* * Multiply u and v, using Karatsuba multiplication with the FNT as the * base case. Returns a pointer to the result or NULL in case of failure * (malloc error). Conditions: ulen >= vlen, ulen >= 4. */ static mpd_uint_t * _mpd_kmul_fnt(const mpd_uint_t *u, const mpd_uint_t *v, mpd_size_t ulen, mpd_size_t vlen, mpd_size_t *rsize) { … } /* Deal with the special cases of multiplying infinities. */ static void _mpd_qmul_inf(mpd_t *result, const mpd_t *a, const mpd_t *b, uint32_t *status) { … } /* * Internal function: Multiply a and b. _mpd_qmul deals with specials but * does NOT finalize the result. This is for use in mpd_fma(). */ static inline void _mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Multiply a and b. */ void mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Multiply a and b. Set NaN/Invalid_operation if the result is inexact. */ static void _mpd_qmul_exact(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* Multiply decimal and mpd_ssize_t. */ void mpd_qmul_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b, const mpd_context_t *ctx, uint32_t *status) { … } /* Multiply decimal and mpd_uint_t. */ void mpd_qmul_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qmul_i32(mpd_t *result, const mpd_t *a, int32_t b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qmul_u32(mpd_t *result, const mpd_t *a, uint32_t b, const mpd_context_t *ctx, uint32_t *status) { … } #ifdef CONFIG_64 void mpd_qmul_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qmul_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { … } #elif !defined(LEGACY_COMPILER) /* Multiply decimal and int64_t. */ void mpd_qmul_i64(mpd_t *result, const mpd_t *a, int64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_i64(&bb, b, &maxcontext, status); mpd_qmul(result, a, &bb, ctx, status); mpd_del(&bb); } /* Multiply decimal and uint64_t. */ void mpd_qmul_u64(mpd_t *result, const mpd_t *a, uint64_t b, const mpd_context_t *ctx, uint32_t *status) { mpd_context_t maxcontext; MPD_NEW_STATIC(bb,0,0,0,0); mpd_maxcontext(&maxcontext); mpd_qset_u64(&bb, b, &maxcontext, status); mpd_qmul(result, a, &bb, ctx, status); mpd_del(&bb); } #endif /* Like the minus operator. */ void mpd_qminus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* Like the plus operator. */ void mpd_qplus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* The largest representable number that is smaller than the operand. */ void mpd_qnext_minus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* The smallest representable number that is larger than the operand. */ void mpd_qnext_plus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * The number closest to the first operand that is in the direction towards * the second operand. */ void mpd_qnext_toward(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } /* * Internal function: Integer power with mpd_uint_t exponent. The function * can fail with MPD_Malloc_error. * * The error is equal to the error incurred in k-1 multiplications. Assuming * the upper bound for the relative error in each operation: * * abs(err) = 5 * 10**-prec * result = x**k * (1 + err)**(k-1) */ static inline void _mpd_qpow_uint(mpd_t *result, const mpd_t *base, mpd_uint_t exp, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * Internal function: Integer power with mpd_t exponent, tbase and texp * are modified!! Function can fail with MPD_Malloc_error. * * The error is equal to the error incurred in k multiplications. Assuming * the upper bound for the relative error in each operation: * * abs(err) = 5 * 10**-prec * result = x**k * (1 + err)**k */ static inline void _mpd_qpow_mpd(mpd_t *result, mpd_t *tbase, mpd_t *texp, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * The power function for integer exponents. Relative error _before_ the * final rounding to prec: * abs(result - base**exp) < 0.1 * 10**-prec * abs(base**exp) */ static void _mpd_qpow_int(mpd_t *result, const mpd_t *base, const mpd_t *exp, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * If the exponent is infinite and base equals one, the result is one * with a coefficient of length prec. Otherwise, result is undefined. * Return the value of the comparison against one. */ static int _qcheck_pow_one_inf(mpd_t *result, const mpd_t *base, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * If abs(base) equals one, calculate the correct power of one result. * Otherwise, result is undefined. Return the value of the comparison * against 1. * * This is an internal function that does not check for specials. */ static int _qcheck_pow_one(mpd_t *result, const mpd_t *base, const mpd_t *exp, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * Detect certain over/underflow of x**y. * ACL2 proof: pow-bounds.lisp. * * Symbols: * * e: EXP_INF or EXP_CLAMP * x: base * y: exponent * * omega(e) = log10(abs(e)) * zeta(x) = log10(abs(log10(x))) * theta(y) = log10(abs(y)) * * Upper and lower bounds: * * ub_omega(e) = ceil(log10(abs(e))) * lb_theta(y) = floor(log10(abs(y))) * * | floor(log10(floor(abs(log10(x))))) if x < 1/10 or x >= 10 * lb_zeta(x) = | floor(log10(abs(x-1)/10)) if 1/10 <= x < 1 * | floor(log10(abs((x-1)/100))) if 1 < x < 10 * * ub_omega(e) and lb_theta(y) are obviously upper and lower bounds * for omega(e) and theta(y). * * lb_zeta is a lower bound for zeta(x): * * x < 1/10 or x >= 10: * * abs(log10(x)) >= 1, so the outer log10 is well defined. Since log10 * is strictly increasing, the end result is a lower bound. * * 1/10 <= x < 1: * * We use: log10(x) <= (x-1)/log(10) * abs(log10(x)) >= abs(x-1)/log(10) * abs(log10(x)) >= abs(x-1)/10 * * 1 < x < 10: * * We use: (x-1)/(x*log(10)) < log10(x) * abs((x-1)/100) < abs(log10(x)) * * XXX: abs((x-1)/10) would work, need ACL2 proof. * * * Let (0 < x < 1 and y < 0) or (x > 1 and y > 0). (H1) * Let ub_omega(exp_inf) < lb_zeta(x) + lb_theta(y) (H2) * * Then: * log10(abs(exp_inf)) < log10(abs(log10(x))) + log10(abs(y)). (1) * exp_inf < log10(x) * y (2) * 10**exp_inf < x**y (3) * * Let (0 < x < 1 and y > 0) or (x > 1 and y < 0). (H3) * Let ub_omega(exp_clamp) < lb_zeta(x) + lb_theta(y) (H4) * * Then: * log10(abs(exp_clamp)) < log10(abs(log10(x))) + log10(abs(y)). (4) * log10(x) * y < exp_clamp (5) * x**y < 10**exp_clamp (6) * */ static mpd_ssize_t _lower_bound_zeta(const mpd_t *x, uint32_t *status) { … } /* * Detect cases of certain overflow/underflow in the power function. * Assumptions: x != 1, y != 0. The proof above is for positive x. * If x is negative and y is an odd integer, x**y == -(abs(x)**y), * so the analysis does not change. */ static int _qcheck_pow_bounds(mpd_t *result, const mpd_t *x, const mpd_t *y, uint8_t resultsign, const mpd_context_t *ctx, uint32_t *status) { … } /* * TODO: Implement algorithm for computing exact powers from decimal.py. * In order to prevent infinite loops, this has to be called before * using Ziv's strategy for correct rounding. */ /* static int _mpd_qpow_exact(mpd_t *result, const mpd_t *base, const mpd_t *exp, const mpd_context_t *ctx, uint32_t *status) { return 0; } */ /* * The power function for real exponents. * Relative error: abs(result - e**y) < e**y * 1/5 * 10**(-prec - 1) */ static void _mpd_qpow_real(mpd_t *result, const mpd_t *base, const mpd_t *exp, const mpd_context_t *ctx, uint32_t *status) { … } /* The power function: base**exp */ void mpd_qpow(mpd_t *result, const mpd_t *base, const mpd_t *exp, const mpd_context_t *ctx, uint32_t *status) { … } /* * Internal function: Integer powmod with mpd_uint_t exponent, base is modified! * Function can fail with MPD_Malloc_error. */ static inline void _mpd_qpowmod_uint(mpd_t *result, mpd_t *base, mpd_uint_t exp, const mpd_t *mod, uint32_t *status) { … } /* The powmod function: (base**exp) % mod */ void mpd_qpowmod(mpd_t *result, const mpd_t *base, const mpd_t *exp, const mpd_t *mod, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qquantize(mpd_t *result, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qreduce(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qrem(mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qrem_near(mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_context_t *ctx, uint32_t *status) { … } static void _mpd_qrescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp, const mpd_context_t *ctx, uint32_t *status) { … } /* * Rescale a number so that it has exponent 'exp'. Does not regard context * precision, emax, emin, but uses the rounding mode. Special numbers are * quietly copied. Restrictions: * * MPD_MIN_ETINY <= exp <= MPD_MAX_EMAX+1 * result->digits <= MPD_MAX_PREC+1 */ void mpd_qrescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp, const mpd_context_t *ctx, uint32_t *status) { … } /* * Same as mpd_qrescale, but with relaxed restrictions. The result of this * function should only be used for formatting a number and never as input * for other operations. * * MPD_MIN_ETINY-MPD_MAX_PREC <= exp <= MPD_MAX_EMAX+1 * result->digits <= MPD_MAX_PREC+1 */ void mpd_qrescale_fmt(mpd_t *result, const mpd_t *a, mpd_ssize_t exp, const mpd_context_t *ctx, uint32_t *status) { … } /* Round to an integer according to 'action' and ctx->round. */ enum { … }; static void _mpd_qround_to_integral(int action, mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qround_to_intx(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qround_to_int(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qtrunc(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qfloor(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qceil(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } int mpd_same_quantum(const mpd_t *a, const mpd_t *b) { … } /* Schedule the increase in precision for the Newton iteration. */ static inline int recpr_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2], mpd_ssize_t maxprec, mpd_ssize_t initprec) { … } /* * Initial approximation for the reciprocal: * k_0 := MPD_RDIGITS-2 * z_0 := 10**(-k_0) * floor(10**(2*k_0 + 2) / floor(v * 10**(k_0 + 2))) * Absolute error: * |1/v - z_0| < 10**(-k_0) * ACL2 proof: maxerror-inverse-approx */ static void _mpd_qreciprocal_approx(mpd_t *z, const mpd_t *v, uint32_t *status) { … } /* * Reciprocal, calculated with Newton's Method. Assumption: result != a. * NOTE: The comments in the function show that certain operations are * exact. The proof for the maximum error is too long to fit in here. * ACL2 proof: maxerror-inverse-complete */ static void _mpd_qreciprocal(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* * Internal function for large numbers: * * q, r = divmod(coeff(a), coeff(b)) * * Strategy: Multiply the dividend by the reciprocal of the divisor. The * inexact result is fixed by a small loop, using at most one iteration. * * ACL2 proofs: * ------------ * 1) q is a natural number. (ndivmod-quotient-natp) * 2) r is a natural number. (ndivmod-remainder-natp) * 3) a = q * b + r (ndivmod-q*b+r==a) * 4) r < b (ndivmod-remainder-<-b) */ static void _mpd_base_ndivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b, uint32_t *status) { … } /* LIBMPDEC_ONLY */ /* * Schedule the optimal precision increase for the Newton iteration. * v := input operand * z_0 := initial approximation * initprec := natural number such that abs(sqrt(v) - z_0) < 10**-initprec * maxprec := target precision * * For convenience the output klist contains the elements in reverse order: * klist := [k_n-1, ..., k_0], where * 1) k_0 <= initprec and * 2) abs(sqrt(v) - result) < 10**(-2*k_n-1 + 2) <= 10**-maxprec. */ static inline int invroot_schedule_prec(mpd_ssize_t klist[MPD_MAX_PREC_LOG2], mpd_ssize_t maxprec, mpd_ssize_t initprec) { … } /* * Initial approximation for the inverse square root function. * Input: * v := rational number, with 1 <= v < 100 * vhat := floor(v * 10**6) * Output: * z := approximation to 1/sqrt(v), such that abs(z - 1/sqrt(v)) < 10**-3. */ static inline void _invroot_init_approx(mpd_t *z, mpd_uint_t vhat) { … } /* * Set 'result' to 1/sqrt(a). * Relative error: abs(result - 1/sqrt(a)) < 10**-prec * 1/sqrt(a) */ static void _mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /* END LIBMPDEC_ONLY */ /* Algorithm from decimal.py */ static void _mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } void mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, uint32_t *status) { … } /******************************************************************************/ /* Base conversions */ /******************************************************************************/ /* Space needed to represent an integer mpd_t in base 'base'. */ size_t mpd_sizeinbase(const mpd_t *a, uint32_t base) { … } /* Space needed to import a base 'base' integer of length 'srclen'. */ static mpd_ssize_t _mpd_importsize(size_t srclen, uint32_t base) { … } static uint8_t mpd_resize_u16(uint16_t **w, size_t nmemb) { … } static uint8_t mpd_resize_u32(uint32_t **w, size_t nmemb) { … } static size_t _baseconv_to_u16(uint16_t **w, size_t wlen, mpd_uint_t wbase, mpd_uint_t *u, mpd_ssize_t ulen) { … } static size_t _coeff_from_u16(mpd_t *w, mpd_ssize_t wlen, const mpd_uint_t *u, size_t ulen, uint32_t ubase, uint32_t *status) { … } /* target base wbase < source base ubase */ static size_t _baseconv_to_smaller(uint32_t **w, size_t wlen, uint32_t wbase, mpd_uint_t *u, mpd_ssize_t ulen, mpd_uint_t ubase) { … } #ifdef CONFIG_32 /* target base 'wbase' == source base 'ubase' */ static size_t _copy_equal_base(uint32_t **w, size_t wlen, const uint32_t *u, size_t ulen) { if (wlen < ulen) { if (!mpd_resize_u32(w, ulen)) { return SIZE_MAX; } } memcpy(*w, u, ulen * (sizeof **w)); return ulen; } /* target base 'wbase' > source base 'ubase' */ static size_t _baseconv_to_larger(uint32_t **w, size_t wlen, mpd_uint_t wbase, const mpd_uint_t *u, size_t ulen, mpd_uint_t ubase) { size_t n = 0; mpd_uint_t carry; assert(wlen > 0 && ulen > 0); assert(ubase < wbase); (*w)[n++] = u[--ulen]; while (--ulen != SIZE_MAX) { carry = _mpd_shortmul_b(*w, *w, n, ubase, wbase); if (carry) { if (n >= wlen) { if (!mpd_resize_u32(w, n+1)) { return SIZE_MAX; } wlen = n+1; } (*w)[n++] = carry; } carry = _mpd_shortadd_b(*w, n, u[ulen], wbase); if (carry) { if (n >= wlen) { if (!mpd_resize_u32(w, n+1)) { return SIZE_MAX; } wlen = n+1; } (*w)[n++] = carry; } } return n; } /* target base wbase < source base ubase */ static size_t _coeff_from_larger_base(mpd_t *w, size_t wlen, mpd_uint_t wbase, mpd_uint_t *u, mpd_ssize_t ulen, mpd_uint_t ubase, uint32_t *status) { size_t n = 0; assert(wlen > 0 && ulen > 0); assert(wbase < ubase); do { if (n >= wlen) { if (!mpd_qresize(w, n+1, status)) { return SIZE_MAX; } wlen = n+1; } w->data[n++] = (uint32_t)_mpd_shortdiv_b(u, u, ulen, wbase, ubase); /* ulen is at least 1. u[ulen-1] can only be zero if ulen == 1. */ ulen = _mpd_real_size(u, ulen); } while (u[ulen-1] != 0); return n; } #endif /* target base 'wbase' > source base 'ubase' */ static size_t _coeff_from_smaller_base(mpd_t *w, mpd_ssize_t wlen, mpd_uint_t wbase, const uint32_t *u, size_t ulen, mpd_uint_t ubase, uint32_t *status) { … } /* * Convert an integer mpd_t to a multiprecision integer with base <= 2**16. * The least significant word of the result is (*rdata)[0]. * * If rdata is NULL, space is allocated by the function and rlen is irrelevant. * In case of an error any allocated storage is freed and rdata is set back to * NULL. * * If rdata is non-NULL, it MUST be allocated by one of libmpdec's allocation * functions and rlen MUST be correct. If necessary, the function will resize * rdata. In case of an error the caller must free rdata. * * Return value: In case of success, the exact length of rdata, SIZE_MAX * otherwise. */ size_t mpd_qexport_u16(uint16_t **rdata, size_t rlen, uint32_t rbase, const mpd_t *src, uint32_t *status) { … } /* * Convert an integer mpd_t to a multiprecision integer with base<=UINT32_MAX. * The least significant word of the result is (*rdata)[0]. * * If rdata is NULL, space is allocated by the function and rlen is irrelevant. * In case of an error any allocated storage is freed and rdata is set back to * NULL. * * If rdata is non-NULL, it MUST be allocated by one of libmpdec's allocation * functions and rlen MUST be correct. If necessary, the function will resize * rdata. In case of an error the caller must free rdata. * * Return value: In case of success, the exact length of rdata, SIZE_MAX * otherwise. */ size_t mpd_qexport_u32(uint32_t **rdata, size_t rlen, uint32_t rbase, const mpd_t *src, uint32_t *status) { … } /* * Converts a multiprecision integer with base <= UINT16_MAX+1 to an mpd_t. * The least significant word of the source is srcdata[0]. */ void mpd_qimport_u16(mpd_t *result, const uint16_t *srcdata, size_t srclen, uint8_t srcsign, uint32_t srcbase, const mpd_context_t *ctx, uint32_t *status) { … } /* * Converts a multiprecision integer with base <= UINT32_MAX to an mpd_t. * The least significant word of the source is srcdata[0]. */ void mpd_qimport_u32(mpd_t *result, const uint32_t *srcdata, size_t srclen, uint8_t srcsign, uint32_t srcbase, const mpd_context_t *ctx, uint32_t *status) { … } /******************************************************************************/ /* From triple */ /******************************************************************************/ #if defined(CONFIG_64) && defined(__SIZEOF_INT128__) static mpd_ssize_t _set_coeff(uint64_t data[3], uint64_t hi, uint64_t lo) { … } #else static size_t _uint_from_u16(mpd_uint_t *w, mpd_ssize_t wlen, const uint16_t *u, size_t ulen) { const mpd_uint_t ubase = 1U<<16; mpd_ssize_t n = 0; mpd_uint_t carry; assert(wlen > 0 && ulen > 0); w[n++] = u[--ulen]; while (--ulen != SIZE_MAX) { carry = _mpd_shortmul_c(w, w, n, ubase); if (carry) { if (n >= wlen) { abort(); /* GCOV_NOT_REACHED */ } w[n++] = carry; } carry = _mpd_shortadd(w, n, u[ulen]); if (carry) { if (n >= wlen) { abort(); /* GCOV_NOT_REACHED */ } w[n++] = carry; } } return n; } static mpd_ssize_t _set_coeff(mpd_uint_t *data, mpd_ssize_t len, uint64_t hi, uint64_t lo) { uint16_t u16[8] = {0}; u16[7] = (uint16_t)((hi & 0xFFFF000000000000ULL) >> 48); u16[6] = (uint16_t)((hi & 0x0000FFFF00000000ULL) >> 32); u16[5] = (uint16_t)((hi & 0x00000000FFFF0000ULL) >> 16); u16[4] = (uint16_t) (hi & 0x000000000000FFFFULL); u16[3] = (uint16_t)((lo & 0xFFFF000000000000ULL) >> 48); u16[2] = (uint16_t)((lo & 0x0000FFFF00000000ULL) >> 32); u16[1] = (uint16_t)((lo & 0x00000000FFFF0000ULL) >> 16); u16[0] = (uint16_t) (lo & 0x000000000000FFFFULL); return (mpd_ssize_t)_uint_from_u16(data, len, u16, 8); } #endif static int _set_uint128_coeff_exp(mpd_t *result, uint64_t hi, uint64_t lo, mpd_ssize_t exp) { … } int mpd_from_uint128_triple(mpd_t *result, const mpd_uint128_triple_t *triple, uint32_t *status) { … } /******************************************************************************/ /* As triple */ /******************************************************************************/ #if defined(CONFIG_64) && defined(__SIZEOF_INT128__) static void _get_coeff(uint64_t *hi, uint64_t *lo, const mpd_t *a) { … } #else static size_t _uint_to_u16(uint16_t w[8], mpd_uint_t *u, mpd_ssize_t ulen) { const mpd_uint_t wbase = 1U<<16; size_t n = 0; assert(ulen > 0); do { if (n >= 8) { abort(); /* GCOV_NOT_REACHED */ } w[n++] = (uint16_t)_mpd_shortdiv(u, u, ulen, wbase); /* ulen is at least 1. u[ulen-1] can only be zero if ulen == 1. */ ulen = _mpd_real_size(u, ulen); } while (u[ulen-1] != 0); return n; } static void _get_coeff(uint64_t *hi, uint64_t *lo, const mpd_t *a) { uint16_t u16[8] = {0}; mpd_uint_t data[5] = {0}; switch (a->len) { case 5: data[4] = a->data[4]; /* fall through */ case 4: data[3] = a->data[3]; /* fall through */ case 3: data[2] = a->data[2]; /* fall through */ case 2: data[1] = a->data[1]; /* fall through */ case 1: data[0] = a->data[0]; break; default: abort(); /* GCOV_NOT_REACHED */ } _uint_to_u16(u16, data, a->len); *hi = (uint64_t)u16[7] << 48; *hi |= (uint64_t)u16[6] << 32; *hi |= (uint64_t)u16[5] << 16; *hi |= (uint64_t)u16[4]; *lo = (uint64_t)u16[3] << 48; *lo |= (uint64_t)u16[2] << 32; *lo |= (uint64_t)u16[1] << 16; *lo |= (uint64_t)u16[0]; } #endif static enum mpd_triple_class _coeff_as_uint128(uint64_t *hi, uint64_t *lo, const mpd_t *a) { … } mpd_uint128_triple_t mpd_as_uint128_triple(const mpd_t *a) { … }
Codebase Browser
cpython