libinum v1.0.1

This is the official documentation for libinum, version 1.0.1.

Table of contents

Description

libinum is a C library that implements several algorithms, intended for applications in numerical and symbolic computations.

libinum defines functions to perform the following tasks:

In addition, libinum defines the following structures:

As a general rule, the functions and structures can be used with several different data types. For instance, the root finding and integration functions can manipulate double precision, extended precision, or multi-precision floating point numbers. See Data types for details.

Data types

One of the guiding principles of libinum is that it should should not push a specific data type on users. If, hypothetically, the computation being carried out requires the precision of the floating point numbers to be large, then multi-precision floating point numbers might be preferable to double precision numbers. If, instead, computation time is more important than precision, then extended precision floats may be preferred. The structures and functions defined in libinum, therefore, support various data types.

We will now briefly discuss the data types that may be used to represent numbers.

Integers

Fixed-precision integers

The standard integer types of the C language are char, short int, int, long int and long long int, as well as their unsigned counterparts. The number of bits, and, in consequence, the available range of integers, of each of these types, is platform-dependent. In cases where this could be problematic, fixed-precision integers can be used: with 8 bits: int8_t, 16: int16_t, 32: int32_t and 64: int64_t.

libinum provides a function, print_datatype_info(), that prints the equivalence table between char, short int, int, long int and long long int and int8_t, int16_t, int32_t and int64_t for the implementation of the C library used to compile libinum.

Multi-precision integers

In cases where the required number of bits of integers exceeds 64, one can use multi-precision integers, which can be as large as will fit in memory. However, computation times can greatly increase when using multi-precision integers. The implementation used in libinum is provided by the GNU GMP library.

Floating point numbers

Let us first start with a description of floating point numbers. A real number \(x\) is approximated, with a precision of \(m\) bits, as a collection of three numbers: the sign \(s\in\{-1,1\}\); mantissa, whose binary expansion is denoted by \(a_1\cdots a_m\); and exponent \(e\in\mathbb Z\). The approximate value of the number \(x\) is obtained from its sign, mantissa and exponent by $$ x=s\times a_1.a_2\cdots a_m\times2^e. $$ The precision of a floating point number is the number of bits allocated to its sign and mantissa (since \(a_1\) is necessarily equal to \(1\), it is not stored, so the precision of \(x\) is \(m\) instead of \(m+1\)).

In libinum, floating point number may either be double, extended or multi-precision numbers.

Double precision floating point numbers

Double precision floats are represented using the double type.

The precision and maximal and minimal values of the exponent of double floats depends on the compiler. Their values can be printed using the libinum function print_datatype_info().

For example, using version 5.3.0 of the GNU GCC compiler on the x86-64 architecture, the precision is of 53 bits (i.e. 15 decimal digits), and the maximal and minimal values of the exponent are 1024 and -1021.

Extended precision floating point numbers

Extended precision floats are represented using the long double type. They require more memory than double precision floats, and, as a consequence, slightly more computation time.

The precision and maximal and minimal values of the exponent of long double floats depends on the compiler. Their values can be printed using the libinum function print_datatype_info().

For example, using version 5.3.0 of the GNU GCC compiler on the x86-64 architecture, the precision is of 64 bits (i.e. 18 decimal digits), and the maximal and minimal values of the exponent are 16384 and -16381.

Multi-precision floating point numbers

A multi-precision floating point number is a floating point number whose precision can be set to an arbitrary value (until the number fills up the entire memory of the computer). In libinum, multi-precision floats are implemented using the GNU MPFR library, and will be called MPFR floats.

The precision of MPFR floats can be set using the function mpfr_set_default_prec (see the MPFR library documentation for details). The default value of the precision is 53 bits (as of version 3.1.4 of the MPFR library).

In addition, the maximal value (in absolute value) of the exponent can be controlled by setting emax using the function mpfr_set_emax (see the MPFR library documentation for details).

Depending on the MPFR implementation, the precision and emax can either be an int or a long int, so its maximal and minimal values are platform dependent. The libinum function print_datatype_info() prints which it is.

Algorithms

What follows is a detailed description of the functions that implement the algorithms supported by libinum.

Root finding

libinum can compute the root of smooth real functions using the Newton-Raphson algorithm. The algorithm converges quickly, provided an adequate approximation of the root in provided, and that the root is not located at a minimum of the function.

Description

