4 GNU MPC Basics

All declarations needed to use GNU MPC are collected in the include file mpc.h. It is designed to work with both C and C++ compilers. You should include that file in any program using the GNU MPC library by adding the line

   #include "mpc.h"

To use, for instance, input and output functions, this line should come after the inclusion of standard C headers such as stdio.h.

4.1 Nomenclature and Types

A complex number is a pair of two arbitrary precision floating-point numbers of type mpfr_t (for the real and imaginary parts). The C data type for such objects is mpc_t. As for GNU MP and GNU MPFR, a variable of type mpc_t is in fact an array with one struct element, which is allocated on the stack (while the mantissas of its real and imaginary parts will eventually be allocated on the heap, usually by a call to mpc_init2, see Initialization Functions). The type mpc_ptr represents a pointer to the corresponding struct, while mpc_srcptr represents a constant pointer.

The precision is the number of bits used to represent the mantissa of the real and imaginary parts; the corresponding C data type is mpfr_prec_t. For more details on the allowed precision range, see Section “Nomenclature and Types” in GNU MPFR. While for most applications it should be appropriate to choose the same precision for the real and imaginary parts, it is possible to create numbers with different precisions for the two parts.

The rounding mode specifies the way to round the result of a complex operation, in case the exact result can not be represented exactly in the destination mantissa; the corresponding C data type is mpc_rnd_t. A complex rounding mode is a pair of two rounding modes: one for the real part, one for the imaginary part.

4.2 Function Classes

The main class of functions in the GNU MPC library consists of functions for complex arithmetic. Their names begin with mpc_, and they operate (mainly) on arguments of type mpc_t (or mpc_ptr and mpc_srcptr to be more precise; for the sake of simplicity, we will use mpc_t in the following, the interested reader may consult the include file mpc.h for full details).

A new class of experimental functions, prefixed with mpcb_, operates on complex balls, and functions starting with mpcr_ operate on their corresponding radii, see Ball Arithmetic.

4.3 GNU MPC Variable Conventions

As a general rule, all GNU MPC functions expect output arguments before input arguments. This notation is based on an analogy with the assignment operator.

GNU MPC allows you to use the same variable for both input and output in the same expression. For example, the main function for floating-point multiplication, mpc_mul, can be used like this: mpc_mul (x, x, x, rnd_mode). This computes the square of x with rounding mode rnd_mode and puts the result back in x.

Before you can assign to an GNU MPC variable, you need to initialise it by calling one of the special initialization functions. When you are done with a variable, you need to clear it out, using one of the functions for that purpose.

A variable should only be initialised once, or at least cleared out between each initialization. After a variable has been initialised, it may be assigned to any number of times.

For efficiency reasons, avoid to initialise and clear out a variable in loops. Instead, initialise it before entering the loop, and clear it out after the loop has exited.

You do not need to be concerned about allocating additional space for GNU MPC variables, since each of its real and imaginary part has a mantissa of fixed size. Hence unless you change its precision, or clear and reinitialise it, a complex variable will have the same allocated space during all its life.

4.4 Rounding Modes

A complex rounding mode is of the form MPC_RNDxy where x and y are one of N (to nearest), Z (towards zero), U (towards plus infinity), D (towards minus infinity), A (away from zero, that is, towards plus or minus infinity depending on the sign of the number to be rounded). The first letter refers to the rounding mode for the real part, and the second one for the imaginary part. For example MPC_RNDZU indicates to round the real part towards zero, and the imaginary part towards plus infinity.

The ‘round to nearest’ mode works as in the IEEE 754 standard: in case the number to be rounded lies exactly in the middle of two representable numbers, it is rounded to the one with the least significant bit set to zero. For example, the number 5, which is represented by (101) in binary, is rounded to (100)=4 with a precision of two bits, and not to (110)=6.

A useful macro is MPC_RND(r1,r2), which from two MPFR rounding modes r1 for the real part and r2 for the imaginary part, constructs the corresponding MPC rounding mode. Similarly, for an MPC rounding mode r, MPC_RND_RE(r) and MPC_RND_IM(r) yield the corresponding MPFR rounding modes for the real and imaginary part, respectively.

4.5 Semantics

Finite precision complex floating-point numbers are in reality elements of Q(i); more precisely, their real and imaginary parts are rational numbers with denominator a power of 2 (such a rational number is also called a dyadic rational, and we may call its complex counterpart a dyadic number). All inputs to functions are thus considered exact, and results of function evaluations are mathematically defined as complex numbers. We consider their real and imaginary parts separately, that is, a function with n complex inputs and one complex result is identified with a function with 2n real inputs and two real results. If such a result is itself a dyadic number representable at the precision of the output variable, then GNU MPC assigns this exact result; otherwise the mathematical result falls between two dyadic numbers that are representable in the target precision, and GNU MPC assigns one of these numbers to the output variable according to the rounding mode. Notice that this uniquely defines the expected output of any function. Since every variable carries its own precision (more precisely, two precisions, for the real and the imaginary part), function calls are made without an additional argument for the target precision; the rounding mode, however, needs in general to be given as an additional argument. Notice that in the case of an inexact output, there is a rounding error of at most the difference between two consecutive representable dyadic numbers. When such an output is used as an input to another function, it is again considered exact, and the GNU MPC library does not track how rounding errors propagate through successive function applications (unless ball arithmetic is used).

