Arithmetic Functions

Basic Arithmetic

Most of these functions are part of the specification, and some descriptions contain direct quotes.

add

void mpd_qadd(mpd_t *result, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_add(mpd_t *result, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to a + b.

void mpd_qadd_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                    const mpd_context_t *ctx, uint32_t *status);
void mpd_qadd_i32(mpd_t *result, const mpd_t *a, int32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qadd_i64(mpd_t *result, const mpd_t *a, int64_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qadd_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                   const mpd_context_t *ctx, uint32_t *status);
void mpd_qadd_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qadd_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                  const mpd_context_t *ctx, uint32_t *status);

void mpd_add_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                   mpd_context_t *ctx);
void mpd_add_i32(mpd_t *result, const mpd_t *a, int32_t b,
                 mpd_context_t *ctx);
void mpd_add_i64(mpd_t *result, const mpd_t *a, int64_t b,
                 mpd_context_t *ctx);
void mpd_add_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                  mpd_context_t *ctx);
void mpd_add_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                 mpd_context_t *ctx);
void mpd_add_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                 mpd_context_t *ctx);

Set result to a + b.

Changed in version 2.4.0: All functions are now available in both the 64-bit and 32-bit builds.

subtract

void mpd_qsub(mpd_t *result, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_sub(mpd_t *result, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to a - b.

void mpd_qsub_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                    const mpd_context_t *ctx, uint32_t *status);
void mpd_qsub_i32(mpd_t *result, const mpd_t *a, int32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qsub_i64(mpd_t *result, const mpd_t *a, int64_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qsub_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                   const mpd_context_t *ctx, uint32_t *status);
void mpd_qsub_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qsub_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                  const mpd_context_t *ctx, uint32_t *status);

void mpd_sub_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                   mpd_context_t *ctx);
void mpd_sub_i32(mpd_t *result, const mpd_t *a, int32_t b,
                 mpd_context_t *ctx);
void mpd_sub_i64(mpd_t *result, const mpd_t *a, int64_t b,
                 mpd_context_t *ctx);
void mpd_sub_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                  mpd_context_t *ctx);
void mpd_sub_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                 mpd_context_t *ctx);
void mpd_sub_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                 mpd_context_t *ctx);

Set result to a - b.

Changed in version 2.4.0: All functions are now available in both the 64-bit and 32-bit builds.

multiply

void mpd_qmul(mpd_t *result, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_mul(mpd_t *result, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to a * b.

void mpd_qmul_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_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_i64(mpd_t *result, const mpd_t *a, int64_t b,
                  const mpd_context_t *ctx, uint32_t *status);
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_u32(mpd_t *result, const mpd_t *a, uint32_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);

void mpd_mul_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                   mpd_context_t *ctx);
void mpd_mul_i32(mpd_t *result, const mpd_t *a, int32_t b,
                 mpd_context_t *ctx);
void mpd_mul_i64(mpd_t *result, const mpd_t *a, int64_t b,
                 mpd_context_t *ctx);
void mpd_mul_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                  mpd_context_t *ctx);
void mpd_mul_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                 mpd_context_t *ctx);
void mpd_mul_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                 mpd_context_t *ctx);

Set result to a * b.

Changed in version 2.4.0: All functions are now available in both the 64-bit and 32-bit builds.

divide