Given a function \(f\) and a real number \(x_0\), the algorithm produces a sequence \((x_n)\) which, provided the algorithm converges, tends to a root of \(f\). The sequence is defined as $$ x_{n+1}=x_n-\frac{f(x_{n})}{f'(x_{n})} $$ where \(f'\) denotes the derivative of \(f\). The following estimate holds: $$ |x_{n+1}-x_n|\leqslant \frac12|x_n-x_{n-1}|^2\sup_{\xi\in[x_{n+1},x_n]}\frac{f''(\xi)}{f'(\xi)} $$ so that, provided \(f\) is smooth and its derivative does not vanish in the intervals \([x_{n+1},x_n]\), the algorithm converges quadratically.

Implementation

This algorithm has been implemented using double, extended and multi-precision floats.

Integrals

libinum can compute definite integrals of smooth real functions using Gauss-Legendre quadratures.

Description

The main idea of Gauss-Legendre quadratures is to find a set of abcissa \(\{t_1,\cdots,t_N\}\in[-1,1]^N\) and weights \(\{w_1,\cdots,w_N\}\in\mathbb R^N\) such that, if \(f\) were a polynomial of degree \(\leqslant2N-1\), then the integral would be equal to a discrete sum: $$ \int_{a}^bdx\ f(x)=\sum_{i=1}^Nw_if\left(\frac{a+b}2+t_i\frac{b-a}2\right). $$ In the general case, if \(f\) is \(\mathcal C^{2N}\), then $$ \left|\int_{a}^bdx\ f(x)-\sum_{i=1}^Nw_if\left(\frac{a+b}2+t_i\frac{b-a}2\right)\right| \leqslant \frac{N!^4}{(2N+1)(2N)!^3}\sup_{x\in[a,b]}\frac{d^{2N}f(x)}{dx^{2N}}. $$ The number \(N\) is called the order of the integration.

As it turns out, the abcissa are the roots of the \(N\)-th Legendre polynomial \(L_N\), defined by the recursive equation $$ (N+1)L_{N+1}(x)=(2N+1)xL_N(x)-NL_{N-1}(x),\quad L_0(x)=1,\quad L_1(x)=x $$ and the weights are $$ w_i=\frac2{(1-x_i^2)L_N'(t_i)}. $$

Implementation

This algorithm has been implemented using double, extended and multi-precision floats

Types

What follows is a detailed description of the C types defined in libinum.

array

Array that can be dynamically resized.

Structure

Arrays of several types of objects are defined. They are structures, with the following keys: where the keys correspond to array_*'s must be initialized before they are used, and freed when they are no longer needed. khen the array is freed, the objects it contains are freed as well.

Functions

Extra functions for arrays of ordered objects

In addition, we define the following ordering for the possible objects in an array: If the ordering is defined, then the following functions are defined:

Extra functions for array_char

In addition, the following functions are defined for array_char.

polynomial

A polynomial with real coefficients. A polynomial is represented as an array of coefficients and an array of exponents. For example \(1+x+2x^2\) is represented as ({1.,1.,2.},{0,1,2}).

Structure

The coefficients of polynomials can be double, extended, or multi-precision floats, and their exponents are unsigned integers. Polynomials are represented as structures, with the following keys: where the keys correspond to polynomial_*'s must be initialized before they are used, and freed when they are no longer needed.

Functions

polynomialMV

Multi-variable polynomials. A polynomial is represented as an array of coefficients and an array of arrays of indices, each of which represents a variable. For example, \(3x_1^2+x_1x_2\) is represented as ({3,1},{{1,1},{1,2}}).

Structure

The coefficients of multi-variable polynomials can be integers or multiprecision integers, and the indices are integers. Multi-variable polynomials are represented as structures, with the following keys: where the keys correspond to polynomialMV_*'s must be initialized before they are used, and freed when they are no longer needed. When the polynomial is freed, its coefficients and factors are freed as well.

Functions

Functions

Functions are not, strictly speaking, a C type. Instead, in this section, we describe how functions whose pointers are to be passed to the various algorithms of libinum, should be formatted.

Functions must take 3 arguments, and return an integer (a return code): int f(TYPE* out, TYPE in, void* extra_args) where TYPE is the appropriate data type, e.g. double or mpfr_t. in is the argument of the function, out is made to point to \(f(\)in\()\), and extra_args is a pointer to void, which may contain extra arguments to be passed to the function, for example a parameter that the function depends on. For example, to implement the function \(f_\alpha(x)=x^\alpha\) using MPFR floats, 

int f(mpfr_t out, mpfr_t in, void* extra_args){
  int alpha=*((*int)extra_args);
  mpfr_pow_ui(out, in, alpha, MPFR_RNDN);
  return(0);
}

Utilities

In addition, the following functions are defined

Examples

Examples of programs using libinum are provided with the source code, in the doc/libinum-examples directory. If libinum is installed on the filesystem, then the examples can also be found at /usr/share/doc/libinum/libinum-examples.

Authors

libinum was written by Ian Jauslin.