This semantics is a generalization of the IEEE-754 standard (also included into recent C standards) for real floating-point numbers of fixed precision to complex floating-point numbers of arbitrary precision. The goal of the GNU MPC library is to provide operations with complex numbers that are well-defined with predictable results.

Since the real and imaginary parts of complex numbers are mpfr_t, besides finite values, they may also contain special values from the GNU MPFR library such as NaN and (signed) infinities; additionally zeroes also have signs. Since GNU MPFR does not distinguish between quiet and signaling NaN, this holds for GNU MPC too. These are handled analogously to GNU MPFR; for special cases, see Special Values, Signed Zeroes and Branch Cuts.

4.6 Return Value

Most GNU MPC functions have a return value of type int, which is used to indicate the position of the rounded real and imaginary parts with respect to the exact (infinite precision) values. If this integer is k, the macros MPC_INEX_RE(k) and MPC_INEX_IM(k) give 0 if the corresponding rounded value is exact, a negative value if the rounded value is less than the exact one, and a positive value if it is greater than the exact one. Similarly, functions computing a result of type mpfr_t return an integer that is 0, positive or negative depending on whether the rounded value is the same, larger or smaller then the exact result.

Some functions, such as mpc_sin_cos, compute two complex results; the macros MPC_INEX1(k) and MPC_INEX2(k), applied to the return value k of such a function, yield the exactness value corresponding to the first and the second computed value, respectively.


4.7 Special Values, Signed Zeroes and Branch Cuts

When evaluating a function in an argument whose real or imaginary part is a special value (that is NaN or a (signed) infinity) or a (signed) zero, our goal is to compute the same result as specified for the corresponding function in the C standard.

More precisely, we implement the following principles for real numbers (which appear to be consistent with the C standard). The value +Inf is considered as the limit of sequences of positive numbers tending to infinity, -Inf as the limit of sequences of negative numbers with absolute values tending to infinity. To handle division by infinities it is then natural to also consider signed zeroes as limits of sequences, namely +0 and -0 as the limits of sequences of positive resp. negative numbers with absolute value tending to zero, so that we obtain for the division of two real numbers that 1/+Inf=+0, 1/-Inf=-0, 1/+0=+Inf and 1/-0=-Inf. On the other hand, +0/+0=NaN and +Inf/+Inf=NaN, since any value is possible depending on how fast the dividend and divisor sequences, respectively, tend to zero or infinity. There is only one, unsigned, NaN value representing not a number; and any formula involving a NaN results in a NaN.

A peculiarity arises when adding signed zeroes; clearly, (+0)+(+0)=+0 and (-0)+(-0)=-0, but what should the result of (+0)+(-0) be? If each operand is replaced by a sequence of positive resp. negative real numbers tending to zero, then the sequence of their sums may or may not have a common sign (it could even be alternating between positive and negative). A logical result would be an unsigned zero, but this is not foreseen by the IEEE-754 standard since the binary format reserves one bit for the sign, so that every number has either a plus or a minus sign (which is ignored for NaN). So in case the sign of a resulting zero cannot be determined, the result is arbitrarily fixed to be +0. One could imagine a fourth type of zero, the exact zero; such a number is also not foreseen by the standard, so that any zero input to a function is considered signed, and as for other numbers, “inexact” output (here, an arbitrarily assigned plus sign) is considered “exact” as an input to a subsequent function call. Beware that although +0 and -0 differ in sign, they nevertheless compare as equal; so equal function inputs may lead to differing results!

Since complex numbers are given through their real and imaginary parts, there are in fact four different complex zeroes depending on the signs (all of which compare as equal).

Signs in parts of complex numbers become meaningful when branch cuts come into play. The complex logarithm, for instance, is defined up to integer multiples of 2 Pi i, since exp(2 Pi i) = 1. More precisely, when trying to define the logarithm as a continuous function, one walk on a circle around the origin leads to values that differ by 2 Pi i. So there is necessarily a discontinuity, which is placed along a branch cut. In GNU MPC, the chosen branch cuts are the same as those specified for the corresponding functions in the C standard. For instance, the complex logarithm has a branch cut along the negative real axis (which in particular leads to results coinciding with GNU MPFR on the positive real axis). Since -1+0i can be seen as above the branch cut, but -1-0i can be seen as below the branch cut, we compute the exact results log(-1+0i)=+0+Pi i and log(-1-0i)=+0-Pi i (and round to a representable dyadic number), where the imaginary parts differ by 2 Pi i and the signs of the real parts cannot be determined and are thus set to plus. (In fact extending the preceding arguments, we can also consider the exact finite value -1 as the limit of a sequence tending to -1; then the real part of the result tends to zero and is positive when -1 is approached from below on the real axis and negative when -1 is approached from above on the real axis.)

This semantics becomes fuzzy when mixed operations between complex and real numbers enter the game. Two different approaches stand out. First of all, a real value may first be promoted to a complex value, and then the previous rules apply. For instance, with this approach, (1-0i)+(+0)=(1-0i)+(+0+0i)=1+0i, where the sign of the imaginary part of the result cannot be determined and is thus set to plus. Alternatively, the zero imaginary part of the real number may be interpreted as an exact zero; then (1-0i)+(+0)=1-0i. The C standard appears to follow the second approach, which we thus also adopt in MPC.

4.8 Exceptions

This version of GNU MPC does not support exceptions (invalid, division by zero, overflow, underflow, inexact). A future version might support them, through the GNU MPFR mechanism for supporting exceptions.