void mpd_qdiv(mpd_t *q, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_div(mpd_t *q, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to a / b.

Changed in version 2.5.0: If MPD_Malloc_error occurs for extremely large precisions, the function now retries the operation in case the result is exact and fits in a much lower precision.

void mpd_qdiv_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                    const mpd_context_t *ctx, uint32_t *status);
void mpd_qdiv_i32(mpd_t *result, const mpd_t *a, int32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qdiv_i64(mpd_t *result, const mpd_t *a, int64_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qdiv_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                   const mpd_context_t *ctx, uint32_t *status);
void mpd_qdiv_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_qdiv_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                  const mpd_context_t *ctx, uint32_t *status);

void mpd_div_ssize(mpd_t *result, const mpd_t *a, mpd_ssize_t b,
                   mpd_context_t *ctx);
void mpd_div_i32(mpd_t *result, const mpd_t *a, int32_t b,
                 mpd_context_t *ctx);
void mpd_div_i64(mpd_t *result, const mpd_t *a, int64_t b,
                 mpd_context_t *ctx);
void mpd_div_uint(mpd_t *result, const mpd_t *a, mpd_uint_t b,
                  mpd_context_t *ctx);
void mpd_div_u32(mpd_t *result, const mpd_t *a, uint32_t b,
                 mpd_context_t *ctx);
void mpd_div_u64(mpd_t *result, const mpd_t *a, uint64_t b,
                 mpd_context_t *ctx);

Set result to a / b.

Changed in version 2.4.0: All functions are now available in both the 64-bit and 32-bit builds.

fused-multiply-add

void mpd_qfma(mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_t *c,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_fma(mpd_t *r, const mpd_t *a, const mpd_t *b, const mpd_t *c,
             mpd_context_t *ctx);

Set r to a * b + c with a single, final rounding.

Integer Division

divide-integer

void mpd_qdivint(mpd_t *q, const mpd_t *a, const mpd_t *b,
                 const mpd_context_t *ctx, uint32_t *status);
void mpd_divint(mpd_t *q, const mpd_t *a, const mpd_t *b,
                mpd_context_t *ctx);

Set q to the integer part of a / b.

remainder

void mpd_qrem(mpd_t *r, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_rem(mpd_t *r, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set r to the remainder of a / b.

remainder-near

void mpd_qrem_near(mpd_t *r, const mpd_t *a, const mpd_t *b,
                   const mpd_context_t *ctx, uint32_t *status);
void mpd_rem_near(mpd_t *r, const mpd_t *a, const mpd_t *b,
                  mpd_context_t *ctx);

Set r to a - b * n, where n is the integer nearest the exact value of a / b. If two integers are equally near then the even one is chosen.

divmod

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_divmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b,
                mpd_context_t *ctx);

Set q and r to the integer part and remainder of a / b. This function is not part of the standard.

Powers and Logarithms

exp

void mpd_qexp(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
              uint32_t *status);
void mpd_exp(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to e ** a. Results are correctly rounded if the allcr field in the context is set to 1.

ln

void mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
             uint32_t *status);
void mpd_ln(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the natural logarithm of a. Results are correctly rounded if the allcr field in the context is set to 1.

Changed in version 2.4.0: This function is now thread-safe.

log10

void mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                uint32_t *status);
void mpd_log10(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the base-10 logarithm of a. Results are correctly rounded if the allcr field in the context is set to 1.

Changed in version 2.4.0: This function is now thread-safe.

power

void mpd_qpow(mpd_t *result, const mpd_t *base, const mpd_t *exp,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_pow(mpd_t *result, const mpd_t *base, const mpd_t *exp,
             mpd_context_t *ctx);

Set result to base ** exp. Integer powers are exact, provided that that the result is finite and fits into prec.

Results are not correctly rounded, even if the allcr context field is set to 1. This might change in a future release. The error of the function is less than 0.5ULP+t (1ULP+t for directed rounding modes), where t has a maximum of 0.1ULP, but is almost always less than 0.01ULP.

Changed in version 2.4.0: This function is now thread-safe.

square-root

void mpd_qsqrt(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
               uint32_t *status);
void mpd_sqrt(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the square root of a. This function is always correctly rounded and always uses the ROUND_HALF_EVEN mode.

Changed in version 2.5.0: If MPD_Malloc_error occurs for extremely large precisions, the function now retries the operation in case the result is exact and fits in a much lower precision.

inverse-square-root

void mpd_qinvroot(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                  uint32_t *status);
void mpd_invroot(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the reciprocal of the square root of a. The function always uses the ROUND_HALF_EVEN mode. Results are not correctly rounded, even if the allcr context field is set to 1. This might change in a future release.

Sign and Absolute

minus

void mpd_qminus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                uint32_t *status);
void mpd_minus(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to - a. The operation is context sensitive.

plus

void mpd_qplus(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
               uint32_t *status);
void mpd_plus(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to + a. The operation is context sensitive.

abs

void mpd_qabs(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
              uint32_t *status);
void mpd_abs(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the absolute value of a. For negative numbers, mpd_qabs is the same as mpd_qminus, for positive numbers the same as mpd_qplus.

Comparisons

compare

int mpd_qcompare(mpd_t *result, const mpd_t *a, const mpd_t *b,
                 const mpd_context_t *ctx, uint32_t *status);
int mpd_compare(mpd_t *result, const mpd_t *a, const mpd_t *b,
                mpd_context_t *ctx);

Set result to -1 if a is less than b, 0 if a is equal to b and 1 if a is greater than b. For convenience, the result is also returned as a C integer. The MPD_Invalid_operation condition is added to status if at least one of the operands is a signaling NaN. In this case, result is set to NaN and INT_MAX is returned.

compare-signal

int mpd_qcompare_signal(mpd_t *result, const mpd_t *a, const mpd_t *b,
                        const mpd_context_t *ctx, uint32_t *status);
int mpd_compare_signal(mpd_t *result, const mpd_t *a, const mpd_t *b,
                       mpd_context_t *ctx);

Same as mpd_qcompare, except that the MPD_Invalid_operation condition is also added to status if at least one of the operands is a quiet NaN.

compare-total

int mpd_compare_total(mpd_t *result, const mpd_t *a, const mpd_t *b);

Compare a and b using a total ordering. This function is always quiet.

compare-total-magnitude

int mpd_compare_total_mag(mpd_t *result, const mpd_t *a, const mpd_t *b);

Compare the magnitude of a and b using a total ordering. This function is always quiet.

cmp

int mpd_qcmp(const mpd_t *a, const mpd_t *b, uint32_t *status);
int mpd_cmp(const mpd_t *a, const mpd_t *b, mpd_context_t *ctx);

Same as mpd_qcompare_signal, but without the decimal result operand.

cmp-total

int mpd_cmp_total(const mpd_t *a, const mpd_t *b);

Same as mpd_compare_total, but without the decimal result operand.

cmp-total-magnitude

int mpd_cmp_total_mag(const mpd_t *a, const mpd_t *b);

Same as mpd_compare_total_mag, but without the decimal result operand.

max

void mpd_qmax(mpd_t *result, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_max(mpd_t *result, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to the maximum of a and b. If one of the operands is a quiet NaN and the other is numeric, the numeric operand is returned.

max-mag

void mpd_qmax_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_max_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
                 mpd_context_t *ctx);

Same as mpd_qmax, but the numerical comparison takes place with the signs of the operands set to 0 (positive).

min

void mpd_qmin(mpd_t *result, const mpd_t *a, const mpd_t *b,
              const mpd_context_t *ctx, uint32_t *status);
void mpd_min(mpd_t *result, const mpd_t *a, const mpd_t *b,
             mpd_context_t *ctx);

Set result to the minimum of a and b. If one of the operands is a quiet NaN and the other is numeric, the numeric operand is returned.

min-mag

void mpd_qmin_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_min_mag(mpd_t *result, const mpd_t *a, const mpd_t *b,
                 mpd_context_t *ctx);

Same as mpd_qmin, but the numerical comparison takes place with the signs of the operands set to 0 (positive).

Closest Numbers

next-minus

void mpd_qnext_minus(mpd_t *result, const mpd_t *a,
                     const mpd_context_t *ctx, uint32_t *status);
void mpd_next_minus(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

The closest representable number that is smaller than a.

next-plus

void mpd_qnext_plus(mpd_t *result, const mpd_t *a,
                    const mpd_context_t *ctx, uint32_t *status);
void mpd_next_plus(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

The closest representable number that is larger than a.

next-toward

void mpd_qnext_toward(mpd_t *result, const mpd_t *a, const mpd_t *b,
                      const mpd_context_t *ctx, uint32_t *status);
void mpd_next_toward(mpd_t *result, const mpd_t *a, const mpd_t *b,
                     mpd_context_t *ctx);

Representable number closest to a that is in the direction towards b.

Quantizing and Normalizing

quantize

void mpd_qquantize(mpd_t *result, const mpd_t *a, const mpd_t *b,
                   const mpd_context_t *ctx, uint32_t *status);
void mpd_quantize(mpd_t *result, const mpd_t *a, const mpd_t *b,
                  mpd_context_t *ctx);

Set result to the number that is equal in value to a, but has the exponent of b.

rescale

void mpd_qrescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp,
                  const mpd_context_t *ctx, uint32_t *status);
void mpd_rescale(mpd_t *result, const mpd_t *a, mpd_ssize_t exp,
                 mpd_context_t *ctx);

Set result to the number that is equal in value to a, but has the exponent exp. Special numbers are copied without signaling. This function is not part of the standard. It is also not equivalent to the rescale function that was removed from the standard.

same-quantum

int mpd_same_quantum(const mpd_t *a, const mpd_t *b);

Return 1 if a and b have the same exponent, 0 otherwise.

reduce

void mpd_qreduce(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                 uint32_t *status);
void mpd_reduce(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

If a is finite after applying rounding and overflow/underflow checks, result is set to the simplest form of a with all trailing zeros removed.

Integral Values

round-to-integral-exact

void mpd_qround_to_intx(mpd_t *result, const mpd_t *a,
                        const mpd_context_t *ctx, uint32_t *status);
void mpd_round_to_intx(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Round to an integer, using the rounding mode of the context. Only a signaling NaN causes an MPD_Invalid_operation condition.

round-to-integral-value

void mpd_qround_to_int(mpd_t *result, const mpd_t *a,
                       const mpd_context_t *ctx, uint32_t *status);
void mpd_round_to_int(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Same as mpd_qround_to_intx, but the MPD_Inexact and MPD_Rounded flags are never set.

floor

void mpd_qfloor(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                uint32_t *status);
void mpd_floor(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the floor of a. Not part of the standard.

Changed in version 2.5.0: The function now sets MPD_Invalid_operation if a is special. Previously it followed mpd_round_to_int, for which the specification allows special values.

ceiling

void mpd_qceil(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
               uint32_t *status);
void mpd_ceil(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

Set result to the ceiling of a. Not part of the standard.

Changed in version 2.5.0: The function now sets MPD_Invalid_operation if a is special. Previously it followed mpd_round_to_int, for which the specification allows special values.

truncate

void mpd_qtrunc(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
                uint32_t *status);
void mpd_trunc(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx);

Set result to the truncated value of a. Not part of the standard.

Changed in version 2.5.0: The function now sets MPD_Invalid_operation if a is special. Previously it followed mpd_round_to_int, for which the specification allows special values.

Scale

logb

void mpd_qlogb(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
               uint32_t *status);
void mpd_logb(mpd_t *result, const mpd_t *a, mpd_context_t *ctx);

If a is NaN, the general rules apply. If a is infinite, result is +Infinity. If a is zero, result is -Infinity and MPD_Invalid_operation is added to status. Otherwise, result is the same as the adjusted exponent of a, or floor(log10(a)).

scaleb

void mpd_qscaleb(mpd_t *result, const mpd_t *a, const mpd_t *b,
                 const mpd_context_t *ctx, uint32_t *status);
void mpd_scaleb(mpd_t *result, const mpd_t *a, const mpd_t *b,
                mpd_context_t *ctx);

b must be an integer with exponent 0. If a is infinite, result is set to a. Otherwise, result is a with the value of b added to the exponent.

Integer Functions

powmod

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_powmod(mpd_t *result, const mpd_t *base, const mpd_t *exp,
                const mpd_t *mod, mpd_context_t *ctx);

Set result to (base ** exp) % mod. All operands must be integers. The function fails if result does not fit in the current prec.