Constant gmp_mpfr_sys::C::MPC::Ball_Arithmetic

source · 
pub const Ball_Arithmetic: ();
Expand description

This constant is a place-holder for documentation; do not use it in code.

Next: , Previous: , Up: GNU MPC   [Index]

6 Ball Arithmetic

Since release 1.3.0, GNU MPC contains a simple and very limited implementation of complex balls (or rather, circles). This part is experimental, its interface may vary and it may be removed completely in future releases.

A complex ball of the new type mpcb_t is defined by a non-zero centre c of the type mpc_t and a relative radius r of the new type mpcr_t, and it represents all complex numbers z = c (1 + ϑ) with |ϑ| ≤ r, or equivalently the closed circle with centre c and radius r |c|. The approach of using a relative error (or radius) instead of an absolute one simplifies error analyses for multiplicative operations (multiplication, division, square roots, and the AGM), at the expense of making them more complicated for additive operations. It has the major drawback of not being able to represent balls centred at 0; in floating point arithmetic, however, 0 is never reached by rounding, but only through operations with exact result, which could be handled at a higher, application level. For more discussion on these issues, see the file algorithms.tex.

6.1 Radius type and functions

The radius type is defined by 

struct {
   int64_t mant;
   int64_t exp;
}

with the usual trick in the GNU multiprecision libraries of defining the main type mpcr_t as a 1-dimensional array of this struct, and variable and constant pointers mpcr_ptr and mpcr_srcptr. It can contain the special values infinity or zero, or floating point numbers encoded as m⋅2e for a positive mantissa m and an arbitrary (usually negative) exponent e. Normalised finite radii use 31 bits for the mantissa, that is, 230≤m≤231 - 1. The special values infinity and 0 are encoded through the sign of m, but should be tested for and set using dedicated functions.

Unless indicated otherwise, the following functions assume radius arguments to be normalised, they return normalised results, and they round their results up, not necessarily to the smallest representable number, although reasonable effort is made to get a tight upper bound: They only guarantee that their outputs are an upper bound on the true results. (There may be a trade-off between tightness of the result and speed of computation. For instance, when a 32-bit mantissa is normalised, an even mantissa should be divided by 2, an odd mantissa should be divided by 2 and 1 should be added, and then in both cases the exponent must be increased by 1. It might be more efficient to add 1 all the time instead of testing the last bit of the mantissa.)

Function: int mpcr_inf_p (mpcr_srcptr r)
Function: int mpcr_zero_p (mpcr_srcptr r)

Test whether r is infinity or zero, respectively, and return a boolean.

Function: int mpcr_lt_half_p (mpcr_srcptr r)

Return true if r<1/2, and false otherwise. (Everywhere in this document, true means any non-zero value, and false means zero.)

Function: int mpcr_cmp (mpcr_srcptr r, mpcr_srcptr s)

Return +1, 0 or -1 depending on whether r is larger than, equal to or less than s, with the natural total order on the compactified non-negative real axis letting 0 be smaller and letting infinity be larger than any finite real number.

Function: void mpcr_set_inf (mpcr_ptr r)
Function: void mpcr_set_zero (mpcr_ptr r)
Function: void mpcr_set_one (mpcr_ptr r)
Function: void mpcr_set (mpcr_ptr r, mpcr_srcptr s)
Function: void mpcr_set_ui64_2si64 (mpcr_ptr r, uint64_t mant, int64_t exp)

Set r to infinity, zero, 1, s or mant⋅2exp, respectively.

Function: void mpcr_max (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t)

Set r to the maximum of s and t.

Function: int64_t mpcr_get_exp (mpcr_srcptr r)

Assuming that r is neither infinity nor 0, return its exponent e when writing r = m⋅2e with 1/2 ≤ m < 1. (Notice that this is not the same as the field exp in the struct representing a radius, but that instead it is independent of the implementation.) Otherwise the behaviour is undefined.

Function: void mpcr_out_str (FILE *f, mpcr_srcptr r)

Output r on f, which may be stdout. Caveat: This function so far serves mainly for debugging purposes, its behaviour will probably change in the future.

Function: void mpcr_add (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t)
Function: void mpcr_sub (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t)
Function: void mpcr_mul (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t)
Function: void mpcr_div (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t)
Function: void mpcr_mul_2ui (mpcr_ptr r, mpcr_srcptr s, unsigned long int t)
Function: void mpcr_div_2ui (mpcr_ptr r, mpcr_srcptr s, unsigned long int t)
Function: void mpcr_sqr (mpcr_ptr r, mpcr_srcptr s)
Function: void mpcr_sqrt (mpcr_ptr r, mpcr_srcptr s)

Set r to the sum, difference, product or quotient of s and t, or to the product of s by 2t or to the quotient of s by 2t, or to the square or the square root of s. If any of the arguments is infinity, or if a difference is negative, the result is infinity.

Function: void mpcr_sub_rnd (mpcr_ptr r, mpcr_srcptr s, mpcr_srcptr t, mpfr_rnd_t rnd)

Set r to the difference of s and t, rounded into direction rnd, which can be one of MPFR_RNDU or MPFR_RNDD. If one of the arguments is infinity or the difference is negative, the result is infinity. Calling the function with MPFR_RNDU is equivalent to calling mpcr_sub.

This is one out of several functions taking a rounding parameter. Rounding down may be useful to obtain an upper bound when dividing by the result.

Function: void mpcr_c_abs_rnd (mpcr_ptr r, mpc_srcptr z, mpfr_rnd_t rnd)

Set r to the absolute value of the complex number z, rounded in direction rnd, which may be one of MPFR_RNDU or MPFR_RNDD.

Function: void mpcr_add_rounding_error (mpcr_ptr r, mpfr_prec_t p, mpfr_rnd_t rnd)

Set r to r + (1 + r) 2-p if rnd equals MPFR_RNDN, and to r + (1 + r) 21-p otherwise. The idea is that if a (potentially not representable) centre of an ideal complex ball of radius r is rounded to a representable complex number at precision p, this shifts the centre by up to 1/2 ulp (for rounding to nearest) or 1 ulp (for directed rounding of at least one of the real or imaginary parts), which increases the radius accordingly. So this function is typically called internally at the end of each operation with complex balls to account for the error made by rounding the centre.

6.2 Ball type and functions

The ball type is defined by 

typedef struct {
  mpc_t  c;
  mpcr_t r;
}

or, more precisely, mpcb_t is again a 1-dimensional array of such a struct, and variable and constant pointer types are defined as mpcb_ptr and mpcb_srcptr, respectively. As usual, the components should only be accessed through corresponding functions.

To understand functions on balls, one needs to consider the balls passed as arguments as sets of complex values, to which a mathematical function is applied; the C function “rounds up” in the sense that it returns a ball containing all possible values of the function in all the possible input values. Reasonable effort is made to return small balls, but again there is no guarantee that the result is the smallest possible one. In the current implementation, the centre of a ball returned as a value is obtained by applying the function to the centres of the balls passed as arguments, and rounding. While this is a natural approach, it is not the only possible one; however, it also simplifies the error analysis as already carried out for functions with regular complex arguments. Whenever the centre of a complex ball has a non-finite real or imaginary part (positive or negative infinity or NaN) the radius is set to infinity; this can be interpreted as the “useless ball”, representing the whole complex plane, whatever the value of the centre is.

Unlike for variables of mpc_t type, where the precision needs to be set explicitly at initialisation, variables of type mpcb_t handle their precision dynamically. Ball centres always have the same precision for their real and their imaginary parts (again this is a choice of the implementation; if they are of very different sizes, one could theoretically reduce the precision of the part that is smaller in absolute value, which is more strongly affected by the common error coded in the radius). When setting a complex ball from a value of a different type, an additional precision parameter is passed, which determines the precision of the centre. Functions on complex balls set the precision of their result depending on the input. In the current implementation, this is the minimum of the argument precisions, so if all balls are initially set to the same precision, this is preserved throughout the computations. (Notice that the exponent of the radius encodes roughly the number of correct binary digits of the ball centre; so it would also make sense to reduce the precision if the radius becomes larger.)

The following functions on complex balls are currently available; the eclectic collection is motivated by the desire to provide an implementation of the arithmetic-geometric mean of complex numbers through the use of ball arithmetic. As for functions taking complex arguments, there may be arbitrary overlaps between variables representing arguments and results; for instance mpcb_mul (z, z, z) is an allowed way of replacing the ball z by its square.

Function: void mpcb_init (mpcb_ptr z)
Function: void mpcb_clear (mpcb_ptr z)

Initialise or free memory for z; mpcb_init must be called once before using a variable, and mpcb_clear must be called once before stopping to use a variable. Unlike its mpc_t counterpart, mpcb_init does not fix the precision of z, but it sets its radius to infinity, so that z represents the whole complex plane.

Function: mpfr_prec_t mpcb_get_prec (mpcb_srcptr z)

Return the (common) precision of the real and the complex parts of the centre of z.

Function: void mpcb_set (mpcb_ptr z, mpcb_srcptr z1)

Set z to z1, preserving the precision of the centre.

Function: void mpcb_set_inf (mpcb_ptr z)

Set z to the whole complex plane. This is intended to be used much in the spirit of an assertion: When a precondition is not satisfied inside a function, it can set its result to this value, which will propagate through further computations.

Function: void mpcb_set_c (mpcb_ptr z, mpc_srcptr c, mpfr_prec_t prec, unsigned long int err_re, unsigned long int err_im)

Set z to a ball with centre c at precision prec. If prec is at least the maximum of the precisions of the real and the imaginary parts of c and err_re and err_im are 0, then the resulting ball is exact with radius zero. Using a larger value for prec makes sense if c is considered exact and a larger target precision for the result is desired, or some leeway for the working precision is to be taken into account. If prec is less than the precision of c, then usually some rounding error occurs when setting the centre, which is taken into account in the radius.

If err_re and err_im are non-zero, the argument c is considered as an inexact complex number, with a bound on the absolute error of its real part given in err_re as a multiple of 1/2 ulp of the real part of c, and a bound on the absolute error of its imaginary part given in err_im as a multiple of 1/2 ulp of the imaginary part of c. (Notice that if the parts of c have different precisions or exponents, the absolute values of their ulp differ.) Then z is created as a ball with centre c and a radius taking these errors on c as well as the potential additional rounding error for the centre into account. If the real part of c is 0, then err_re must be 0, since ulp of 0 makes no sense; otherwise the radius is set to infinity. The same remark holds for the imaginary part.

Using err_re and err_im different from 0 is particularly useful in two settings: If c is itself the result of a call to an mpc_ function with exact input and rounding mode MPC_RNDNN of both parts to nearest, then its parts are known with errors of at most 1/2 ulp, and setting err_re and err_im to 1 yields a ball which is known to contain the exact result (this motivates the strange unit of 1/2 ulp); if directed rounding was used, err_re and err_im can be set to 2 instead.

And if c is the result of a sequence of calls to mpc_ functions for which some error analysis has been carried out (as is frequently the case internally when implementing complex functions), again the resulting ball z is known to contain the exact result when using appropriate values for err_re and err_im.

Function: void mpcb_set_ui_ui (mpcb_ptr z, unsigned long int re, unsigned long int im, mpfr_prec_t prec)

Set z to a ball with centre re+I*im at precision prec or the size of an unsigned long int, whatever is larger.

Function: void mpcb_neg (mpcb_ptr z, mpcb_srcptr z1)
Function: void mpcb_add (mpcb_ptr z, mpcb_srcptr z1, mpcb_srcptr z2)
Function: void mpcb_mul (mpcb_ptr z, mpcb_srcptr z1, mpcb_srcptr z2)
Function: void mpcb_sqr (mpcb_ptr z, mpcb_srcptr z1)
Function: void mpcb_pow_ui (mpcb_ptr z, mpcb_srcptr z1, unsigned long int e)
Function: void mpcb_sqrt (mpcb_ptr z, mpcb_srcptr z1)
Function: void mpcb_div (mpcb_ptr z, mpcb_srcptr z1, mpcb_srcptr z2)
Function: void mpcb_div_2ui (mpcb_ptr z, mpcb_srcptr z1, unsigned long int e)

These are the exact counterparts of the corresponding functions mpc_neg, mpc_add and so on, but on complex balls instead of complex numbers.

Function: int mpcb_can_round (mpcb_srcptr z, mpfr_prec_t prec_re, mpfr_prec_t prec_im, mpc_rnd_t rnd)

If the function returns true (a non-zero number), then rounding any of the complex numbers in the ball to a complex number with precision prec_re of its real and precision prec_im of its imaginary part and rounding mode rnd yields the same result and rounding direction value, cf. return-value. If the function returns false (that is, 0), then it could not conclude, or there are two numbers in the ball which would be rounded to a different complex number or in a different direction. Notice that the function works in a best effort mode and errs on the side of caution by potentially returning false on a roundable ball; this is consistent with computational functions not necessarily returning the smallest enclosing ball.

If z contains the result of evaluating some mathematical function through a sequence of calls to mpcb functions, starting with exact complex numbers, that is, balls of radius 0, then a return value of true indicates that rounding any value in the ball (its centre is readily available) in direction rnd yields the correct result of the function and the correct rounding direction value with the usual MPC semantics.

Notice that when the precision of z is larger than prec_re or prec_im, the centre need not be representable at the desired precision, and in fact the ball need not contain a representable number at all to be “roundable”. Even worse, when rnd is a directed rounding mode for the real or the imaginary part and the ball of non-zero radius contains a representable number, the return value is necessarily false. Even worse, when the rounding mode for one part is to nearest, the corresponding part of the centre of the ball is representable and the ball has a non-zero radius, then the return value is also necessarily false, since even if rounding may be possible, the rounding direction value cannot be determined.

Function: int mpcb_round (mpc_ptr c, mpcb_srcptr z, mpc_rnd_t rnd)

Set c to the centre of z, rounded in direction rnd, and return the corresponding rounding direction value. If mpcb_can_round, called with z, the precisions of c and the rounding mode rnd returns true, then this function does what is expected, it “correctly rounds the ball” and returns a rounding direction value that is valid for all of the ball. As explained above, the result is then not necessarily (in the presence of directed rounding with radius different from 0, it is rather necessarily not) an element of the ball.

Next: , Previous: , Up: GNU MPC   [Index]