Multiple Precision GNU Scientific Library


Next: , Up: (dir)

Multiple Precision GNU Scientific Library

This document describes version 0.1a3 of MPGSL, an attempt to rewrite some of the GSL modules using the MPFR primitives. At present this package is an experiment to see “what happens”, that is to verify if there is value in such a port.

The package is distributed under the terms of the GNU General Public License (GPL) and can be downloaded from:

http://github.com/marcomaggi/mpgsl/tree/master

Copyright © 2009 by Marco Maggi marcomaggi@gna.org
Copyright © 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 The GSL Team.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.3 or any later version published by the Free Software Foundation; with Invariant Sections being “GNU Free Documentation License” and “GNU General Public License”, no Front–Cover Texts, and no Back–Cover Texts. A copy of the license is included in the section entitled “GNU Free Documentation License”.

Appendices

Indexes


Next: , Previous: Top, Up: Top

1 Overview of the package

MPGSL is an attempt to rewrite some of the GSL modules using the MPFR primitives. At present this package is an experiment to see “what happens”, that is to verify if there is value in such a port.

The GNU Scientific Library (GSL) is a collection of routines for numerical computing, it has been written from scratch in the C language and it presents a modern API for C programmers, allowing wrappers to be written for high level languages.

The Multiple Precision Floating–Point Reliable library (MPFR) implements arithmetics on floating-point numbers; it is based on the GNU MP library.

MPGSL is developed and tested under the GNU+Linux system and officially it supports only the GNU infrastructure. It requires the GNU C library and compiles fine with the GNU C compiler (-std=c99 language). See the README.build file for requirements of the building infrastructure.


Next: , Up: overview

1.1 Using the library

In the following we assume that the library is installed under the /usr/local hierarchy, the version identifier is 1.0.0 and the interface version is 1.2.

A shell script is installed on the system to provide informations about version numbers and installation directories; it should be /usr/local/bin/mpgsl-config. Run mpgsl-config without arguments for the usage screen. An example of output for the important options is:

     $ mpgsl-config --package-version
     1.0.0
     
     $ mpgsl-config --library-interface-version
     1.2
     
     $ mpgsl-config --pkgincludedir
     /usr/local/include/mpgsl/1.0.0
     
     $ mpgsl-config --libdir
     /usr/local/lib
     
     $ mpgsl-config --library-id
     mpgsl1.2

There is a single header file that must be included in C sources: mpgsl.h; in our source code we can put:

     #include <mpgsl.h>

and use the output of mpgsl-config --pkgincludedir as argument to the -I option of the C preprocessor.

The libraries are installed under the directory obtained by running mpgsl-config --libdir; we can use its output as argument of the -L option to the linker.

If we want to link with the shared or static library we take the output of mpgsl-config --library-id and use it as value for the -l option to the linker.


Next: , Previous: overview using, Up: overview

1.2 Using GNU Autoconf to load MPGSL

If we use GNU Autoconf to configure our application, we can embed in our distribution the file mpgsl.m4, installed under $prefix/share/aclocal, and load it in our configure.ac by putting:

     m4_include(mpgsl.m4)

in aclocal.m4. This macro file defines a function that can be invoked with:

     MPGSL_LIBRARY(1,2)

it finds and invokes mpgsl-config to acquire installation directories and library names. The macro accepts two arguments: the major and minor interface version numbers that are needed; both of the arguments are optional, but it is recommended to specify them.

MPGSL_LIBRARY defines the following symbols:

@MPGSL_INCLUDEDIR@
The directory under which header files are installed, example:
          /usr/local/include/mpgsl/1.0.0

@MPGSL_LIBDIR@
The directory under which shared and static libraries are installed, example:
          /usr/local/lib

@MPGSL_CPPFLAGS@
The preprocessor option -I with the value of @MPGSL_INCLUDEDIR@ appended, example:
          -I/usr/local/include/mpgsl/1.0.0

@MPGSL_LDFLAGS@
The linker option -L with the value of @MPGSL_LIBDIR@ appended, example:
          -L/usr/local/lib

@MPGSL_LIBRARY_ID@
The identifier of the shared or static library, example:
          mpgsl1.2

@MPGSL_LIBS@
The linker option -l with the value of @MPGSL_LIBRARY_ID@ appended, example:
          -lmpgsl1.2


Next: , Previous: overview autoconf, Up: overview

1.3 Using pkgconfig to load MPGSL

pkg-config is a program (a compiled executable) that inspects a database of package meta informations and prints to its standard output channel informations about installed packages and libraries. The MPGSL installs a meta data file for use with pkg-config, it should be:

     ${libdir}/pkgconfig/mpgsl.pc

For the full list of available informations look in the file itself, and remember that the value of all the variables set in the meta file can be printed with:

     $ pkg-config mpgsl --variable=<VARNAME>

an example of output for the important variables is:

     $ pkg-config mpgsl --variable=PACKAGE_VERSION
     1.0.0
     
     $ pkg-config mpgsl --variable=library_interface_version
     1.2
     
     $ pkg-config mpgsl --variable=pkgincludedir
     /usr/local/include/mpgsl/1.0.0
     
     $ pkg-config mpgsl --variable=libdir
     /usr/local/lib
     
     $ pkg-config mpgsl --variable=library_id
     mpgsl1.2

The --libs option outputs something like:

     $ pkg-config mpgsl --libs
     -L/usr/local/lib -lmpgsl1.2

while the list of preprocessor options is:

     $ pkg-config mpgsl --cflags
     -I/usr/local/include/mpgsl/1.0.0

Notice that:

     $ pkg-config mpgsl --modversion
     1.2

because what matters is the version of the library interface, not the package distribution version.


Next: , Previous: overview pkgconfig, Up: overview

1.4 Error handling

Errors are handled with the GSL built in system. In particular the GSL_ERROR macro is used to invoke the selected error handler and report errors to the caller. Error Handling.


Previous: overview errors, Up: overview

1.5 Memory allocation

All memory allocation performed by MPGSL is handled with a custom allocator. A default allocator is preselected; it makes use of the standard malloc(), realloc() and free() functions. Custom memory allocation.

The MPFR library relies on the GMP memory allocation functions; GMP makes use of a customisable allocator whose protocol dictates that no memory allocation errors are to be reported to the caller: in case of error the allocation function must not return.

The MPGSL library takes the same approach: in case of error the custom allocator must not return to the caller.


Next: , Previous: overview, Up: Top

2 Custom memory allocation

MPGSL allocates all the memory through a customisable allocator. By default: the library uses an allocator that simply reverts to the standard malloc(), realloc() and free() functions.


Next: , Up: memory

2.1 Data types

The following API has two main purposes:

  1. To implement the classic alloc, realloc and free operations.
  2. To define an allocator data structure that is small enough to be used directly as function argument (that is: passed by value).
— Data Type: mpgsl_memory_allocator_t

Type of the memory allocator. We may build structures of this type to select a memory allocation policy. Fields description follows.

void * data
Pointer to an allocator data structure holding the state of the allocator; it is used as first argument to the allocator function, it can be NULL.
mpgsl_memory_alloc_fun_t * alloc
Pointer to a function used to allocate/reallocate/free a memory block.

— Function Prototype: void mpgsl_memory_alloc_fun_t (void * data, void * pp, size_t dim)

Allocate, reallocate or free memory blocks.

data
The value of the data field in the allocator data structure.
pp
Internally cast to void **, it is a pointer to the variable that holds the pointer to handle.
dim
The dimension of the memory block. It is expressed in bytes for the MPGSL functions.

The protocol:

This function must work like malloc(), realloc() and free() with the fundamental difference that: in case of error it must not return. What it does instead of returning is not the business of MPGSL; the default allocator will terminate the process with an invocation to exit() with code EXIT_FAILURE.

Example of new memory block allocation:

     mpgsl_memory_allocator_t allocator;
     void * p = NULL;
     
     allocator.alloc(allocator.data, &p, 4096);

example of memory block reallocation:

     mpgsl_memory_allocator_t allocator;
     void * p = ...;
     
     allocator.alloc(allocator.data, &p, 4096);

example of memory release:

     mpgsl_memory_allocator_t allocator;
     void * p = ...;
     
     allocator.alloc(allocator.data, &p, 0);


Previous: memory typedefs, Up: memory

2.2 Public interface

The state of the MPGSL library includes two memory allocators: the default one, which is immutable; the current one, which we can customise.

— Function: mpgsl_memory_allocator_t mpgsl_memory_allocator (void)

Return the current memory allocator.

— Function: mpgsl_memory_allocator_t mpgsl_memory_default_allocator (void)

Return the default memory allocator.

— Function: mpgsl_memory_allocator_t mpgsl_memory_set_allocator (mpgsl_memory_allocator_t new)

Select new as current allocator and return the old one.

— Function: void * mpgsl_malloc (size_t dim)
— Function: void * mpgsl_realloc (void * p, size_t dim)
— Function: void mpgsl_free (void * p)

These functions implement the classic alloc, realloc and free operations using the allocator returned by mpgsl_memory_allocator().

— Function: void * mpgsl_malloc_a (mpgsl_memory_allocator_t allocator, size_t dim)
— Function: void * mpgsl_realloc_a (mpgsl_memory_allocator_t allocator, void * p, size_t dim)
— Function: void mpgsl_free_a (mpgsl_memory_allocator_t allocator, void * p)

These functions implement the classic alloc, realloc and free operations using allocator.

— Function: void mpgsl_memory_alloc (void * data, void * q, size_t dim)

Allocate, reallocate or free a block of memory using the standard calloc(), realloc() and free() functions. data is ignored: it is perfectly correct to invoke this function with data set to NULL.

This function can be used as alloc function in a mpgsl_memory_allocator_t, and it is used in the default MPGSL allocator.

The implementation is:

          void
          mpgsl_memory_alloc (void * dummy, void * q, size_t dim)
          {
            void **  pp = q;
            void *   p;
          
          
            if (0 == dim)
              {
                if (NULL != *pp)
                  {
                    free(*pp);
                    *pp = NULL;
                  }
              }
            else
              {
                p = (NULL == *pp)?
                      calloc(1, dim) : realloc(*pp, dim);
                if (NULL == p)
                  {
                    perror(strerror(errno));
                    exit(EXIT_FAILURE);
                  }
                *pp = p;
              }
          }


Next: , Previous: memory, Up: Top

3 Comparing numbers

Special functions are needed to compare MPFR numbers; this chapter describes a number of comparison algorithms. While there is no doubt in determining when a number is greater or exactly equal to another, when performing rounded computations a criterion is needed to determine if two numbers are equal according to a selected tolerance.


Next: , Up: comparison

3.1 Floating–point specific comparison

The code of mpgsl_fcmp() is a rewriting of the original gsl_fcmp() from the GSL. The original code is by Gert Van den Eynde, and the algorithm is by D. E. Knuth “Seminumerical Algorithms”, Section 4.2.2 (3rd Edition). Approximate Comparison of Floating Point Numbers.

— Function: int mpgsl_fcmp (mpfr_t a, mpfr_t b, mpfr_t epsilon)

Compare a and b and return a standard ternary result to indicate if the two numbers are equal according to the selected tolerance epsilon.

The return value is: 0 if the numbers are equal; -1 if a is less than b; +1 if a is greater than b.

The algorithm is implemented with the following pseudo–code:

          abs_max = max(abs(a), abs(b))
          if (0 == abs_max)
            return 0
          else
            delta      = epsilon * 2^get_exponent(abs_max)
            difference = a - b
          
            if (delta < difference)
              return +1
            else if (difference < -delta)
              return -1
            else
              return 0

that is: the numbers are equal if difference is inside the interval [-delta, delta].

delta and difference are MPFR numbers initialised to have the same precision of epsilon.

delta is computed with rounding set to GMP_RNDD.

difference is computed with rounding set to GMP_RNDU.


Previous: comparison fcmp, Up: comparison

3.2 Equality criteria

The following functions define the absdiff–equality and reldiff–equality criteria.

— Function: int mpgsl_absdiff_equal_p (mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr epsilon)

Return non–zero if a and b are equal according to the following criterion:

          |a - b| < |epsilon|

else return zero.

— Function: int mpgsl_reldiff_equal_p (mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr epsilon)

Return non–zero if a and b are equal according to the following criterion:

          |(a - b)/a| < |epsilon|

else return zero.


Next: , Previous: comparison, Up: Top

4 Iterating over values

The MPGSL implements a number of ways to iterate over the results of an operation. The API described here allows us to unify the iterations, with some limitations. As an example, here is how we iterator over the elements of a vector view:

     mpgsl_vector_view           view;
     mpgsl_vector_view_iterator  iter;
     mpfr_ptr                    D;
     
     for (mpgsl_vector_view_iterator_init(&iter, &view);
          mpgsl_iterator_more(&iter);
          mpgsl_iterator_next(&iter))
       {
         D = mpgsl_iterator_ptr(&iter);
       }

and an operation on vector views of the same size can be performed with:

     mpgsl_vector_view           Rview[1], Aview[1], Bview[1];
     mpgsl_vector_view_iterator  Riter[1], Aiter[1], Biter[1];
     
     for (mpgsl_vector_view_iterator_init(Riter, Rview),
          mpgsl_vector_view_iterator_init(Aiter, Aview),
          mpgsl_vector_view_iterator_init(Biter, Bview);
     
          mpgsl_iterator_more(Riter);
     
          mpgsl_iterator_next(Riter),
          mpgsl_iterator_next(Aiter),
          mpgsl_iterator_next(Biter))
       {
         mpfr_add(mpgsl_iterator_ptr(Riter),
                  mpgsl_iterator_ptr(Aiter),
                  mpgsl_iterator_ptr(Biter),
                  GMP_RNDN);
       }

The limitation of this API is that it does not allow the “more” and “next” operations to report an error to the caller. This is not a problem when iterating over vector or matrix views, because they cannot fail.

The iterator initialisation functions are specific for each iterator.

— Data Type: mpgsl_iterator

Type of iterator structure. It must be the first field in the iterator specific structure.

— Function: int mpgsl_iterator_more (void * iterator)

Return non–zero if the iterator has at least one more element.

— Function: void mpgsl_iterator_next (void * iterator)

Advance the iterator to the next element.

— Function: void * mpgsl_iterator_ptr (void * iterator)


Next: , Previous: iterators, Up: Top

5 Vectors and matrices

The API described in this chapter provides a vector and matrix interface to ordinary C arrays. The memory management of these arrays is implemented using a single underlying type, known as a block. By writing our functions in terms of vectors and matrices we can pass a single structure containing both data and dimensions as an argument without needing additional function parameters.


Next: , Up: compound

5.1 Conventions

In the documentation of block, vector and matrix operations the following argument names are used:

dim
A value of type size_t specifying the number of elements in a block or the length of a vector.
b
A pointer to a mpgsl_block structure.
v
A pointer to a mpgsl_vector_view or mpgsl_vector structure.
m
A pointer to a mpgsl_matrix_view or mpgsl_matrix structure.
a
b
c
Pointers to vector or matrix views to be used as operands of an operation.
r
Pointer to a vector or matrix view to be used as result of an operation.
p
A value of type mp_prec_t, the precision to be used to initialise MPFR numbers.
n
A value of type mp_rnd_t, the rounding mode to be used to compoute a result with the MPFR library.
i
j
Values of type size_t to be used as indexes in blocks, vectors and matrices.
k
Usually a value of type long int to be used as parameter in some special functions.


Next: , Previous: compound conv, Up: compound

5.2 Operations on compound values

The following prototype definitions represent possible arguments for functions that iterate over blocks, vectors and matrices. As example, here is a possible implementation of a mapping function for vector views, and its usage to implement the trigonometric sine function:

     int
     vector_map (mpgsl_r_rn_fun_t *          doit,
                 mpgsl_vector_view *         R,
                 const mpgsl_vector_view *   A,
                 mp_rnd_t                    rnd)
     {
       if (R->size != A->size)
         GSL_ERROR("vector views must have same length",
                   GSL_EBADLEN);
       else
         {
           mpgsl_vector_view_iterator  Riter[1], Aiter[1];
           for (mpgsl_vector_view_iterator_init(Riter, R),
                mpgsl_vector_view_iterator_init(Aiter, A);
                mpgsl_iterator_more(Riter);
                mpgsl_iterator_next(Riter),
                mpgsl_iterator_next(Aiter))
             doit(Riter[0].iterator, Aiter[0].iterator, rnd);
           return GSL_SUCCESS;
         }
     }
     
     int
     vector_sin (mpgsl_vector_view *         R,
                 const mpgsl_vector_view *   A,
                 mp_rnd_t                    rnd)
     {
       return vector_map(mpfr_sin, r, a, rnd);
     }

Each prototype name (example mpgsl_r_rrn_fun_t) has:

All the functions with the following types must return a GSL error code.

— Function Signature: int mpgsl_r_r_fun_t (mpfr_ptr r, mpfr_srcptr a)
— Function Signature: int mpgsl_r_rn_fun_t (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t n)
— Function Signature: int mpgsl_r_ro_fun_t (mpfr_ptr r, mpfr_srcptr a, void * custom_context)
— Function Signature: int mpgsl_r_rno_fun_t (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t n, void * custom_context)

Accept a as operand, perform an operation with or without rounding mode n and store the result in r. r and a are allowed to reference the same number.

These are for functions like mpfr_floor() and mpfr_sin().

— Function Signature: int mpgsl_r_rr_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b)
— Function Signature: int mpgsl_r_rrn_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mp_rnd_t n)
— Function Signature: int mpgsl_r_rro_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, void * custom_context)
— Function Signature: int mpgsl_r_rrno_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mp_rnd_t n, void * custom_context)

Accept a and b as operands, perform an operation with or without rounding mode n and store the result in r. r, a and b are allowed to reference the same number.

The rounding mode version is for functions like mpfr_add() and mpfr_mul().

— Function Signature: int mpgsl_r_rrr_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c)
— Function Signature: int mpgsl_r_rrrn_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t n)
— Function Signature: int mpgsl_r_rrro_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, void * custom_context)
— Function Signature: int mpgsl_r_rrrno_fun_t (mpfr_ptr r, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, mp_rnd_t n, void * custom_context)

Accept a, b and c as operands, perform an operation with or without rounding mode n and store the result in r. r, a, b and c are allowed to reference the same number.

This is for functions like mpfr_fma().


Next: , Previous: compound operations, Up: compound

5.3 Blocks

Vectors and matrices are made by slicing an underlying block which references an array of MPFR numbers. A slice is a set of elements formed from an initial offset and a combination of indices and step–sizes (stride).


Next: , Up: compound block

5.3.1 Data types

In the following descriptions we have to remember that:

— Data Type: mpgsl_block

Data type for arrays of MPFR numbers. It represents the “data part” of both vectors and matrices, and it is responsible of its memory ownership.

MPGSL implements no functions to allocate and release mpgsl_block structures, because this structure is normally embedded in other structures or directly allocated on the stack; but it provides mpgsl_block_alloc() and mpgsl_block_free() (along with a lot of others) to allocate and release the array of numbers.

If the data structures are correctly managed: only one mpgsl_block structure references an allocated array of MPFR numbers.

The structure is defined as follows:

          typedef struct {
            mpfr_ptr      data;
            size_t        size;
          } mpgsl_block;

where data is a pointer to the allocated memory block and size is the number of mpfr_t slots in the array. It is fine to set both the fields to zero to represent an empty block; a block structure is interpreted as empty if the data field is set to NULL.

We can initialise a block from a raw array with:

          #define SIZE            1024
          mpfr_t      array[SIZE];
          mpfr_block  block = {
            .data  = &(array[0]),
            .size  = SIZE
          };


Next: , Previous: compound block types, Up: compound block

5.3.2 Iterations

We can iterate over all the elements in a block with the following code, which, as example, sets all the elements to 3:

     mpgsl_block     block
     mpgsl_block_alloc(&block, 12);
     
     mpfr_ptr  D     = block.data;
     mpfr_ptr  Dover = D + block.size;
     
     for (; D < Dover; ++D)
       {
         mpfr_set_si(D, 3, GMP_RNDN);
       }

or with a compact format:

     mpgsl_block     block
     mpfr_ptr        D, Dover;
     
     mpgsl_block_alloc(&block, 12);
     
     for (Dover=((D=block.data))+block.size; D<Dover; ++D)
       {
         mpfr_set_si(D, 3, GMP_RNDN);
       }

We can also use the following generic iterator API:

     mpgsl_block           block[1];
     mpgsl_block_iterator  iterator[1];
     
     mpgsl_block_alloc(block, 12);
     
     for (mpgsl_block_iterator_init(iterator, block);
          mpgsl_iterator_more(iterator);
          mpgsl_iterator_next(iterator))
       {
         mpfr_set_si(mpgsl_iterator_ptr(iterator), 3, GMP_RNDN);
       }
— Data Type: mpgsl_block_iterator

Type of structure for block iterators.

— Function: void mpgsl_block_iterator_init (mpgsl_block_iterator * iter, mpgsl_block * b)

Initialise an element–by–element block iterator.


Next: , Previous: compound block iter, Up: compound block

5.3.3 Block allocation

The functions for allocating blocks make use of the custom allocator: allocation errors are not reported to the caller. Custom memory allocation. See the MPFR documentation for details on initializing MPFR numbers. Initialization Functions.

All the following functions fill the data field of an already allocated mpgsl_block structure, with a pointer to an newly allocated memory block; it is our responsibility to initialise the size field with the number of slots. To allocate a block of 100 MPFR numbers we do:

     mpgsl_block  block = { .size = 100 };
     
     mpgsl_block_alloc(&block);
— Function: void mpgsl_block_free (mpgsl_block * b)

Finalise the MPFR numbers with mpfr_clear(), then free the array memory with mpgsl_free(), finally set the fields of the block to zero. The block structure itself is not released.

If b references an empty block (data pointer set to NULL) nothing happens.

— Function: void mpgsl_block_alloc (mpgsl_block * b)
— Function: void mpgsl_block_alloc_prec (mpgsl_block * b, mp_prec_t p)

Allocate memory for an array of MPFR numbers, reading the requested size from the size field of the block structure. Store a pointer to the allocated memory in the data field of the block.

The MPFR numbers are initialised but left unset, so their value is NaN. The first function applies mpfr_init() to the numbers, while the second applies mpfr_init2() using the given precision argument.

— Function: void mpgsl_block_alloc_zero (mpgsl_block * b)
— Function: void mpgsl_block_alloc_prec_zero (mpgsl_block * b, mp_prec_t p)

Like the above functions but initialize all the elements to zero.

— Function: void mpgsl_block_alloc_value (mpgsl_block * b, mpfr_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_ui (mpgsl_block * b, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_si (mpgsl_block * b, long int x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_uj (mpgsl_block * b, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_sj (mpgsl_block * b, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_d (mpgsl_block * b, double x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_ld (mpgsl_block * b, long double x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_z (mpgsl_block * b, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_q (mpgsl_block * b, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_f (mpgsl_block * b, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_value_decimal64 (mpgsl_block * b, _Decimal64 x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value (mpgsl_block * b, mp_prec_t p, mpfr_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_ui (mpgsl_block * b, mp_prec_t p, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_si (mpgsl_block * b, mp_prec_t p, long int x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_uj (mpgsl_block * b, mp_prec_t p, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_sj (mpgsl_block * b, mp_prec_t p, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_d (mpgsl_block * b, mp_prec_t p, double x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_ld (mpgsl_block * b, mp_prec_t p, long double x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_z (mpgsl_block * b, mp_prec_t p, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_q (mpgsl_block * b, mp_prec_t p, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_f (mpgsl_block * b, mp_prec_t p, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_block_alloc_prec_value_decimal64 (mpgsl_block * b, mp_prec_t p, _Decimal64 x, mp_rnd_t r)

Like the above functions but set all the numbers to x.

The functions with x of type _Decimal64 are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: void mpgsl_block_pick_array_ui (mpgsl_block * b, unsigned long int * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_si (mpgsl_block * b, long int * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_uj (mpgsl_block * b, uintmax_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_sj (mpgsl_block * b, intmax_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_d (mpgsl_block * b, double * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_ld (mpgsl_block * b, long double * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_z (mpgsl_block * b, mpz_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_q (mpgsl_block * b, mpq_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_f (mpgsl_block * b, mpf_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_array_decimal64 (mpgsl_block * b, _Decimal64 * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_ui (mpgsl_block * b, mp_prec_t p, unsigned long int * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_si (mpgsl_block * b, mp_prec_t p, long int * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_uj (mpgsl_block * b, mp_prec_t p, uintmax_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_sj (mpgsl_block * b, mp_prec_t p, intmax_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_d (mpgsl_block * b, mp_prec_t p, double * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_ld (mpgsl_block * b, mp_prec_t p, long double * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_z (mpgsl_block * b, mp_prec_t p, mpz_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_q (mpgsl_block * b, mp_prec_t p, mpq_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_f (mpgsl_block * b, mp_prec_t p, mpf_t * array, mp_rnd_t r)
— Function: void mpgsl_block_pick_prec_array_decimal64 (mpgsl_block * b, mp_prec_t p, _Decimal64 * array, mp_rnd_t r)

Allocate a new array of MPFR numbers and initialise it with values from the array referenced by array; the size of the block is read from the size field of the structure referenced by b, and a pointer to the allocated array is stored in the data field.

The functions with array of type _Decimal64 * are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

Usage example:

          #define DIM            5
          
          const double    array[DIM] = {
            1.1, 2.2, 3.3, 4.4, 5.5
          };
          mpgsl_block     block = { .size = DIM };
          
          mpgsl_block_pick_array_d(&block, array, GMP_RNDN);
          {
            ...
          }
          mpgsl_block_free(&block);


Next: , Previous: compound block alloc, Up: compound block

5.3.4 Initialisation and finalisation

See the MPFR documentation for details on initializing MPFR numbers. Initialization Functions.

— Function: void mpgsl_block_clear (mpgsl_block * b)

Apply mpfr_clear() to all the numbers in b.

— Function: void mpgsl_block_init (mpgsl_block * b)
— Function: void mpgsl_block_init_prec (mpgsl_block * b, mp_prec_t p)

Apply mpfr_init() or mpfr_init2() to all the numbers in b.

— Function: void mpgsl_block_init_zero (mpgsl_block * b)
— Function: void mpgsl_block_init_prec_zero (mpgsl_block * b, mp_prec_t p)

Apply mpfr_init() or mpfr_init2() to all the numbers in b and set them to zero.

— Function: void mpgsl_block_init_value (mpgsl_block * b, mpfr_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_ui (mpgsl_block * b, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_si (mpgsl_block * b, long int x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_uj (mpgsl_block * b, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_sj (mpgsl_block * b, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_d (mpgsl_block * b, double x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_ld (mpgsl_block * b, long double x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_z (mpgsl_block * b, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_q (mpgsl_block * b, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_f (mpgsl_block * b, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_value_decimal64 (mpgsl_block * b, _Decimal64 x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value (mpgsl_block * b, mp_prec_t p, mpfr_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_ui (mpgsl_block * b, mp_prec_t p, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_si (mpgsl_block * b, mp_prec_t p, long int x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_uj (mpgsl_block * b, mp_prec_t p, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_sj (mpgsl_block * b, mp_prec_t p, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_d (mpgsl_block * b, mp_prec_t p, double x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_ld (mpgsl_block * b, mp_prec_t p, long double x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_z (mpgsl_block * b, mp_prec_t p, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_q (mpgsl_block * b, mp_prec_t p, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_f (mpgsl_block * b, mp_prec_t p, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_block_init_prec_value_decimal64 (mpgsl_block * b, mp_prec_t p, _Decimal64 x, mp_rnd_t r)

Apply mpfr_init() or mpfr_init2() to all the numbers in b and set them to x.

The functions with x of type _Decimal64 are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.


Next: , Previous: compound block init, Up: compound block

5.3.5 Accessing blocks

The following accessors are defined mainly to make it easier to write bindings to MPGSL for high level languages, which may use a foreign functions interface.

— Function: void mpgsl_block_set_size (mpgsl_block * b, size_t size)

Store size in the size field of the block.

— Function: size_t mpgsl_block_get_size (mpgsl_block * b)

Return the size field of the structure referenced by b.

— Function: mpfr_ptr mpgsl_block_data (mpgsl_block * b)

Return the data field of the structure referenced by b.

— Function: int mpgsl_block_set (mpgsl_block * b, size_t i, mpfr_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_ui (mpgsl_block * b, size_t i, unsigned long int x, mp_rnd_t r)
— Function: int mpgsl_block_set_si (mpgsl_block * b, size_t i, long int x, mp_rnd_t r)
— Function: int mpgsl_block_set_uj (mpgsl_block * b, size_t i, uintmax_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_sj (mpgsl_block * b, size_t i, intmax_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_d (mpgsl_block * b, size_t i, double x, mp_rnd_t r)
— Function: int mpgsl_block_set_ld (mpgsl_block * b, size_t i, long double x, mp_rnd_t r)
— Function: int mpgsl_block_set_z (mpgsl_block * b, size_t i, mpz_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_q (mpgsl_block * b, size_t i, mpq_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_f (mpgsl_block * b, size_t i, mpf_t x, mp_rnd_t r)
— Function: int mpgsl_block_set_decimal64 (mpgsl_block * b, size_t i, _Decimal64 x, mp_rnd_t r)

Set to x the MPFR number at index i in the array. Return a GSL error code.

The function with x of type _Decimal64 is available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: int mpgsl_block_get (mpgsl_block * b, size_t i, mpfr_ptr * p)

Store in the variable referenced by p a pointer to the MPFR number at index i in the array. Return a GSL error code.


Next: , Previous: compound block access, Up: compound block

5.3.6 Predicates

— Function: int mpgsl_block_same_dimension_p (const mpgsl_block * a, const mpgsl_block * b)

Return true if the two referenced blocks have the same dimensions.


Next: , Previous: compound block pred, Up: compound block

5.3.7 Comparison

The following functions compare the MPFR numbers in the blocks. Comparing numbers, for details on the equality criteria.

Exact and approximate equality
— Function: int mpgsl_block_equal_p (const mpgsl_block * a, const mpgsl_block * b)

Return non–zero if the two blocks have equal values according to mpfr_cmp(); else return zero.

— Function: int mpgsl_block_fcmp_equal_p (const mpgsl_block * a, const mpgsl_block * b, mpfr_t epsilon)

Return non–zero if the two blocks have equal values according to mpgsl_fcmp(), using epsilon tolerance; else return zero.

— Function: int mpgsl_block_absdiff_equal_p (const mpgsl_block * a, const mpgsl_block * b, mpfr_t epsilon)

Return non–zero if the two blocks have numbers satisfying the absdiff–equality criterion; else return zero.

— Function: int mpgsl_block_reldiff_equal_p (const mpgsl_block * a, const mpgsl_block * b, mpfr_t epsilon)

Return non–zero if the two blocks have numbers satisfying the reldiff–equality criterion; else return zero.

Other predicates
— Function: int mpgsl_block_zero_p (const mpgsl_block * a)
— Function: int mpgsl_block_absdiff_zero_p (const mpgsl_block * a, mpfr_t epsilon)

Return non–zero if all the elements of a are zero, respectively, according to mpfr_zero_p() and the absdiff–equality criterion.

The fcmp–equality and reldiff–equality criteria are not suitable to test for approximate zero values.


Previous: compound block cmp, Up: compound block

5.3.8 Input and output

The following functions are simple wrappers for mpfr_fprintf(); they are meant moslty to produce debugging messages. Formatted Output Functions. No fancy formatting is supported by this API.

The return values are the numbers of characters written to the stream, or a negative value if an error occurred.

— Function: int mpgsl_block_fprintf (FILE * f, const char * s, const mpgsl_block * b)

Print the elements from b to f using s for each number. Usage example:

          mpgsl_block     block;
          
          mpgsl_block_fprintf(stderr, "%Rf", &block);
— Function: int mpgsl_block_printf (const char * s, const mpgsl_block * b)

Print the elements from b to stdout using s.

— Function: int mpgsl_block_print (const mpgsl_block * b)

Print the elements from b to stdout using %Rf as format.

— Function: int mpgsl_block_print_stderr (const mpgsl_block * b)

Print the elements from b to stderr using %Rf as format.


Next: , Previous: compound block, Up: compound

5.4 Vectors

This chapter describes the vector API of MPGSL; vectors are one–dimensional views over arrays of MPFR numbers. Vectors are included in MPGSL only as a convenience, because for algebra operations we better use matrices.

All the functions described in this section do bounds checking on vector indexes.


Next: , Up: compound vector

5.4.1 Data types

In the following descriptions we have to remember that:

— Data Type: mpgsl_vector_view

Data type of one–dimensional views over arrays of MPFR numbers. It describes how to visit an array of mpfr_t numbers, but it is not responsible of its memory ownership.

It is defined as follows:

          typedef struct {
            mpfr_ptr      data;
            int           stride;
            size_t        size;
          } mpgsl_vector_view;

where:

data
Is a pointer to the first element in the view.
stride
The offset (for mpfr_ptr pointers) between two adjacent elements in the view. If mpfr_ptr p references an element in the view: p + stride is the next element towards the logical end, p - stride is the next element towards the logical beginning. stride can be positive or negative, but not zero (unless the view is empty).
size
Is the number of elements in the view. If size is zero the view is empty.

— Data Type: mpgsl_vector

Data type of vectors of MPFR numbers. It is defined as follows:

          typedef struct {
            mpgsl_vector_view  view;
            mpgsl_block        block;
          } mpgsl_vector;

structures of this type are allocated with mpgsl_vector_alloc() and released with mpgsl_vector_free(), and:

Specialised allocators exist for mpgsl_vector values that not owning the data block.

Notice that a pointer to the vector is a pointer to the view; this is exploited by the functions acting upon vectors which accept void * arguments to represent mpgsl_vector_view * values.


Next: , Previous: compound vector types, Up: compound vector

5.4.2 Iterations

In the following example we iterate over the elements of an array through a vector view of size 9, stride 3 starting at offset 6:

     mpfr_t  array[1024];
     mpgsl_vector_view view = {
       .data   = &(array[6]),
       .stride = 3,
       .size   = 9
     }
     
     mpfr_ptr   D     = view.data;
     mpfr_ptr   Dover = D + view.size * view.stride;
     int increasing = (0 < view.stride);
     
     for (; increasing? D<Dover : Dover<D; D += view.stride)
       {
         ...
       }

We can also use the following generic iterator API:

     mpgsl_vector_view *         view;
     mpgsl_vector_view_iterator  iterator[1];
     
     view = ...;
     
     for (mpgsl_vector_view_iterator_init(iterator, view);
          mpgsl_iterator_more(iterator);
          mpgsl_iterator_next(iterator))
       {
         mpfr_set_si(mpgsl_iterator_ptr(iterator), 3, GMP_RNDN);
       }
— Data Type: mpgsl_vector_view_iterator

Type of structure for vector view iterators.

— Function: void mpgsl_vector_view_iterator_init (mpgsl_vector_view_iterator * iter, void * view)

Initialise an element–by–element vector view iterator. view is interpreted as pointer to a mpgsl_vector_view structure.


Next: , Previous: compound vector iter, Up: compound vector

5.4.3 Vector allocation

The vector allocator functions make use of the custom allocator: allocation errors are not reported to the caller, so the returned pointer is always non–NULL. Custom memory allocator. See the MPFR documentation for details on initializing MPFR numbers. Initialization Functions.

— Function: void mpgsl_vector_free (mpgsl_vector * v)

Free a previously allocated vector structure.

If v->block.data is non–NULL: the MPFR numbers in the underlying array are cleared with mpfr_clear() and the array itself is released.

— Function: mpgsl_vector * mpgsl_vector_alloc_empty (void)

Allocate a mpgsl_vector structure with all the fields initialised to zero. It represents an empty vector.

— Function: mpgsl_vector * mpgsl_vector_alloc (void)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec (size_t dim, mp_prec_t p)

Allocate a new structure for a vector of length dim.

A new array of MPFR numbers is allocated for the elements of the vector, and stored in the block field of the vector structure. The block is “owned” by the vector, and it will be deallocated when the vector is deallocated.

All the numbers are initialised but left unset, so their value is NaN. The first function makes use of mpfr_init(); the second function makes use of mpfr_init2() and sets to p the precision of all the numbers.

— Function: mpgsl_vector * mpgsl_vector_alloc_zero (size_t dim)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_zero (size_t dim, mp_prec_t p)

Like the above functions but set all the numbers to zero.

— Function: mpgsl_vector * mpgsl_vector_alloc_value (size_t dim, mpfr_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_ui (size_t dim, unsigned long int x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_si (size_t dim, long int x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_uj (size_t dim, uintmax_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_sj (size_t dim, intmax_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_d (size_t dim, double x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_ld (size_t dim, long double x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_z (size_t dim, mpz_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_q (size_t dim, mpq_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_f (size_t dim, mpf_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_value_decimal64 (size_t dim, _Decimal64 x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value (size_t dim, mp_prec_t p, mpfr_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_ui (size_t dim, mp_prec_t p, unsigned long int x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_si (size_t dim, mp_prec_t p, long int x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_uj (size_t dim, mp_prec_t p, uintmax_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_sj (size_t dim, mp_prec_t p, intmax_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_d (size_t dim, mp_prec_t p, double x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_ld (size_t dim, mp_prec_t p, long double x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_z (size_t dim, mp_prec_t p, mpz_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_q (size_t dim, mp_prec_t p, mpq_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_f (size_t dim, mp_prec_t p, mpf_t x, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_alloc_prec_value_decimal64 (size_t dim, mp_prec_t p, _Decimal64 x, mp_rnd_t r)

Like the above functions but set all the numbers to x.

The functions with x of type _Decimal64 are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: mpgsl_vector * mpgsl_vector_pick_array_ui (size_t dim, const unsigned long int * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_si (size_t dim, const long int * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_uj (size_t dim, const uintmax_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_sj (size_t dim, const intmax_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_d (size_t dim, const double * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_ld (size_t dim, const long double * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_z (size_t dim, const mpz_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_q (size_t dim, const mpq_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_f (size_t dim, const mpf_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_array_decimal64 (size_t dim, const _Decimal64 * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_ui (size_t dim, const mp_prec_t p, unsigned long int * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_si (size_t dim, const mp_prec_t p, long int * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_uj (size_t dim, const mp_prec_t p, uintmax_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_sj (size_t dim, const mp_prec_t p, intmax_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_d (size_t dim, const mp_prec_t p, double * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_ld (size_t dim, const mp_prec_t p, long double * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_z (size_t dim, const mp_prec_t p, mpz_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_q (size_t dim, const mp_prec_t p, mpq_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_f (size_t dim, const mp_prec_t p, mpf_t * array, mp_rnd_t r)
— Function: mpgsl_vector * mpgsl_vector_pick_prec_array_decimal64 (size_t dim, const mp_prec_t p, _Decimal64 * array, mp_rnd_t r)

Allocate a new vector and initialise it with values from the array referenced by array, of size dim.

The functions with array of type _Decimal64 * are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

Usage example:

          #define SIZE            5
          
          const double    array[5] = {
            1.1, 2.2, 3.3, 4.4, 5.5
          };
          mpgsl_vector *  vector;
          
          vector = mpgsl_vector_pick_array_d(SIZE, array, GMP_RNDN);
          {
            ...
          }
          mpgsl_vector_free(vector);


Next: , Previous: compound vector alloc, Up: compound vector

5.4.4 Vector views


Next: , Up: compound vector view
5.4.4.1 Vector views over blocks and arrays

The following functions can be used to visit an array of MPFR numbers, provided that we allocate a mpgsl_block to describe it. Usage example:

     #define LEN    100
     
     mpgsl_vector_view  view;
     int                e;
     mpfr_t             array[LEN];
     mpgsl_block        block = {
       .data = &(array[0]),
       .size = LEN
     };
     
     e = mpgsl_vector_init_view_from_block(&view, &block,
                                           3, 9, 3);
     if (GSL_SUCCESS != e)
       {
         /* handle the error */
       }
— Function: int mpgsl_vector_init_view_from_block (mpgsl_vector_view * dst, mpgsl_block * src, size_t offset, size_t dim, size_t stride)

Initialise an already allocated vector view referenced by dst, with data from the raw block referenced by src. Return a GSL error code.

The view will start at offset from the beginning of the block, will have dim elements and will visit the block with stride.

— Function: int mpgsl_vector_alloc_view_from_block (mpgsl_vector_view ** dstp, mpgsl_block * src, size_t offset, size_t dim, size_t stride)

Wrapper for mpgsl_vector_init_view_from_block() that allocates a new vector view structure using mpgsl_malloc(). A pointer to the new structure is stored in the variable referenced by dstp.


Next: , Previous: compound vector view block, Up: compound vector view
5.4.4.2 Vector views over vector views

The following functions can be used to visit a vector view using a vector view with different offset and/or stride. Usage example:

     #define LEN     100
     
     mpgsl_vector_view  one, two;
     int                e;
     mpfr_t             array[LEN];
     mpgsl_block        block = {
       .data = &(array[0]),
       .size = LEN
     };
     
     e = mpgsl_vector_init_view_from_block(&one, &block, 3, 9, 3);
     if (GSL_SUCCESS != e) { /* handle the error */ }
     
     e = mpgsl_vector_init_view_from_view(&two, &one, 1, 3, 3);
     if (GSL_SUCCESS != e) { /* handle the error */ }
— Function: int mpgsl_vector_init_view_from_view (mpgsl_vector_view * dst, void * src, size_t offset, size_t dim, size_t stride)

Initialise an already allocated vector view referenced by dst, with data from the vector view referenced by src. Return a GSL error code.

src is of type void * so that it can reference both a mpgsl_vector and a mpgsl_vector_view structure.

The view will start at offset from the beginning of the source, will have dim elements and will visit the source with stride.

Notice that offset, dim and stride are relative to the source view, not the underlying array of MPFR numbers. So, for example, if the source view has stride 3 over the array and the destination view has stride 4 over the source view: the destination view has stride 4 * 3 = 12 over the array.

— Function: int mpgsl_vector_alloc_view_from_view (mpgsl_vector_view ** dstp, void * src, size_t offset, size_t dim, size_t stride)

Wrapper for mpgsl_vector_init_view_from_view() that allocates a new vector view structure using mpgsl_malloc(). A pointer to the new structure is stored in the variable referenced by dstp.


Previous: compound vector view view, Up: compound vector view
5.4.4.3 Vectors over vector views

The following functions can be used to visit a raw block, or a vector view, using a full vector structure rather than a simple vector view. The vector structure is initialised to not own the underlying block of MPFR numbers.

From the operations point of view: there is nothing that can be done with a full vector structure, that cannot be done with a simple vector view. But there are unfortunate application architectures in which we cannot distinguish a pointer to a vector view from a pointer to a full vector; in these cases, when using dynamically allocated memory, we cannot determine the correct destructor to invoke (mpgsl_free() or mpgsl_vector_free()), so the solution is to use only full vectors so that we know that the destructor is mpgsl_vector_free().

Usage example:

     #define LEN     100
     
     mpgsl_vector  vector;
     int           e;
     mpfr_t        array[LEN];
     mpgsl_block   block = {
       .data = &(array[0]),
       .size = LEN
     };
     
     e = mpgsl_vector_init_vector_from_block(&vector, &block,
                                             3, 9, 3);
     if (GSL_SUCCESS != e) { /* handle the error */ }
— Function: int mpgsl_vector_init_vector_from_block (mpgsl_vector * dst, mpgsl_block * src, size_t offset, size_t dim, size_t stride)

Initialise an already allocated vector structure referenced by dst, with data from the raw block referenced by src. Return a GSL error code.

The view will start at offset from the beginning of the block, will have dim elements and will visit the block with stride.

— Function: int mpgsl_vector_alloc_vector_from_block (mpgsl_vector ** dstp, mpgsl_block * src, size_t offset, size_t dim, size_t stride)

Wrapper for mpgsl_vector_init_vector_from_block() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by dstp.

— Function: int mpgsl_vector_init_vector_from_view (mpgsl_vector * dst, void * src, size_t offset, size_t dim, size_t stride)

Initialise an already allocated vector structure referenced by dst, with data from the vector view referenced by src. Return a GSL error code.

src is of type void * so that it can reference both a mpgsl_vector and a mpgsl_vector_view structure.

The view will start at offset from the beginning of the source, will have dim elements and will visit the source with stride.

Notice that offset, dim and stride are relative to the source view, not the underlying array of MPFR numbers. So, for example, if the source view has stride 3 over the array and the destination view has stride 4 over the source view: the destination view has stride 4 * 3 = 12 over the array.

— Function: int mpgsl_vector_alloc_vector_from_view (mpgsl_vector ** dstp, void * src, size_t offset, size_t dim, size_t stride)

Wrapper for mpgsl_vector_init_vector_from_view() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by dstp.


Next: , Previous: compound vector view, Up: compound vector

5.4.5 Copying vectors

— Function: mpgsl_vector * mpgsl_vector_clone (void * v, mp_rnd_t n)

Interpret v as a pointer to a vector view and build a clone of the represented vector. mpgsl_vector_alloc() is used to allocate a new vector structure and underlying array of MPFR numbers, then the new vector is initialised with values from the view.

The precision of all the numbers in the destination vector matches the one in the source view. n is used as rounding mode for all the set operations.

— Function: int mpgsl_vector_copy (void * dst, const void * src, mp_rnd_t n)

Interpret both dst and src as pointers to vector views and copy the elements from src to dst. The two views must have the same length. n is used as rounding mode for all the set operations. Return a GSL error code.


Next: , Previous: compound vector copy, Up: compound vector

5.4.6 Dimension related functions

All the following functions interpret their void * arguments as pointers to vector views.

— Function: int mpgsl_vector_same_dimension_p (const void * a, const void * b)
— Function: int mpgsl_vector_same_dimension3_p (const void * a, const void * b, const void * c)
— Function: int mpgsl_vector_same_dimension4_p (const void * a, const void * b, const void * c, const void * d)

Return true if the referenced vector views have the same dimensions.


Next: , Previous: compound vector dim, Up: compound vector

5.4.7 Comparison functions

The following functions compare the MPFR numbers in vector views. Comparing numbers, for details on the equality criteria. All the following functions interpret their void * arguments as pointers to vector views.

Exact and approximate equality
— Function: int mpgsl_vector_equal_p (const void * a, const void * b)

Return non–zero if the two vector views have equal values according to mpfr_cmp(); else return zero.

— Function: int mpgsl_vector_fcmp_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two vector views have equal values according to mpgsl_fcmp(), using epsilon tolerance; else return zero.

— Function: int mpgsl_vector_absdiff_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two vector views have numbers satisfying the absdiff–equality criterion; else return zero.

— Function: int mpgsl_vector_reldiff_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two vector views have numbers satisfying the reldiff–equality criterion; else return zero.

Other predicates
— Function: int mpgsl_vector_zero_p (const void * v)
— Function: int mpgsl_vector_absdiff_zero_p (const void * v, mpfr_srcptr epsilon)

Return non–zero if all the elements of v are zero, respectively, according to mpfr_zero_p() and the absdiff–equality criterion.

The fcmp–equality and reldiff–equality criteria are not suitable to test for approximate zero values.

— Function: int mpgsl_vector_ispos_p (const void * v)
— Function: int mpgsl_vector_isneg_p (const void * v)
— Function: int mpgsl_vector_isnonpos_p (const void * v)
— Function: int mpgsl_vector_isnonneg_p (const void * v)

These functions return non–zero if all the elements of the vector v are strictly positive, strictly negative, non–positive or non–negative respectively, and 0 otherwise.


Next: , Previous: compound vector cmp, Up: compound vector

5.4.8 Accessing vector fields

All the v arguments are interpreted as pointers to vector views.

— Function: int mpgsl_vector_get (const void * v, size_t i, mpfr_ptr * dstp)

Store in the variable referenced by dstp a pointer to the i-th element of v. Return a GSL error code.

— Function: int mpgsl_vector_set (void * v, size_t i, mpfr_srcptr x, mp_rnd_t r)
— Function: int mpgsl_vector_set_ui (void * v, size_t i, unsigned long int x, mp_rnd_t r)
— Function: int mpgsl_vector_set_si (void * v, size_t i, long int x, mp_rnd_t r)
— Function: int mpgsl_vector_set_uj (void * v, size_t i, uintmax_t x, mp_rnd_t r)
— Function: int mpgsl_vector_set_sj (void * v, size_t i, intmax_t x, mp_rnd_t r)
— Function: int mpgsl_vector_set_d (void * v, size_t i, double x, mp_rnd_t r)
— Function: int mpgsl_vector_set_ld (void * v, size_t i, long double x, mp_rnd_t r)
— Function: int mpgsl_vector_set_z (void * v, size_t i, mpz_t x, mp_rnd_t r)
— Function: int mpgsl_vector_set_q (void * v, size_t i, mpq_t x, mp_rnd_t r)
— Function: int mpgsl_vector_set_f (void * v, size_t i, mpf_t x, mp_rnd_t r)
— Function: int mpgsl_vector_set_decimal64 (void * v, size_t i, _Decimal64 x, mp_rnd_t r)

Interpret v as a pointer to a vector view structure and set the i-th element to x. Return a GSL error code.

The function with x of type _Decimal64 is available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: void mpgsl_vector_set_all (void * v, mpfr_srcptr x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_ui (void * v, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_si (void * v, long int x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_uj (void * v, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_sj (void * v, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_d (void * v, double x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_ld (void * v, long double x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_z (void * v, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_q (void * v, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_f (void * v, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_vector_set_all_decimal64 (void * v, _Decimal64 x, mp_rnd_t r)

Interpret v as a pointer to a vector view structure and set all the values of v to x.

The function with x of type _Decimal64 is available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: void mpgsl_vector_set_all_zero (void * v)

Interpret v as a pointer to a vector view structure and set all the elements to zero.

— Function: int mpgsl_vector_set_basis (void * v, size_t i)

Interpret v as a pointer to a vector view structure and make it a basis vector by setting all the elements to zero except for the i-th element which is set to one. Return a GSL error code.


Next: , Previous: compound vector access, Up: compound vector

5.4.9 Finding maximum and minimum elements

— Function: int mpgsl_vector_min (void * v, mpfr_ptr * minp)
— Function: int mpgsl_vector_max (void * v, mpfr_ptr * maxp)
— Function: int mpgsl_vector_minmax (void * v, mpfr_ptr * minp, mpfr_ptr * maxp)

Interpret v as pointer to a vector view and store: in the variable referenced by minp a pointer to the minimum value, in the variable referenced by maxp a pointer to the maximum value. Return a GSL error code.

It is an error if the view holds a NaN element.

— Function: int mpgsl_vector_abs_min (void * v, mpfr_ptr * minp)
— Function: int mpgsl_vector_abs_max (void * v, mpfr_ptr * maxp)
— Function: int mpgsl_vector_abs_minmax (void * v, mpfr_ptr * minp, mpfr_ptr * maxp)

Interpret v as pointer to a vector view and store: in the variable referenced by minp a pointer to the number whose absolute value is minimum, in the variable referenced by maxp a pointer to the number whose absolute value is maximum. Return a GSL error code.

It is an error if the view holds a NaN element.

— Function: int mpgsl_vector_min_index (void * v, size_t * minp)
— Function: int mpgsl_vector_max_index (void * v, size_t * maxp)
— Function: int mpgsl_vector_minmax_index (void * v, size_t * minp, size_t * maxp)

Interpret v as pointer to a vector view and store: in the variable referenced by minp the index of the minimum value, in the variable referenced by maxp the index of the maximum value. Return a GSL error code.

It is an error if the view holds a NaN element.

— Function: int mpgsl_vector_abs_min_index (void * v, size_t * minp)
— Function: int mpgsl_vector_abs_max_index (void * v, size_t * maxp)
— Function: int mpgsl_vector_abs_minmax_index (void * v, size_t * minp, size_t * maxp)

Interpret v as pointer to a vector view and store: in the variable referenced by minp the index of the number whose absolute value is minimum, in the variable referenced by maxp the index of the number whose absolute value is maximum. Return a GSL error code.

It is an error if the view holds a NaN element.


Next: , Previous: compound vector minmax, Up: compound vector

5.4.10 Swapping elements and permutations

— Function: int mpgsl_vector_swap (void * v, void * w)

Interpret both the arguments as pointers to vector views and exchange the elements of the vectors. The two vectors must have the same length. Return a GSL error code.

— Function: int mpgsl_vector_swap_elements (void * v, size_t i, size_t j)

Interpret v as pointer to a vector view and exchange the i-th and j-th elements. Return a GSL error code.

— Function: void mpgsl_vector_reverse (void * v)

Interpret v as pointer to a vector view and reverse the order of the elements.

— Function: int mpgsl_vector_permute (void * v, gsl_permutation * perm)

Apply the permutation to the vector view.

— Function: int mpgsl_vector_permute_inverse (void * v, gsl_permutation * perm)

Apply the inverse of the permutation to the vector view. This function is the inverse of mpgsl_vector_permute().


Next: , Previous: compound vector swap, Up: compound vector

5.4.11 Operations on vectors

All the functions described in this section interpret their void* arguments as pointers to vector views and return a GSL error code. The n argument is used as rounding mode for the operation. For all the functions:


Next: , Up: compound vector operations
5.4.11.1 Mapping functions

All the following functions map the function argument f over the elements of the operands to produce a vector result. r is interpreted as pointer to a vector view in which the result must be stored. a, b and c are interpreted as pointers to vector views in which the operands are stored. n is the rounding mode. custom_context is a pointer to a custom value that can be used by f.

All the map function names have prefix mpgsl_map_ and as suffix a signature identifying the arguments; the signature has the following parts:

All the map functions return a GSL error code.

Data types for operations, for the definitions of the function arguments.

— Function: int mpgsl_map_vr__vr (mpgsl_r_r_fun_t * f, void * r, const void * a)
— Function: int mpgsl_map_vr__vr_n (mpgsl_r_rn_fun_t * f, void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_map_vr__vr_o (mpgsl_r_ro_fun_t * f, void * r, const void * a, void * custom_context)
— Function: int mpgsl_map_vr__vr_no (mpgsl_r_rno_fun_t * f, void * r, const void * a, mp_rnd_t n, void * custom_context)

Apply an operation to a single vector operand.

— Function: int mpgsl_map_vr__vr_vr (mpgsl_r_rr_fun_t * f, void * r, const void * a, const void * b)
— Function: int mpgsl_map_vr__vr_vr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * a, const void * b, mp_rnd_t n)
— Function: int mpgsl_map_vr__vr_vr_o (mpgsl_r_rro_fun_t * f, void * r, const void * a, const void * b, void * custom_context)
— Function: int mpgsl_map_vr__vr_vr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * a, const void * b, mp_rnd_t n, void * custom_context)

Apply an operation to two vector operands.

— Function: int mpgsl_map_vr__vr_vr_vr (mpgsl_r_rrr_fun_t * f, void * r, const void * a, const void * b, const void * c)
— Function: int mpgsl_map_vr__vr_vr_vr_n (mpgsl_r_rrrn_fun_t * f, void * r, const void * a, const void * b, const void * c, mp_rnd_t n)
— Function: int mpgsl_map_vr__vr_vr_vr_o (mpgsl_r_rrro_fun_t * f, void * r, const void * a, const void * b, const void * c, void * custom_context)
— Function: int mpgsl_map_vr__vr_vr_vr_no (mpgsl_r_rrrno_fun_t * f, void * r, const void * a, const void * b, const void * c, mp_rnd_t n, void * custom_context)

Apply an operation to three vector operands.

— Function: int mpgsl_map_vr__vr_lr (mpgsl_r_rr_fun_t * f, void * r, const void * a, mpfr_srcptr lambda)
— Function: int mpgsl_map_vr__vr_lr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * a, mpfr_srcptr lambda, mp_rnd_t n)
— Function: int mpgsl_map_vr__vr_lr_o (mpgsl_r_rro_fun_t * f, void * r, const void * a, mpfr_srcptr lambda, void * custom_context)
— Function: int mpgsl_map_vr__vr_lr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * a, mpfr_t lambda, mp_rnd_t n, void * custom_context)

Apply an operation a vector operand and a scalar operand.


Next: , Previous: compound vector map, Up: compound vector operations
5.4.11.2 Operation prototypes

All the following function prototypes can be used to declare pointers to functions for the operations described in this chapter. All the functions are meant to:

— Function Signature: int mpgsl_vr__vr_fun_t (void * r, const void * a)
— Function Signature: int mpgsl_vr__vr_n_fun_t (void * r, const void * a, mp_rnd_t n)

Perform an operation on a single vector operand. This is for functions like mpgsl_vector_round() and mpgsl_vector_sin().

— Function Signature: int mpgsl_vr__vr_vr_fun_t (void * r, const void * a, const void * b)
— Function Signature: int mpgsl_vr__vr_vr_n_fun_t (void * r, const void * a, const void * b, mp_rnd_t n)

Perform an operation on two vector operands. This is for functions like mpgsl_vector_add().

— Function Signature: int mpgsl_vr__vr_vr_vr_fun_t (void * r, const void * a, const void * b, const void * c)
— Function Signature: int mpgsl_vr__vr_vr_vr_n_fun_t (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)

Perform an operation on three vector operands. This is for functions like mpgsl_vector_fma().

— Function Signature: int mpgsl_vr__vr_lr_fun_t (void * r, const void * a, mpfr_srcptr lambda)
— Function Signature: int mpgsl_vr__vr_lr_n_fun_t (void * r, const void * a, mpfr_srcptr lambda, mp_rnd_t n)

Perform an operation on a vector operand and a scalar operand. This is for functions like mpgsl_vector_add_scalar().


Next: , Previous: compound vector op types, Up: compound vector operations
5.4.11.3 Arithmetics between vectors

Basic Arithmetic Functions.

— Function: int mpgsl_vector_add (void * r, const void * a, const void * b, mp_rnd_t n)

Add the elements from b to the elements from a.

— Function: int mpgsl_vector_sub (void * r, const void * a, const void * b, mp_rnd_t n)

Subtract the elements from b to the elements from a.

— Function: int mpgsl_vector_mul (void * r, const void * a, const void * b, mp_rnd_t n)

Multiply the elements from b to the elements from a.

— Function: int mpgsl_vector_div (void * r, const void * a, const void * b, mp_rnd_t n)

Divide the elements from a by the elements from b.

— Function: int mpgsl_vector_fma (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)

Fused multiply and add. Multiply a and b element by element, then add the elements from c.

— Function: int mpgsl_vector_fms (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)

Fused multiply and subtract. Multiply a and b element by element, then subtract the elements from c.

— Function: int mpgsl_vector_agm (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the arithmetic–geometric mean between the elements from a and the elements from b.

— Function: int mpgsl_vector_hypot (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the Euclidean norm of the elements from a and the elements from b. This is the square root of the sum between the square of a and the square of b.

— Function: int mpgsl_vector_sqr (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_sqrt (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_rec_sqrt (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_cbrt (void * r, const void * a, mp_rnd_t n)

Compute the square, square root, reciprocal square root and cubic square root.

— Function: int mpgsl_vector_root (void * r, const void * a, unsigned long int k, mp_rnd_t rnd)

Compute the k-th root.

— Function: int mpgsl_vector_pow (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the power of the elements from a raised to the elements of b.

— Function: int mpgsl_vector_neg (void * r, const void * a, mp_rnd_t n)

Compute the opposite.

— Function: int mpgsl_vector_abs (void * r, const void * a, mp_rnd_t n)

Compute the absolute value.


Next: , Previous: compound vector arithmetics vector, Up: compound vector operations
5.4.11.4 Arithmetics between a vector and a scalar
— Function: int mpgsl_vector_add_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Add the scalar x to all the elements from a.

— Function: int mpgsl_vector_sub_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Subtract the scalar x from all the elements from a.

— Function: int mpgsl_vector_mul_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Multiply the scalar x to all the elements from a.

— Function: int mpgsl_vector_div_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Divide all the elements from a by the scalar x.

— Function: int mpgsl_vector_pow_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Compute the power of the elements from a raised to the scalar x.

— Function: int mpgsl_vector_swapdiv_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Divide the scalar x by all the elements from a.

— Function: int mpgsl_vector_swappow_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Compute the power of x raised to the elements from a.


Next: , Previous: compound vector arithmetics scalar, Up: compound vector operations
5.4.11.5 Logarithms and exponentials

Special Functions.

— Function: int mpgsl_vector_log (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_log2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_log10 (void * r, const void * a, mp_rnd_t n)

Compute the natural logarithm, the base 2 logarithm and the base 10 logarithm.

— Function: int mpgsl_vector_log1p (void * r, const void * a, mp_rnd_t n)

Compute the sum between 1 and the elements from a, then apply the natural logarithm to the result.

— Function: int mpgsl_vector_exp (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_exp2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_exp10 (void * r, const void * a, mp_rnd_t n)

Compute the exponential, the power of 2 and the power of 10.

— Function: int mpgsl_vector_expm1 (void * r, const void * a, mp_rnd_t n)

Compute the exponential of the elements from a, then subtract 1 from the result.


Next: , Previous: compound vector func log, Up: compound vector operations
5.4.11.6 Trigonometric functions

Special Functions.

— Function: int mpgsl_vector_sin (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_cos (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_tan (void * r, const void * a, mp_rnd_t n)

Compute the sine, cosine and tangent functions.

— Function: int mpgsl_vector_sec (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_csc (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_cot (void * r, const void * a, mp_rnd_t n)

Compute the secant, cosecant and cotangent functions.

— Function: int mpgsl_vector_asin (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_acos (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_atan (void * r, const void * a, mp_rnd_t n)

Compute the arc sine, arc cosine and arc tangent functions.

— Function: int mpgsl_vector_atan2 (void * r, const void * y, const void * x, mp_rnd_t n)

Compute the second form of arc tangent of the operands. Special Functions, for details on the expression.


Next: , Previous: compound vector func trigo, Up: compound vector operations
5.4.11.7 Hyperbolic functions

Special Functions.

— Function: int mpgsl_vector_sinh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_cosh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_tanh (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic sine, cosine and tangent functions.

— Function: int mpgsl_vector_sech (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_csch (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_coth (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic secant, cosecant and cotangent functions.

— Function: int mpgsl_vector_asinh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_acosh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_atanh (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic arc sine, arc cosine and arc tangent functions.


Next: , Previous: compound vector func hyper, Up: compound vector operations
5.4.11.8 Special functions

Special Functions.

— Function: int mpgsl_vector_eint (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_li2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_gamma (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_zeta (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_erf (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_erfc (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_j0 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_j1 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_y0 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_y1 (void * r, const void * a, mp_rnd_t n)

Compute miscellaneous special functions implemented by MPFR. Special Functions, for details.

— Function: int mpgsl_vector_jn (void * r, long k, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_yn (void * r, long k, const void * a, mp_rnd_t n)

Compute the first and second kinds of Bessel functions of order k.


Previous: compound vector func special, Up: compound vector operations
5.4.11.9 Integer related functions

Integer Related Functions

— Function: int mpgsl_vector_rint (void * r, const void * a, mp_rnd_t n)

Round to the nearest representable integer in the given rounding mode.

— Function: int mpgsl_vector_rint_ceil (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_rint_floor (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_rint_round (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_vector_rint_trunc (void * r, const void * a, mp_rnd_t n)

Round in various directions. See the original MPFR documentation for details.

— Function: int mpgsl_vector_ceil (void * r, const void * a)
— Function: int mpgsl_vector_floor (void * r, const void * a)
— Function: int mpgsl_vector_round (void * r, const void * a)
— Function: int mpgsl_vector_trunc (void * r, const void * a)

Round in various directions. See the original MPFR documentation for details.

— Function: int mpgsl_vector_frac (void * r, const void * a, mp_rnd_t n)

Extract the fractional part.


Previous: compound vector operations, Up: compound vector

5.4.12 Input and output functions

The following functions are simple wrappers for mpfr_fprintf(); they are meant moslty to produce debugging messages. Formatted Output Functions. No fancy formatting is supported by this API.

The void * arguments are interpreted as pointers to vector views. The return values are the numbers of characters output to the channel.

— Function: int mpgsl_vector_fprintf (FILE * file, const char * format, const void * v)

Print the elements from v to file using format for each number. Usage example:

          mpgsl_vector     vector;
          
          mpgsl_vector_fprintf(stderr, "%Rf", &vector);
— Function: int mpgsl_vector_printf (const char * format, const void * v)

Print the elements from v to stdout using format.

— Function: int mpgsl_vector_print (const void * v)

Print the elements from v to stdout using %Rf as format.

— Function: int mpgsl_vector_print_stderr (const void * v)

Print the elements from v to stderr using %Rf as format.


Previous: compound vector, Up: compound

5.5 Matrices

This chapter describes the matrix API of MPGSL; matrices are bidimensional views over arrays of MPFR numbers. MPGSL defines data types for the following objects:

matrix view descriptors
Parameters to constructors of proper matrix views over arrays of MPFR numbers; view descriptors have the purpose to be “easy” to declare by humans.
matrix views
Describe bidimensional views over arrays of MPFR numbers to be used for matrix operations. They have the purpose to provide efficent access to elements in row–major order. They do not own the underlying array memory.
matrices
Full matrices with ownership of underlying array memory.


Next: , Up: compound matrix

5.5.1 Data types

In the following descriptions we have to remember that:

Row and column indexes are zero based.

— Data Type: mpgsl_matrix_view

Data type of bidimensional views over arrays of MPFR numbers. It is not responsible of its memory ownership.

It is defined as follows:

          typedef struct {
            mpfr_ptr   data;
            int        row_stride;
            int        col_step;
            int        crlf;
            size_t     num_of_rows;
            size_t     num_of_cols;
          } mpgsl_matrix_view;

where:

data
Is a pointer to the first element in the view, that is element (0, 0) in the matrix view.
num_of_rows
It is the number of rows in the view and the number of elements in each column.
num_of_cols
It is the number of columns in the view and the number of elements in each row.
row_stride
The offset (for mpfr_ptr pointers) between two adjacent elements on the same row. It can be positive or negative, but not zero.
col_step
It is the offset (for mpfr_ptr pointers) between two adjacent elements in the view on the same column. It can be positive or negative, but not zero.
crlf
A precomputed value defined as:
               crlf = col_step - num_of_cols * row_stride

and stored in the matrix view for efficiency.

— Data Type: mpgsl_matrix

Data type of matrices of MPFR numbers. It is defined as follows:

          typedef struct {
            mpgsl_matrix_view  view;
            mpgsl_block        block;
          } mpgsl_matrix;

structures of this type are allocated with mpgsl_matrix_alloc() and released with mpgsl_matrix_free(), and:

Specialised allocators exist to allocate mpgsl_matrix values not owning the data block.

Notice that a pointer to the matrix is a pointer to the view; this is exploited by the functions acting upon matrices which accept void * arguments to represent mpgsl_matrix_view * values.

— Data Type: mpgsl_matrix_view_over_block

Data type of structures describing a matrix view on top of a block. It is used as parameter for matrix view constructors. It is defined as:

          typedef struct {
            size_t     num_of_rows;
            size_t     num_of_cols;
            int        row_stride;
            int        col_step;
            int        transpose;
            size_t     offset;
          } mpgsl_matrix_view_over_block;

with:

offset
The offset (for mpfr_ptr pointers) from the beginning of the underlying block to element (0, 0) in the view.
num_of_rows
The number of rows in the view.
num_of_cols
The number of columns in the view.
row_stride
The offset (for mpfr_ptr pointers) between two adjacent elements on the same row. It can be positive or negative, but not zero. It has the same value of the homonymous field in mpgsl_matrix_view.
col_step
It is the offset (for mpfr_ptr pointers) between two adjacent elements in the view on the same column. It can be positive or negative, but not zero. It has the same value of the homonymous field in mpgsl_matrix_view.
transpose
If non–zero the described matrix is transposed when building the view structure.

— Data Type: mpgsl_matrix_view_over_view

Data type of structures describing a matrix view (upper) on top of another matrix view (lower). It is used as parameter for matrix view constructors. It is defined as:

          typedef struct {
            size_t     num_of_rows;
            size_t     num_of_cols;
            int        row_stride;
            int        col_stride;
            int        transpose;
            size_t     first_row;
            size_t     first_col;
          } mpgsl_matrix_view_over_view;

with:

first_row
first_col
The row and column indexes in the underlying matrix view of element (0, 0) of the upper view.
num_of_rows
The number of rows in the upper view.
num_of_cols
The number of columns in the upper view.
row_stride
How many columns in the lower view to step between adjacent columns in the upper view.
col_stride
How many rows in the lower view to step between adjacent rows in the upper view.
transpose
If non–zero the described view is transposed when building the real view structure.


Next: , Previous: compound matrix types, Up: compound matrix

5.5.2 Discussion on matrix views

The representation of matrix views is a bit complex so, with the purpose of understanding fields usage, we discuss it here by examples. A matrix view structure is defined as:

     typedef struct {
       mpfr_ptr  data;
       int       row_stride;
       int       col_step;
       int       crlf;
       size_t    num_of_rows;
       size_t    num_of_cols;
     } mpgsl_matrix_view;

the fact that row_stride has stride in it, while col_step has step in it has a purpose; another data type has a col_stride field which is similar to row_stride in the MPGSL modeling of matrices. The fields are meant to allow:

The rules of moving a pointer around always are:


Next: , Up: compound matrix intro view
5.5.2.1 Plain direct views

Let's have an array of 100 MPFR numbers, initialised with values from 0 to 99 in increasing order:

     mpgsl_block  block = { .size = 100 };
     
     mpgsl_block_alloc(&block);
     for (size_t i=0; i<block.size; ++i)
       mpfr_set_ui(&(block.data[i]), i, GMP_RNDN);

we want to interpret it as a matrix with 10 columns and 10 rows, like the following:

      0  1  2  3  4  5  6  7  8  9
     10 11 12 13 14 15 16 17 18 19
     20 21 22 23 24 25 26 27 28 29
     30 31 32 33 34 35 36 37 38 39
     40 41 42 43 44 45 46 47 48 49
     50 51 52 53 54 55 56 57 58 59
     60 61 62 63 64 65 66 67 68 69
     70 71 72 73 74 75 76 77 78 79
     80 81 82 83 84 85 86 87 88 89
     90 91 92 93 94 95 96 97 98 99

for this we set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = block.data,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = DIM,
       .crlf         = 0
     };

this is straightforward (no?); we can explicitly note that:

Now we want a view over the following submatrix:

      0  1  2  3  4  5  6
     10 11 12 13 14 15 16
     20 21 22 23 24 25 26
     30 31 32 33 34 35 36
     40 41 42 43 44 45 46

so we set:

     mpgsl_matrix_view  view = {
       .data         = block.data,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_step     = DIM,
       .crlf         = DIM - 7 * 1
                    /* col_step - num_of_cols * row_stride */
     };

The col_step is the offset to add to a mpfr_ptr pointer to move from an element to the next on the same column: it is unchanged from the previous full view.

The crlf field is a little complicated, think of it like this: if mpfr_ptr p references 26, by adding col_step == DIM we move to number 36, and by subtracting num_of_cols * row_stride we move to a row_stride before the number 30. The good news is that the formula to compute crlf from col_step and row_stride is always the same in all the views.

Now we can understand how to build a view for the following matrix which skips columns:

      2  4  6
     12 14 16
     22 24 26
     32 34 36
     42 44 46

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[2]),
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = DIM,
       .crlf         = DIM - 3 * 2
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

      2  4  6
     22 24 26
     42 44 46

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[2]),
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = DIM * 2
       .crlf         = DIM * 2 - 3 * 2
                    /* col_step - num_of_cols * row_stride */
     };


Next: , Previous: compound matrix intro view plain direct, Up: compound matrix intro view
5.5.2.2 Plain transposed views

Let's have the same array of MPFR numbers defined in Plain direct views. Let's build a view for the following matrix:

     0 10 20 30 40 50 60 70 80 90
     1 11 21 31 41 51 61 71 81 91
     2 12 22 32 42 52 62 72 82 92
     3 13 23 33 43 53 63 73 83 93
     4 14 24 34 44 54 64 74 84 94
     5 15 25 35 45 55 65 75 85 95
     6 16 26 36 46 56 66 76 86 96
     7 17 27 37 47 57 67 77 87 97
     8 18 28 38 48 58 68 78 88 98
     9 19 29 39 49 59 69 79 89 99

which is the transposition of the full matrix used before; we set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = block.data,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = DIM,
       .col_step     = 1,
       .crlf         = 1 - DIM * DIM
                    /* col_step - num_of_cols * row_stride */
                    /* it is negative */
     };

we see that row_stride and col_step in the transposed view are just the swapped values they had in the direct view. For crlf we use the same formula as before.

Now we want a view over the following transposed submatrix:

     0 10 20 30 40 50 60
     1 11 21 31 41 51 61
     2 12 22 32 42 52 62
     3 13 23 33 43 53 63
     4 14 24 34 44 54 64

so we set:

     mpgsl_matrix_view  view = {
       .data         = block.data,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = DIM,
       .col_step     = 1,
       .crlf         = 1 - 7 * DIM
                    /* col_step - num_of_cols * row_stride */
     };

Now we can understand how to build a view for the following matrix which skips rows:

     2 12 22 32 42 52 62
     4 14 24 34 44 54 64
     6 16 26 36 46 56 66

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[2]),
       .num_of_rows  = 3,
       .num_of_cols  = 7,
       .row_stride   = DIM,
       .col_step     = 2,
       .crlf         = 2 - 7 * DIM
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

     2 22 42
     4 24 44
     6 26 46

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[2]),
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2 * DIM,
       .col_step     = 2,
       .crlf         = 2 - 3 * (2 * DIM)
                    /* col_step - num_of_cols * row_stride */
     };


Next: , Previous: compound matrix intro view plain transpose, Up: compound matrix intro view
5.5.2.3 Reversed columns direct views

Let's have the same array of MPFR numbers defined in Plain direct views. Let's build a view for the following matrix:

      9  8  7  6  5  4  3  2  1  0
     19 18 17 16 15 14 13 12 11 10
     29 28 27 26 25 24 23 22 21 20
     39 38 37 36 35 34 33 32 31 30
     49 48 47 46 45 44 43 42 41 40
     59 58 57 56 55 54 53 52 51 50
     69 68 67 66 65 64 63 62 61 60
     79 78 77 76 75 74 73 72 71 70
     89 88 87 86 85 84 83 82 81 80
     99 98 97 96 95 94 93 92 91 90

which is the full direct matrix used before with the columns reversed; we set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = &(block.data[9]),
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -1,
       .col_step     = DIM
       .crlf         = DIM - DIM * (-1)
                    /* col_step - num_of_cols * row_stride */
     };

the only difference from the direct view is that the row stride is negative.

Now we want a view over the following submatrix:

      9  8  7  6  5  4  3
     19 18 17 16 15 14 13
     29 28 27 26 25 24 23
     39 38 37 36 35 34 33
     49 48 47 46 45 44 43

so we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[9]),
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -1,
       .col_step     = DIM
       .crlf         = DIM - 7 * (-1)
                    /* col_step - num_of_cols * row_stride */
     };

Now we can understand how to build a view for the following matrix which skips columns:

      6  4  2
     16 14 12
     26 24 22
     36 34 32
     46 44 42

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[6]),
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = DIM
       .crlf         = DIM - 3 * (-2)
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

      6  4  2
     26 24 22
     46 44 42

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[6]),
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = 2 * DIM
       .crlf         = (2 * DIM) - 3 * (-2)
                    /* col_step - num_of_cols * row_stride */
     };


Next: , Previous: compound matrix intro view revcols direct, Up: compound matrix intro view
5.5.2.4 Reversed columns transposed views

Let's have the same array of MPFR numbers defined in Plain direct views. Let's build a view for the following matrix:

     9 19 29 39 49 59 69 79 89 99
     8 18 28 38 48 58 68 78 88 98
     7 17 27 37 47 57 67 77 87 97
     6 16 26 36 46 56 66 76 86 96
     5 15 25 35 45 55 65 75 85 95
     4 14 24 34 44 54 64 74 84 94
     3 13 23 33 43 53 63 73 83 93
     2 12 22 32 42 52 62 72 82 92
     1 11 21 31 41 51 61 71 81 91
     0 10 20 30 40 50 60 70 80 90

which is the full direct matrix used before with the columns reversed then transposed; notice that the sequence column reversing then transposition, is equivalent to transposition then row reversing. We set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = &(block.data[9]),
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = DIM,
       .col_step     = -1
       .crlf         = -1 - DIM * DIM
                    /* col_step - num_of_cols * row_stride */
     };

Now we want a view over the following submatrix:

     9 19 29 39 49
     8 18 28 38 48
     7 17 27 37 47
     6 16 26 36 46
     5 15 25 35 45
     4 14 24 34 44
     3 13 23 33 43

so we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[9]),
       .num_of_rows  = 7,
       .num_of_cols  = 5,
       .row_stride   = DIM,
       .col_step     = -1,
       .crlf         = -1 - 5 * DIM
                    /* col_step - num_of_cols * row_stride */
     };

Now we can understand how to build a view for the following matrix which skips rows:

     6 16 26 36 46
     4 14 24 34 44
     2 12 22 32 42

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[6]),
       .num_of_rows  = 3,
       .num_of_cols  = 5,
       .row_stride   = DIM,
       .col_step     = -2
       .crlf         = -2 - 5 * DIM
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

     6 26 46
     4 24 44
     2 22 42

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[6]),
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2 * DIM,
       .col_step     = -2,
       .crlf         = -2 - 3 * (2 * DIM)
                    /* col_step - num_of_cols * row_stride */
     };


Next: , Previous: compound matrix intro view revcols trans, Up: compound matrix intro view
5.5.2.5 Reversed rows direct views

Let's have the same array of MPFR numbers defined in Plain direct views. Let's build a view for the following matrix:

     90 91 92 93 94 95 96 97 98 99
     80 81 82 83 84 85 86 87 88 89
     70 71 72 73 74 75 76 77 78 79
     60 61 62 63 64 65 66 67 68 69
     50 51 52 53 54 55 56 57 58 59
     40 41 42 43 44 45 46 47 48 49
     30 31 32 33 34 35 36 37 38 39
     20 21 22 23 24 25 26 27 28 29
     10 11 12 13 14 15 16 17 18 19
      0  1  2  3  4  5  6  7  8  9

which is the full direct matrix used before with the rows reversed; we set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = &(block.data[90]),
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = -DIM,
       .crlf         = -DIM - DIM * 1
                    /* col_step - num_of_cols * row_stride */
     };

Now we want a view over the following submatrix:

     90 91 92 93 94 95 96
     80 81 82 83 84 85 86
     70 71 72 73 74 75 76
     60 61 62 63 64 65 66
     50 51 52 53 54 55 56

so we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[90]),
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_step     = -DIM,
       .crlf         = -DIM - 7 * 1
                    /* col_step - num_of_cols * row_stride */
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 94 96
     82 84 86
     72 74 76
     62 64 66
     52 54 56

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[92]),
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = -DIM,
       .crlf         = -DIM - 3 * 2
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

     92 94 96
     72 74 76
     52 54 56

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[92]),
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = -2 * DIM,
       .crlf         = (-2 * DIM) - 3 * 2
                    /* col_step - num_of_cols * row_stride */
     };


Previous: compound matrix intro view revrows direct, Up: compound matrix intro view
5.5.2.6 Reversed rows transposed views

Let's have the same array of MPFR numbers defined in Plain direct views. Let's build a view for the following matrix:

     90 80 70 60 50 40 30 20 10  0
     91 81 71 61 51 41 31 21 11  1
     92 82 72 62 52 42 32 22 12  2
     93 83 73 63 53 43 33 23 13  3
     94 84 74 64 54 44 34 24 14  4
     95 85 75 65 55 45 35 25 15  5
     96 86 76 66 56 46 36 26 16  6
     97 87 77 67 57 47 37 27 17  7
     98 88 78 68 58 48 38 28 18  8
     99 89 79 69 59 49 39 29 19  9

which is the full direct matrix used before with the rows reversed then transposed; notice that the sequence row reversing then transposition, is equivalent to transposition then column reversing. We set:

     #define DIM     10
     
     mpgsl_matrix_view  view = {
       .data         = &(block.data[90]),
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -DIM,
       .col_step     = 1,
       .crlf         = 1 - DIM * (-DIM)
                    /* col_step - num_of_cols * row_stride */
     };

Now we want a view over the following submatrix:

     90 80 70 60 50 40 30
     91 81 71 61 51 41 31
     92 82 72 62 52 42 32
     93 83 73 63 53 43 33
     94 84 74 64 54 44 34

so we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[90]),
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -DIM,
       .col_step     = 1,
       .crlf         = 1 - 7 * (-DIM)
                    /* col_step - num_of_cols * row_stride */
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 72 52
     93 73 53
     94 74 54
     95 75 55
     96 76 56
     97 77 57
     98 78 58

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[92]),
       .num_of_rows  = 7,
       .num_of_cols  = 3,
       .row_stride   = -2 * DIM,
       .col_step     = 1,
       .crlf         = 1 - 3 * (-2 * DIM)
                    /* col_step - num_of_cols * row_stride */
     };

For the following matrix which skips both rows and columns:

     92 72 52
     94 74 54
     96 76 56
     98 78 58

we set:

     mpgsl_matrix_view  view = {
       .data         = &(block.data[92]),
       .num_of_rows  = 4,
       .num_of_cols  = 3,
       .row_stride   = -2 * DIM,
       .col_step     = 2,
       .crlf         = 2 - 3 * (-2 * DIM)
                    /* col_step - num_of_cols * row_stride */
     };


Next: , Previous: compound matrix intro view, Up: compound matrix

5.5.3 Discussion on matrix view–over–block descriptors

Here we discuss by examples the representation of matrix view–over–block descriptors; they should be a little easier to use than direct matrix views. The view–over–block descriptor structure is defined as:

     typedef struct {
       size_t        num_of_rows;
       size_t        num_of_cols;
       int           row_stride;
       int           col_step;
       int           transpose;
       size_t        offset;
     } mpgsl_matrix_view_over_block;

all the fields but transpose are used to slice the underlying array of MPFR numbers into a rectangular matrix; if transpose is non–zero the matrix is transposed by the view.

View–over–block descriptors are used with the functions mpgsl_matrix_init_view_over_block() and mpgsl_matrix_alloc_view_over_block().


Next: , Up: compound matrix block desc
5.5.3.1 Plain direct view over block descriptors

Let's have an array of 100 MPFR numbers, initialised with values from 0 to 99 in increasing order:

     mpgsl_block  block = { .size = 100 };
     
     mpgsl_block_alloc(&block);
     for (size_t i=0; i<block.size; ++i)
       mpfr_set_ui(&(block.data[i]), i, GMP_RNDN);

we want to interpret it as a matrix with 10 rows and 10 columns, like the following:

      0  1  2  3  4  5  6  7  8  9
     10 11 12 13 14 15 16 17 18 19
     20 21 22 23 24 25 26 27 28 29
     30 31 32 33 34 35 36 37 38 39
     40 41 42 43 44 45 46 47 48 49
     50 51 52 53 54 55 56 57 58 59
     60 61 62 63 64 65 66 67 68 69
     70 71 72 73 74 75 76 77 78 79
     80 81 82 83 84 85 86 87 88 89
     90 91 92 93 94 95 96 97 98 99

for this we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 0,
       .transpose    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = DIM
     };

this is straightforward once we understand that:

offset is the offset of the matrix element (0, 0) from the beginning the block. transpose is zero because the view is a direct slicing of the block.

Now we want a view over the following submatrix:

      0  1  2  3  4  5  6
     10 11 12 13 14 15 16
     20 21 22 23 24 25 26
     30 31 32 33 34 35 36
     40 41 42 43 44 45 46

so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 0,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_step     = DIM
     };

Now we can understand how to build a view for the following matrix which skips columns:

      2  4  6
     12 14 16
     22 24 26
     32 34 36
     42 44 46

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 2,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = DIM
     };

For the following matrix which skips both rows and columns:

      2  4  6
     22 24 26
     42 44 46

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 2,
       .transpose    = 0,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = 2 * DIM
     };


Next: , Previous: compound matrix block desc plain direct, Up: compound matrix block desc
5.5.3.2 Plain transposed view over block descriptors

Let's have the same array of MPFR numbers defined in Plain direct view over block descriptors. Let's build a view for the following matrix:

     0 10 20 30 40 50 60 70 80 90
     1 11 21 31 41 51 61 71 81 91
     2 12 22 32 42 52 62 72 82 92
     3 13 23 33 43 53 63 73 83 93
     4 14 24 34 44 54 64 74 84 94
     5 15 25 35 45 55 65 75 85 95
     6 16 26 36 46 56 66 76 86 96
     7 17 27 37 47 57 67 77 87 97
     8 18 28 38 48 58 68 78 88 98
     9 19 29 39 49 59 69 79 89 99

which is the transposition of the full matrix used before; we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 0,
       .transpose    = 1,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = DIM
     };

that is: we slice the block for a direct view, as before, but set transpose to 1. Again: all the fields but transpose are used to slice the array for a “direct” view, then transpose selects a direct or transposed view.

Now we want a view over the following transposed submatrix:

     0 10 20 30 40 50 60
     1 11 21 31 41 51 61
     2 12 22 32 42 52 62
     3 13 23 33 43 53 63
     4 14 24 34 44 54 64

we have to slice the block for the following direct submatrix:

      0  1  2  3  4
     10 11 12 13 14
     20 21 22 23 24
     30 31 32 33 34
     40 41 42 43 44
     50 51 52 53 54
     60 61 62 63 64

and set transpose to 1; so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 0,
       .transpose    = 1,
       .num_of_rows  = 7,
       .num_of_cols  = 5,
       .row_stride   = 1,
       .col_step     = DIM
     };

Now we can understand how to build a view for the following matrix which skips rows:

     2 12 22 32 42 52 62
     4 14 24 34 44 54 64
     6 16 26 36 46 56 66

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 2,
       .transpose    = 1,
       .num_of_rows  = 7,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = DIM
     };

For the following matrix which skips both rows and columns:

     2 22 42
     4 24 44
     6 26 46

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 2,
       .transpose    = 1,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = 2 * DIM
     };


Next: , Previous: compound matrix block desc plain trans, Up: compound matrix block desc
5.5.3.3 Reversed columns direct view over block descriptors

Let's have the same array of MPFR numbers defined in Plain direct view over block descriptors. Let's build a view for the following matrix:

      9  8  7  6  5  4  3  2  1  0
     19 18 17 16 15 14 13 12 11 10
     29 28 27 26 25 24 23 22 21 20
     39 38 37 36 35 34 33 32 31 30
     49 48 47 46 45 44 43 42 41 40
     59 58 57 56 55 54 53 52 51 50
     69 68 67 66 65 64 63 62 61 60
     79 78 77 76 75 74 73 72 71 70
     89 88 87 86 85 84 83 82 81 80
     99 98 97 96 95 94 93 92 91 90

which is the full direct matrix used before with the columns reversed. To do it we have to apply the rules of usage for the row_stride and col_step fields:

So we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 9,
       .transpose    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -1,
       .col_step     = DIM

Now we want a view over the following submatrix:

      9  8  7  6  5  4  3
     19 18 17 16 15 14 13
     29 28 27 26 25 24 23
     39 38 37 36 35 34 33
     49 48 47 46 45 44 43

so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 9,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -1,
       .col_step     = DIM
     };

Now we can understand how to build a view for the following matrix which skips columns:

      6  4  2
     16 14 12
     26 24 22
     36 34 32
     46 44 42

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 6,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = DIM
     };

For the following matrix which skips both rows and columns:

      6  4  2
     26 24 22
     46 44 42

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 6,
       .transpose    = 0,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = 2 * DIM
     };


Next: , Previous: compound matrix block desc revcols direct, Up: compound matrix block desc
5.5.3.4 Reversed columns transposed view over block descriptors

Let's have the same array of MPFR numbers defined in Plain direct view over block descriptors. Let's build a view for the following matrix:

     9 19 29 39 49 59 69 79 89 99
     8 18 28 38 48 58 68 78 88 98
     7 17 27 37 47 57 67 77 87 97
     6 16 26 36 46 56 66 76 86 96
     5 15 25 35 45 55 65 75 85 95
     4 14 24 34 44 54 64 74 84 94
     3 13 23 33 43 53 63 73 83 93
     2 12 22 32 42 52 62 72 82 92
     1 11 21 31 41 51 61 71 81 91
     0 10 20 30 40 50 60 70 80 90

which is the full direct matrix used before with the columns reversed then transposed; notice that the sequence column reversing then transposition, is equivalent to transposition then row reversing.

We have to use all the fields but transpose to slice the block as a direct matrix with the columns reversed, additionally we set transpose to non–zero. So we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 9,
       .transpose    = 1,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -1,
       .col_step     = DIM
     };

Now we want a view over the following submatrix:

     9 19 29 39 49
     8 18 28 38 48
     7 17 27 37 47
     6 16 26 36 46
     5 15 25 35 45
     4 14 24 34 44
     3 13 23 33 43

so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 9,
       .transpose    = 1,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -1,
       .col_step     = DIM
     };

Now we can understand how to build a view for the following matrix which skips rows:

     6 16 26 36 46
     4 14 24 34 44
     2 12 22 32 42

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 6,
       .transpose    = 1,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = DIM
     };

For the following matrix which skips both rows and columns:

     6 26 46
     4 24 44
     2 22 42

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 6,
       .transpose    = 1,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_step     = 2*DIM
     };


Next: , Previous: compound matrix block desc revcols trans, Up: compound matrix block desc
5.5.3.5 Reversed rows direct view over block descriptors

Let's have the same array of MPFR numbers defined in Plain direct view over block descriptors. Let's build a view for the following matrix:

     90 91 92 93 94 95 96 97 98 99
     80 81 82 83 84 85 86 87 88 89
     70 71 72 73 74 75 76 77 78 79
     60 61 62 63 64 65 66 67 68 69
     50 51 52 53 54 55 56 57 58 59
     40 41 42 43 44 45 46 47 48 49
     30 31 32 33 34 35 36 37 38 39
     20 21 22 23 24 25 26 27 28 29
     10 11 12 13 14 15 16 17 18 19
      0  1  2  3  4  5  6  7  8  9

which is the full direct matrix used before with the rows reversed. To do it we have to apply the rules of usage for the row_stride and col_step fields:

So we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 90,
       .transpose    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = -DIM
     };

Now we want a view over the following submatrix:

     90 91 92 93 94 95 96
     80 81 82 83 84 85 86
     70 71 72 73 74 75 76
     60 61 62 63 64 65 66
     50 51 52 53 54 55 56

so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 90,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_step     = -DIM
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 94 96
     82 84 86
     72 74 76
     62 64 66
     52 54 56

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 92,
       .transpose    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = -DIM
     };

For the following matrix which skips both rows and columns:

     92 94 96
     72 74 76
     52 54 56

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 92,
       .transpose    = 0,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_step     = -2 * DIM
     };


Previous: compound matrix block desc revrows direct, Up: compound matrix block desc
5.5.3.6 Reversed rows transposed view over block descriptors

Let's have the same array of MPFR numbers defined in Plain direct view over block descriptors. Let's build a view for the following matrix:

     90 80 70 60 50 40 30 20 10  0
     91 81 71 61 51 41 31 21 11  1
     92 82 72 62 52 42 32 22 12  2
     93 83 73 63 53 43 33 23 13  3
     94 84 74 64 54 44 34 24 14  4
     95 85 75 65 55 45 35 25 15  5
     96 86 76 66 56 46 36 26 16  6
     97 87 77 67 57 47 37 27 17  7
     98 88 78 68 58 48 38 28 18  8
     99 89 79 69 59 49 39 29 19  9

which is the full direct matrix used before with the rows reversed then transposed; notice that the sequence row reversing then transposition, is equivalent to transposition then column reversing.

We have to use all the fields but transpose to slice the block as a direct matrix with the rows reversed, additionally we set transpose to non–zero. We set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block  desc = {
       .offset       = 90,
       .transpose    = 1,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = -DIM
     };

Now we want a view over the following submatrix:

     90 80 70 60 50 40 30
     91 81 71 61 51 41 31
     92 82 72 62 52 42 32
     93 83 73 63 53 43 33
     94 84 74 64 54 44 34

so we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 90,
       .transpose    = 1,
       .num_of_rows  = 7,
       .num_of_cols  = 5,
       .row_stride   = 1,
       .col_step     = -DIM
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 72 52
     93 73 53
     94 74 54
     95 75 55
     96 76 56
     97 77 57
     98 78 58

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 92,
       .transpose    = 1,
       .num_of_rows  = 3,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_step     = -2 * DIM
     };

For the following matrix which skips both rows and columns:

     92 72 52
     94 74 54
     96 76 56
     98 78 58

we set:

     mpgsl_matrix_view_over_block  desc = {
       .offset       = 92,
       .transpose    = 1,
       .num_of_rows  = 3,
       .num_of_cols  = 4,
       .row_stride   = 2,
       .col_step     = -2 * DIM
     };


Next: , Previous: compound matrix block desc, Up: compound matrix

5.5.4 Discussion on matrix view–over–view descriptors

Here we discuss by examples the representation of matrix views over views. The view–over–view descriptor structure is defined as:

     typedef struct {
       size_t        num_of_rows;
       size_t        num_of_cols;
       int           row_stride;
       int           col_stride;
       int           transpose;
       size_t        first_row;
       size_t        first_col;
     } mpgsl_matrix_view_over_view;

View–over–view descriptors are used with the functions mpgsl_matrix_init_view_over_view() and mpgsl_matrix_alloc_view_over_view().


Next: , Up: compound matrix view desc
5.5.4.1 Plain direct view–over–view descriptors

Let's have an array of 100 MPFR numbers, initialised with values from 0 to 99 in increasing order:

     mpgsl_block  block = { .size = 100 };
     
     mpgsl_block_alloc(&block);
     for (size_t i=0; i<block.size; ++i)
       mpfr_set_ui(&(block.data[i]), i, GMP_RNDN);

we want to interpret it as a matrix with 10 columns and 10 rows, like the following:

      0  1  2  3  4  5  6  7  8  9
     10 11 12 13 14 15 16 17 18 19
     20 21 22 23 24 25 26 27 28 29
     30 31 32 33 34 35 36 37 38 39
     40 41 42 43 44 45 46 47 48 49
     50 51 52 53 54 55 56 57 58 59
     60 61 62 63 64 65 66 67 68 69
     70 71 72 73 74 75 76 77 78 79
     80 81 82 83 84 85 86 87 88 89
     90 91 92 93 94 95 96 97 98 99

for this we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_block over_block = {
       .offset       = 0,
       .transpose    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_step     = DIM
     };

On top of this view we want another view which changes nothing: the upper view defines the same matrix of the lower view. For this we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_stride   = 1,
       .transpose    = 0
     };

we see that:

Now we want a view over the following submatrix of over_block:

      0  1  2  3  4  5  6
     10 11 12 13 14 15 16
     20 21 22 23 24 25 26
     30 31 32 33 34 35 36
     40 41 42 43 44 45 46

so we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_stride   = 1,
       .transpose    = 0
     };

Now we can understand how to build a view for the following matrix which skips columns of over_block:

      2  4  6
     12 14 16
     22 24 26
     32 34 36
     42 44 46

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 2,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = 1,
       .transpose    = 0
     };

For the following matrix which skips both rows and columns of over_block:

      2  4  6
     22 24 26
     42 44 46

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 2,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = 2,
       .transpose    = 0
     };


Next: , Previous: compound matrix view desc plain direct, Up: compound matrix view desc
5.5.4.2 Plain transposed view over view descriptors

Let's have the same array of MPFR numbers defined in Plain direct view–over–view descriptors. We also consider the same direct view over_block defined there.

Let's build a view for the following matrix:

     0 10 20 30 40 50 60 70 80 90
     1 11 21 31 41 51 61 71 81 91
     2 12 22 32 42 52 62 72 82 92
     3 13 23 33 43 53 63 73 83 93
     4 14 24 34 44 54 64 74 84 94
     5 15 25 35 45 55 65 75 85 95
     6 16 26 36 46 56 66 76 86 96
     7 17 27 37 47 57 67 77 87 97
     8 18 28 38 48 58 68 78 88 98
     9 19 29 39 49 59 69 79 89 99

which is the transposition of over_block; we set the fields describing over_block and additionally set transpose to 1:

     #define DIM     10
     
     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_stride   = 1,
       .transpose    = 1
     };

Now we want a view over the following transposed submatrix of over_block:

     0 10 20 30 40 50 60
     1 11 21 31 41 51 61
     2 12 22 32 42 52 62
     3 13 23 33 43 53 63
     4 14 24 34 44 54 64

we slice over_block for the following submatrix:

      0  1  2  3  4  5  6
     10 11 12 13 14 15 16
     20 21 22 23 24 25 26
     30 31 32 33 34 35 36
     40 41 42 43 44 45 46

then set transpose to 1; we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 0,
       .num_of_rows  = 7,
       .num_of_cols  = 5,
       .row_stride   = 1,
       .col_stride   = 1,
       .transpose    = 1
     };

Now we can understand how to build a view for the following matrix which skips rows of the transposed of over_block:

     2 12 22 32 42 52 62
     4 14 24 34 44 54 64
     6 16 26 36 46 56 66

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 2,
       .num_of_rows  = 7,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = 1,
       .transpose    = 1
     };

For the following matrix which skips both rows and columns of the transposed of over_block:

     2 22 42
     4 24 44
     6 26 46

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 2,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = 2,
       .transpose    = 1
     };


Next: , Previous: compound matrix view desc plain trans, Up: compound matrix view desc
5.5.4.3 Reversed columns direct view over view descriptors

Let's have the same array of MPFR numbers defined in Plain direct view–over–view descriptors. We also consider the same direct view over_block defined there.

Let's build a view for the following matrix:

      9  8  7  6  5  4  3  2  1  0
     19 18 17 16 15 14 13 12 11 10
     29 28 27 26 25 24 23 22 21 20
     39 38 37 36 35 34 33 32 31 30
     49 48 47 46 45 44 43 42 41 40
     59 58 57 56 55 54 53 52 51 50
     69 68 67 66 65 64 63 62 61 60
     79 78 77 76 75 74 73 72 71 70
     89 88 87 86 85 84 83 82 81 80
     99 98 97 96 95 94 93 92 91 90

which is the view described by over_block with the columns reversed; we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 9,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -1,
       .col_stride   = 1,
       .transpose    = 0
     }

Now we want a view over the following submatrix:

      9  8  7  6  5  4  3
     19 18 17 16 15 14 13
     29 28 27 26 25 24 23
     39 38 37 36 35 34 33
     49 48 47 46 45 44 43

so we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 9,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -1,
       .col_stride   = 1,
       .transpose    = 0
     };

Now we can understand how to build a view for the following matrix which skips columns:

      6  4  2
     16 14 12
     26 24 22
     36 34 32
     46 44 42

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 6,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_stride   = 1,
       .transpose    = 0
     };

For the following matrix which skips both rows and columns:

      6  4  2
     26 24 22
     46 44 42

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 6,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_stride   = 2,
       .transpose    = 0
     };


Next: , Previous: compound matrix view desc revcols direct, Up: compound matrix view desc
5.5.4.4 Reversed columns transposed view over view descriptors

Let's have the same array of MPFR numbers defined in Plain direct view–over–view descriptors. We also consider the same direct view over_block defined there.

Let's build a view for the following matrix:

     9 19 29 39 49 59 69 79 89 99
     8 18 28 38 48 58 68 78 88 98
     7 17 27 37 47 57 67 77 87 97
     6 16 26 36 46 56 66 76 86 96
     5 15 25 35 45 55 65 75 85 95
     4 14 24 34 44 54 64 74 84 94
     3 13 23 33 43 53 63 73 83 93
     2 12 22 32 42 52 62 72 82 92
     1 11 21 31 41 51 61 71 81 91
     0 10 20 30 40 50 60 70 80 90

which is over_block with the columns reversed then transposed; notice that the sequence column reversing then transposition, is equivalent to transposition then row reversing. We set:

     #define DIM     10
     
     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 9,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = -1,
       .col_stride   = 1,
       .transpose    = 1
     };

Now we want a view over the following submatrix:

     9 19 29 39 49
     8 18 28 38 48
     7 17 27 37 47
     6 16 26 36 46
     5 15 25 35 45
     4 14 24 34 44
     3 13 23 33 43

so we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 9,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = -1,
       .col_stride   = 1,
       .transpose    = 1
     };

Now we can understand how to build a view for the following matrix which skips rows:

     6 16 26 36 46
     4 14 24 34 44
     2 12 22 32 42

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 6,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_stride   = 1,
       .transpose    = 1
     };

For the following matrix which skips both rows and columns:

     6 26 46
     4 24 44
     2 22 42

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 0,
       .first_col    = 6,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = -2,
       .col_stride   = 2,
       .transpose    = 1
     };


Next: , Previous: compound matrix view desc revcols trans, Up: compound matrix view desc
5.5.4.5 Reversed rows direct view over view descriptors

Let's have the same array of MPFR numbers defined in Plain direct view–over–view descriptors. We also consider the same direct view over_block defined there.

Let's build a view for the following matrix:

     90 91 92 93 94 95 96 97 98 99
     80 81 82 83 84 85 86 87 88 89
     70 71 72 73 74 75 76 77 78 79
     60 61 62 63 64 65 66 67 68 69
     50 51 52 53 54 55 56 57 58 59
     40 41 42 43 44 45 46 47 48 49
     30 31 32 33 34 35 36 37 38 39
     20 21 22 23 24 25 26 27 28 29
     10 11 12 13 14 15 16 17 18 19
      0  1  2  3  4  5  6  7  8  9

which is over_block with the rows reversed; we set:

     #define DIM     10
     
     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_stride   = -1,
       .transpose    = 0
     };

Now we want a view over the following submatrix:

     90 91 92 93 94 95 96
     80 81 82 83 84 85 86
     70 71 72 73 74 75 76
     60 61 62 63 64 65 66
     50 51 52 53 54 55 56

so we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 0,
       .num_of_rows  = 5,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_stride   = -1,
       .transpose    = 0
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 94 96
     82 84 86
     72 74 76
     62 64 66
     52 54 56

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 2,
       .num_of_rows  = 5,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = -1,
       .transpose    = 0
     };

For the following matrix which skips both rows and columns:

     92 94 96
     72 74 76
     52 54 56

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 2,
       .num_of_rows  = 3,
       .num_of_cols  = 3,
       .row_stride   = 2,
       .col_stride   = -2,
       .transpose    = 0
     };


Previous: compound matrix view desc revrows direct, Up: compound matrix view desc
5.5.4.6 Reversed rows transposed view over view descriptors

Let's have the same array of MPFR numbers defined in Plain direct view–over–view descriptors. We also consider the same direct view over_block defined there.

Let's build a view for the following matrix:

     90 80 70 60 50 40 30 20 10  0
     91 81 71 61 51 41 31 21 11  1
     92 82 72 62 52 42 32 22 12  2
     93 83 73 63 53 43 33 23 13  3
     94 84 74 64 54 44 34 24 14  4
     95 85 75 65 55 45 35 25 15  5
     96 86 76 66 56 46 36 26 16  6
     97 87 77 67 57 47 37 27 17  7
     98 88 78 68 58 48 38 28 18  8
     99 89 79 69 59 49 39 29 19  9

which is over_block with the rows reversed then transposed; notice that the sequence row reversing then transposition, is equivalent to transposition then column reversing. We set:

     #define DIM     10
     
     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 0,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM,
       .row_stride   = 1,
       .col_stride   = -1,
       .transpose    = 1
     };

Now we want a view over the following submatrix:

     90 80 70 60 50 40 30
     91 81 71 61 51 41 31
     92 82 72 62 52 42 32
     93 83 73 63 53 43 33
     94 84 74 64 54 44 34

so we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 0,
       .num_of_rows  = 7,
       .num_of_cols  = 5,
       .row_stride   = 1,
       .col_stride   = -1,
       .transpose    = 1
     };

Now we can understand how to build a view for the following matrix which skips columns:

     92 72 52
     93 73 53
     94 74 54
     95 75 55
     96 76 56
     97 77 57
     98 78 58

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 2,
       .num_of_rows  = 3,
       .num_of_cols  = 7,
       .row_stride   = 1,
       .col_stride   = -2,
       .transpose    = 1
     };

For the following matrix which skips both rows and columns:

     92 72 52
     94 74 54
     96 76 56
     98 78 58

we set:

     mpgsl_matrix_view_over_view  over_view = {
       .first_row    = 9,
       .first_col    = 2,
       .num_of_rows  = 3,
       .num_of_cols  = 4,
       .row_stride   = 2,
       .col_stride   = -2,
       .transpose    = 1
     };


Next: , Previous: compound matrix view desc, Up: compound matrix

5.5.5 Iterators

In this section we review the methods of iteration over the elements of a matrix. There are two interfaces: direct usage of pointers, generic function–based API.


Next: , Up: compound matrix iter
5.5.5.1 Row iteration

In the following example we iterate over the elements of row 4 (zero based) of a 12 by 12 array of MPFR numbers:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpfr_ptr  D = Dview.data + 4 * Dview.col_step;
     
     for (size_t col=0; col<Dview.num_of_cols; ++col, D+=Dview.stride)
       {
         mpfr_printf("%Rf\n", D);
       }

Alternatively the generic iterator interface can be used:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpgsl_matrix_view_row_iterator   Diter;
     mpfr_ptr  D;
     
     for (mpgsl_matrix_view_row_iterator_init(&Diter, &Dview, 4);
          mpgsl_iterator_more(&Diter);
          mpgsl_iterator_next(&Diter))
       {
         D = mpgsl_iterator_ptr(&Diter);
         mpfr_printf("%Rf\n", D);
       }
— Data Type: mpgsl_matrix_view_row_iterator

Type of data structure representing the state of a row iterator.

— Function: void mpgsl_matrix_view_row_iterator_init (mpgsl_matrix_view_row_iterator * iterator, void * vp, size_t row)

Initialise a row iterator. vp is interpreted as pointer to a matrix view structure.


Next: , Previous: compound matrix iter row, Up: compound matrix iter
5.5.5.2 Column iterations

In the following example we iterate over the elements of column 5 (zero based) of a 12 by 12 array of MPFR numbers:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpfr_ptr  D = Dview.data + 5 * Dview.stride;
     
     for (size_t row=0; row<Dview.num_of_rows; ++row, D+=Dview.col_step)
       {
         mpfr_printf("%Rf\n", D);
       }

Alternatively the generic iterator interface can be used:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpgsl_matrix_view_column_iterator  Diter;
     mpfr_ptr  D;
     
     for (mpgsl_matrix_view_column_iterator_init(&Diter, &Dview, 5);
          mpgsl_iterator_more(&Diter);
          mpgsl_iterator_next(&Diter))
       {
         D = mpgsl_iterator_ptr(&Diter);
         mpfr_printf("%Rf\n", D);
       }
— Data Type: mpgsl_matrix_view_column_iterator

Type of data structure representing the state of a column iterator.

— Function: void mpgsl_matrix_view_column_iterator_init (mpgsl_matrix_view_column_iterator * iterator, void * vp, size_t column)

Initialise a column iterator. vp is interpreted as pointer to a matrix view structure.


Previous: compound matrix iter col, Up: compound matrix iter
5.5.5.3 Element by element iterations

In the following example we iterate of the elements of a 12 by 12 array of MPFR numbers with the classic row–by–row, column–by–column nested loops:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpfr_ptr  D = Dview.data;
     
     for (size_t row=0; row<Dview.num_of_rows; ++row, D+=Dview.crlf)
       for (size_t col=0; col<Dview.num_of_cols; ++col, D+=Dview.stride)
         {
           mpfr_printf("%Rf\n", D);
         }

and to get a single element at coordinates (i, j) equal to (8, 7), we do:

     size_t    i=8, j=7;
     mpfr_ptr  D = Dview.data + i * Dview.col_step + j * Dview.stride;

Alternatively the generic iterator interface can be used:

     #define DIM     12
     mpfr_t             matrix[DIM * DIM];
     mpgsl_matrix_view  Dview = {
       .data         = &(matrix[0]),
       .stride       = 1,
       .crlf         = 0,
       .col_step     = DIM,
       .num_of_rows  = DIM,
       .num_of_cols  = DIM
     };
     
     mpgsl_matrix_view_iterator Diter;
     mpfr_ptr  D;
     
     for (mpgsl_matrix_view_iterator_init(&Diter, &Dview, 4);
          mpgsl_iterator_more(&Diter);
          mpgsl_iterator_next(&Diter))
       {
         D = mpgsl_iterator_ptr(&Diter);
         mpfr_printf("%Rf\n", D);
       }
— Data Type: mpgsl_matrix_view_iterator

Type of data structure representing the state of an element by element iterator.

— Function: void mpgsl_matrix_view_iterator_init (mpgsl_matrix_view_iterator * iterator, void * vp)

Initialise an element by element iterator. vp is interpreted as pointer to a matrix view structure.


Next: , Previous: compound matrix iter, Up: compound matrix

5.5.6 Matrix allocation

The matrix allocator functions make use of the custom allocator: allocation errors are not reported to the caller, so the returned pointer is always non–NULL. Custom memory allocator. See the MPFR documentation for details on initializing MPFR numbers. Initialization Functions.

— Function: void mpgsl_matrix_free (mpgsl_matrix * v)

Free a previously allocated matrix structure.

If v->block.data is non–NULL: the MPFR numbers in the underlying array are cleared with mpfr_clear() and the array itself is released.

— Function: mpgsl_matrix * mpgsl_matrix_alloc_empty (void)

Allocate a mpgsl_matrix structure with all the fields initialised to zero. It represents an empty matrix.

— Function: mpgsl_matrix * mpgsl_matrix_alloc (size_t num_of_rows, size_t num_of_cols)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec (size_t num_of_rows, size_t num_of_cols, mp_prec_t p)

Allocate a new structure for a matrix.

A new array of MPFR numbers is allocated for the elements of the matrix, and stored in the block field of the matrix structure. The block is “owned” by the matrix, and it will be deallocated when the matrix is deallocated.

All the numbers are initialised but left unset, so their value is NaN. The first function makes use of mpfr_init(); the second function makes use of mpfr_init2() and sets to p the precision of all the numbers.

— Function: mpgsl_matrix * mpgsl_matrix_alloc_zero (size_t num_of_rows, size_t num_of_cols)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_zero (size_t num_of_rows, size_t num_of_cols, mp_prec_t p)

Like the above functions but set all the numbers to zero.

— Function: mpgsl_matrix * mpgsl_matrix_alloc_value (size_t num_of_rows, size_t num_of_cols, mpfr_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_ui (size_t num_of_rows, size_t num_of_cols, unsigned long int x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_si (size_t num_of_rows, size_t num_of_cols, long int x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_uj (size_t num_of_rows, size_t num_of_cols, uintmax_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_sj (size_t num_of_rows, size_t num_of_cols, intmax_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_d (size_t num_of_rows, size_t num_of_cols, double x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_ld (size_t num_of_rows, size_t num_of_cols, long double x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_z (size_t num_of_rows, size_t num_of_cols, mpz_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_q (size_t num_of_rows, size_t num_of_cols, mpq_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_f (size_t num_of_rows, size_t num_of_cols, mpf_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_value_decimal64 (size_t num_of_rows, size_t num_of_cols, _Decimal64 x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, mpfr_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_ui (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, unsigned long int x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_si (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, long int x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_uj (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, uintmax_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_sj (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, intmax_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_d (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, double x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_ld (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, long double x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_z (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, mpz_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_q (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, mpq_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_f (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, mpf_t x, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_alloc_prec_value_decimal64 (size_t num_of_rows, size_t num_of_cols, mp_prec_t p, _Decimal64 x, mp_rnd_t r)

Like the above functions but set all the numbers to x.

The functions with x of type _Decimal64 are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: mpgsl_matrix * mpgsl_matrix_pick_array_ui (size_t num_of_rows, size_t num_of_cols, const unsigned long int * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_si (size_t num_of_rows, size_t num_of_cols, const long int * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_uj (size_t num_of_rows, size_t num_of_cols, const uintmax_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_sj (size_t num_of_rows, size_t num_of_cols, const intmax_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_d (size_t num_of_rows, size_t num_of_cols, const double * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_ld (size_t num_of_rows, size_t num_of_cols, const long double * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_z (size_t num_of_rows, size_t num_of_cols, const mpz_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_q (size_t num_of_rows, size_t num_of_cols, const mpq_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_f (size_t num_of_rows, size_t num_of_cols, const mpf_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_array_decimal64 (size_t num_of_rows, size_t num_of_cols, const _Decimal64 * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_ui (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, unsigned long int * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_si (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, long int * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_uj (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, uintmax_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_sj (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, intmax_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_d (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, double * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_ld (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, long double * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_z (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, mpz_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_q (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, mpq_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_f (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, mpf_t * array, mp_rnd_t r)
— Function: mpgsl_matrix * mpgsl_matrix_pick_prec_array_decimal64 (size_t num_of_rows, size_t num_of_cols, const mp_prec_t p, _Decimal64 * array, mp_rnd_t r)

Allocate a new matrix and initialise it with values from the array referenced by array, which must hold at least num_of_rows times num_of_cols elements.

The functions with array of type _Decimal64 * are available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

Usage example:

          #define NUM_OF_ROWS     3
          #define NUM_OF_COLS     4
          
          const long int  array[NUM_OF_ROWS * NUM_OF_COLS] = {
            1,  2,  3,  4,
            5,  6,  7,  8,
            9, 10, 11, 12
          };
          mpgsl_matrix *  matrix;
          
          matrix = mpgsl_matrix_pick_array_d(NUM_OF_ROWS, NUM_OF_COLS,
                                             array, GMP_RNDN);
          {
            ...
          }
          mpgsl_matrix_free(matrix);


Next: , Previous: compound matrix alloc, Up: compound matrix

5.5.7 Allocating and initialising matrix and vector views


Next: , Up: compound matrix view
5.5.7.1 Views over blocks

The following functions can be used to visit an array of MPFR numbers, provided that we allocate a mpgsl_block to describe it.

— Function: int mpgsl_matrix_init_view_over_block (mpgsl_matrix_view * v, mpgsl_block * b, mpgsl_matrix_view_over_block * descr)

Initialise an already allocated matrix view referenced by v, over numbers from the block referenced by b. The view is described by descr. Return a GSL error code.

— Function: int mpgsl_matrix_alloc_view_over_block (mpgsl_matrix_view ** vp, mpgsl_block * b, mpgsl_matrix_view_over_block * descr)

Wrapper for mpgsl_matrix_init_view_over_block() that allocates a new matrix view structure using mpgsl_malloc(). A pointer to the new structure is stored in the variable referenced by vp.


Next: , Previous: compound matrix view block, Up: compound matrix view
5.5.7.2 Views over views

The following functions can be used to visit a matrix view using a matrix view.

— Function: int mpgsl_matrix_init_view_over_view (mpgsl_matrix_view * dst, void * src, mpgsl_matrix_view_over_view * descr)

Initialise an already allocated matrix view referenced by dst, on top of the matrix view referenced by src. The view is described by desc. Return a GSL error code.

src is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_view_over_view (mpgsl_matrix_view ** dstp, void * src, mpgsl_matrix_view_over_view * descr)

Wrapper for mpgsl_matrix_init_view_over_view() that allocates a new matrix view structure using mpgsl_malloc(). A pointer to the new structure is stored in the variable referenced by dstp.


Next: , Previous: compound matrix view view, Up: compound matrix view
5.5.7.3 Full matrix over blocks and vectors

The following functions can be used to visit a raw block, or a matrix view, using a full matrix structure rather than a simple matrix view. The matrix structure is initialised to not own the underlying block of MPFR numbers.

From the operations point of view: there is nothing that can be done with a full matrix structure, that cannot be done with a simple matrix view. But there are unfortunate application architectures in which we cannot distinguish a pointer to a matrix view from a pointer to a full matrix; in these cases, when using dynamically allocated memory, we cannot determine the correct destructor to invoke (mpgsl_free() or mpgsl_matrix_free()), so the solution is to use only full matrices so that we know that the destructor is mpgsl_matrix_free().

— Function: int mpgsl_matrix_init_matrix_over_block (mpgsl_matrix * v, mpgsl_block * b, mpgsl_matrix_view_over_block * descr)

Initialise an already allocated matrix structure referenced by v, over numbers in the block referenced by b. The view is described by descr. Return a GSL error code.

— Function: int mpgsl_matrix_alloc_matrix_over_block (mpgsl_matrix ** vp, mpgsl_block * b, mpgsl_matrix_view_over_block * descr)

Wrapper for mpgsl_matrix_init_matrix_over_block() that allocates a new matrix structure. A pointer to the new structure is stored in the variable referenced by vp.

— Function: int mpgsl_matrix_init_matrix_over_view (mpgsl_matrix * dst, void * src, mpgsl_matrix_view_over_view * descr)

Initialise an already allocated matrix structure referenced by dst, with data from the matrix view referenced by src. The view is described by descr. Return a GSL error code.

src is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_matrix_over_view (mpgsl_matrix ** dstp, void * src, mpgsl_matrix_view_over_view * descr)

Wrapper for mpgsl_matrix_init_matrix_over_view() that allocates a new matrix structure. A pointer to the new structure is stored in the variable referenced by dstp.


Next: , Previous: compound matrix view full, Up: compound matrix view
5.5.7.4 Vector views over matrices
Vector views over rows
— Function: int mpgsl_matrix_init_vector_view_over_row (mpgsl_vector_view * v, void * m, size_t row_index)

Initialise a vector view to reference elements from row row_index of matrix view m. m is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_vector_view_over_row (mpgsl_vector_view ** vp, void * m, size_t row_index)

Wrapper for mpgsl_matrix_init_vector_view_over_row() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by vp.

Vector views over columns
— Function: int mpgsl_matrix_init_vector_view_over_column (mpgsl_vector_view * v, void * m, size_t column_index)

Initialise a vector view to reference elements from column column_index of matrix view m. m is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_vector_view_over_column (mpgsl_vector_view ** vp, void * m, size_t column_index)

Wrapper for mpgsl_matrix_init_vector_view_over_column() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by vp.

Vector views over diagonals

Given the following matrix:

      0  1  2  3  4  5
     10 11 12 13 14 15
     20 21 22 23 24 25
     30 31 32 33 34 35

the diagonal of index 0 is the vector 0 11 22 33; the superdiagonals have positive index and are:

     superdiag index +1      1 12 23 34
     superdiag index +2      2 13 34 25
     superdiag index +3      3 14 25
     superdiag index +4      4 15
     superdiag index +5      5

subdiagonals have negative index and are:

     subdiag index -1      10 21 32
     subdiag index -2      20 31
     subdiag index -3      30
— Function: int mpgsl_matrix_init_vector_view_over_diagonal (mpgsl_vector_view * v, void * m, int diagonal_index)

Initialise a vector view to reference elements from diagonal diagonal_index of matrix view m.

m is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_vector_view_over_diagonal (mpgsl_vector_view ** vp, void * m, int diagonal_index)

Wrapper for mpgsl_matrix_init_vector_view_over_diagonal() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by vp.

Vector views over antidiagonals

Given the following matrix:

      0  1  2  3  4  5
     10 11 12 13 14 15
     20 21 22 23 24 25
     30 31 32 33 34 35

the antidiagonal of index 0 is the vector 5 14 23 32; the superantidiagonals have positive index and are:

     superantidiag index +1      4 13 22 31
     superantidiag index +2      3 12 21 30
     superantidiag index +3      2 11 20
     superantidiag index +4      1 10
     superantidiag index +5      0

subantidiagonals have negative index and are:

     subdiag index -1      15 24 33
     subdiag index -2      25 34
     subdiag index -3      35
— Function: int mpgsl_matrix_init_vector_view_over_antidiagonal (mpgsl_vector_view * v, void * m, int antidiagonal_index)

Initialise a vector view to reference elements from antidiagonal antidiagonal_index of matrix view m.

m is of type void * so that it can reference both a mpgsl_matrix and a mpgsl_matrix_view structure.

— Function: int mpgsl_matrix_alloc_vector_view_over_antidiagonal (mpgsl_vector_view ** vp, void * m, int antidiagonal_index)

Wrapper for mpgsl_matrix_init_vector_view_over_antidiagonal() that allocates a new vector structure. A pointer to the new structure is stored in the variable referenced by vp.


Previous: compound matrix view vector, Up: compound matrix view
5.5.7.5 Matrix views over vectors

In all the following functions the parameter v is meant to be a pointer to a mpgsl_vector_view or mpgsl_vector structure.

— Function: void mpgsl_matrix_init_row_view_over_vector (mpgsl_matrix_view * m, void * v)

Initialise the view m to represent the vector view v. The matrix has one row and a number of columns equal to the size of the vector.

— Function: void mpgsl_matrix_alloc_row_view_over_vector (mpgsl_matrix_view ** mp, void * v)

Wrapper for mpgsl_matrix_init_row_view_over_vector() that allocates a new matrix structure. A pointer to the new structure is stored in the variable referenced by mp.

— Function: void mpgsl_matrix_init_column_view_over_vector (mpgsl_matrix_view * m, void * v)

Initialise the view m to represent the vector view v. The matrix has one column and a number of rows equal to the size of the vector.

— Function: void mpgsl_matrix_alloc_column_view_over_vector (mpgsl_matrix_view ** mp, void * v)

Wrapper for mpgsl_matrix_init_column_view_over_vector() that allocates a new matrix structure. A pointer to the new structure is stored in the variable referenced by mp.


Next: , Previous: compound matrix view, Up: compound matrix

5.5.8 Copying matricess

— Function: mpgsl_matrix * mpgsl_matrix_clone (const void * v, mp_rnd_t n)

Interpret v as a pointer to a matrix view and build a clone of the represented matrix. mpgsl_matrix_alloc() is used to allocate a new matrix structure and underlying array of MPFR numbers, then the new matrix is initialised with values from the view. The precision of all the numbers in the destination matrix matches the one in the source view. n is used as rounding mode for all the set operations.

— Function: int mpgsl_matrix_copy (void * dst, const void * src, mp_rnd_t n)

Interpret both dst and src as pointers to matrix views and copy the elements from src to dst. The two views must have the same length. The precision of the numbers in dst is set to the one of the numbers from src. n is used as rounding mode for all the set operations. Return a GSL error code.


Next: , Previous: compound matrix copy, Up: compound matrix

5.5.9 Dimension related functions

All the following functions interpret their void * arguments as pointers to matrix views.

— Function: int mpgsl_matrix_same_dimension_p (const void * a, const void * b)
— Function: int mpgsl_matrix_same_dimension3_p (const void * a, const void * b, const void * c)
— Function: int mpgsl_matrix_same_dimension4_p (const void * a, const void * b, const void * c, const void * d)

Return true if the referenced matrix views have the same dimensions.


Next: , Previous: compound matrix dim, Up: compound matrix

5.5.10 Comparison functions

The following functions compare the MPFR numbers in matrix views. Comparing numbers, for details on the equality criteria. All the following functions interpret their void * arguments as pointers to matrix views.

Exact and approximate equality
— Function: int mpgsl_matrix_equal_p (const void * a, const void * b)

Return non–zero if the two matrix views have equal values according to mpfr_cmp(); else return zero.

— Function: int mpgsl_matrix_fcmp_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two matrix views have equal values according to mpgsl_fcmp(), using epsilon tolerance; else return zero.

— Function: int mpgsl_matrix_absdiff_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two matrix views have numbers satisfying the absdiff–equality criterion; else return zero.

— Function: int mpgsl_matrix_reldiff_equal_p (const void * a, const void * b, mpfr_srcptr epsilon)

Return non–zero if the two matrix views have numbers satisfying the reldiff–equality criterion; else return zero.

Other predicates
— Function: int mpgsl_matrix_zero_p (const void * m)
— Function: int mpgsl_matrix_absdiff_zero_p (const void * m, mpfr_srcptr epsilon)

Return non–zero if all the elements of m are zero, respectively, according to mpfr_zero_p() and the absdiff–equality criterion.

The fcmp–equality and reldiff–equality criteria are not suitable to test for approximate zero values.

— Function: int mpgsl_matrix_ispos_p (const void * m)
— Function: int mpgsl_matrix_isneg_p (const void * m)
— Function: int mpgsl_matrix_isnonpos_p (const void * m)
— Function: int mpgsl_matrix_isnonneg_p (const void * m)

These functions return non–zero if all the elements of the matrix m are strictly positive, strictly negative, non–positive or non–negative respectively, and 0 otherwise.


Next: , Previous: compound matrix cmp, Up: compound matrix

5.5.11 Accessing matrix fields

All the m arguments are interpreted as pointers to matrix views. The row and col arguments are row and column indices inside the matrix view.

— Function: int mpgsl_matrix_get (const void * m, size_t row, size_t col, mpfr_ptr * dstp)

Store in the variable referenced by dstp a pointer to the selected element of m. Return a GSL error code.

— Function: int mpgsl_matrix_set (void * m, size_t row, size_t col, mpfr_srcptr x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_ui (void * m, size_t row, size_t col, unsigned long int x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_si (void * m, size_t row, size_t col, long int x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_uj (void * m, size_t row, size_t col, uintmax_t x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_sj (void * m, size_t row, size_t col, intmax_t x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_d (void * m, size_t row, size_t col, double x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_ld (void * m, size_t row, size_t col, long double x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_z (void * m, size_t row, size_t col, mpz_t x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_q (void * m, size_t row, size_t col, mpq_t x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_f (void * m, size_t row, size_t col, mpf_t x, mp_rnd_t r)
— Function: int mpgsl_matrix_set_decimal64 (void * m, size_t row, size_t col, _Decimal64 x, mp_rnd_t r)

Set the selected element to x. Return a GSL error code.

The function with x of type _Decimal64 is available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: void mpgsl_matrix_set_all (void * m, mpfr_srcptr x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_ui (void * m, unsigned long int x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_si (void * m, long int x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_uj (void * m, uintmax_t x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_sj (void * m, intmax_t x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_d (void * m, double x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_ld (void * m, long double x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_z (void * m, mpz_t x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_q (void * m, mpq_t x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_f (void * m, mpf_t x, mp_rnd_t r)
— Function: void mpgsl_matrix_set_all_decimal64 (void * m, _Decimal64 x, mp_rnd_t r)

Set all the values of m to x.

The function with x of type _Decimal64 is available only if the preprocessor symbol MPFR_WANT_DECIMAL_FLOATS is defined by MPFR.

— Function: void mpgsl_matrix_set_all_zero (void * m)

Set all the elements to zero.

— Function: int mpgsl_matrix_set_eye (void * m)

Set to zero all the elements of m, with the exception of the main diagonal elements which are set to one. The matrix must be square. Return a GSL error code.


Next: , Previous: compound matrix access, Up: compound matrix

5.5.12 Finding maximum and minimum elements

— Function: int mpgsl_matrix_min (void * m, mpfr_ptr * minp)
— Function: int mpgsl_matrix_max (void * m, mpfr_ptr * maxp)
— Function: int mpgsl_matrix_minmax (void * m, mpfr_ptr * minp, mpfr_ptr * maxp)

Interpret m as pointer to a matrix view and store: in the variable referenced by minp a pointer to the minimum value, in the variable referenced by maxp a pointer to the maximum value. Return a GSL error code. It is an error if the view holds a NaN element.

— Function: int mpgsl_matrix_abs_min (void * m, mpfr_ptr * minp)
— Function: int mpgsl_matrix_abs_max (void * m, mpfr_ptr * maxp)
— Function: int mpgsl_matrix_abs_minmax (void * m, mpfr_ptr * minp, mpfr_ptr * maxp)

Interpret m as pointer to a matrix view and store: in the variable referenced by minp a pointer to the number whose absolute value is minimum, in the variable referenced by maxp a pointer to the number whose absolute value is maximum. Return a GSL error code. It is an error if the view holds a NaN element.

— Function: int mpgsl_matrix_min_index (void * m, size_t * minrowp, size_t * mincolp)
— Function: int mpgsl_matrix_max_index (void * m, size_t * maxrowp, size_t * maxcolp)
— Function: int mpgsl_matrix_minmax_index (void * m, size_t * minrowp, size_t * mincolp, size_t * maxrowp, size_t * maxcolp)

Interpret m as pointer to a matrix view and store in the variable referenced by:

minrowp
The row index of the minimum value.
mincolp
The column index of the minimum value.
maxrowp
The row index of the maximum value.
maxcolp
The column index of the maximum value.

Return a GSL error code. It is an error if the view holds a NaN element.

— Function: int mpgsl_matrix_abs_min_index (void * m, size_t * minrowp, size_t * mincolp)
— Function: int mpgsl_matrix_abs_max_index (void * m, size_t * maxrowp, size_t * maxcolp)
— Function: int mpgsl_matrix_abs_minmax_index (void * m, size_t * minrowp, size_t * mincolp, size_t * maxrowp, size_t * maxcolp)

Interpret m as pointer to a matrix view and store in the variable referenced by:

minrowp
The row index of the number whose absolute value is minimum.
mincolp
The column index of the number whose absolute value is minimum.
maxrowp
The row index of the number whose absolute value is maximum.
maxcolp
The column index of the number whose absolute value is maximum.

Return a GSL error code. It is an error if the view holds a NaN element.


Next: , Previous: compound matrix minmax, Up: compound matrix

5.5.13 Swapping elements and permutations

All the following functions interpret the arguments m, ma and mb as pointers to matrix views.

— Function: void mpgsl_matrix_transpose (void * m)

Mutate the view to represent the same matrix transposed.

— Function: void mpgsl_matrix_transpose2 (void * dst, void * src)

Initialise the matrix view referenced by dst to represent the tranposition of the matrix view represented by src. The data block is left untouched; the two views references the same data block.

— Function: int mpgsl_matrix_swap (void * ma, void * mb)

Swap the elements of the matrices. The two matrices must have the same length. Return a GSL error code.

— Function: int mpgsl_matrix_swap_elements (void * m, size_t row1, size_t col1, size_t row2, size_t col2)

Swap the element at coordinates (row1, col1) with the element at coordinates (row2, col2). Return a GSL error code.

— Function: int mpgsl_matrix_swap_rows (void * m, size_t row1, size_t row2)
— Function: int mpgsl_matrix_swap_cols (void * m, size_t col1, size_t col2)

Swap the selected rows or columns.

— Function: void mpgsl_matrix_reverse_rows (void * m)
— Function: void mpgsl_matrix_reverse_cols (void * m)

Mutate the matrix view to describe the same matrix with rows or columns reversed.

— Function: int mpgsl_matrix_permute_rows (void * m, gsl_permutation * perm)
— Function: int mpgsl_matrix_permute_cols (void * m, gsl_permutation * perm)

Apply the permutation to the rows or columns of the matrix view. Return a GSL error code.

— Function: int mpgsl_matrix_permute_inverse_rows (void * m, gsl_permutation * perm)
— Function: int mpgsl_matrix_permute_inverse_cols (void * m, gsl_permutation * perm)

Apply the inverse of the permutation to the rows or columns of the matrix view. Return a GSL error code.


Next: , Previous: compound matrix swap, Up: compound matrix

5.5.14 Mapping functions


Next: , Up: compound matrix map
5.5.14.1 Introduction to mapping functions

All the mapping functions map their function parameter f over the elements of the operands to produce a vector or matrix result. The signature of f must be one among the predefined prototypes described in Operations on compound values.

Among the other typical parameters:

All the mapping function names have prefix mpgsl_map_ and as suffix a signature identifying the arguments; examples:

     mpgsl_map_mr__mr_n
     mpgsl_map_mr__mr_vr
     mpgsl_map_mr__mr_mr_no

the signature has the following parts:

             mpgsl_map_mr__mr_mr_no
               ^       ^   ^  ^   ^
               |       |   |  |   |
        -------    ----     --+    --------
       |          |           |            |
     prefix     result     operands     ancillary arguments

The result and operands identifiers can be:

vr
For vectors of real numbers; the corresponding function parameter is interpreted as a vector view.
mr
For matrices of real numbers; the corresponding function parameter is interpreted as a matrix view.
lr
For real number scalars; the corresponding function parameter is a mpfr_ptr or mpfr_srcptr.

The fields in the suffix are:

All the map functions return a GSL error code.


Next: , Previous: compound matrix map intro, Up: compound matrix map
5.5.14.2 Mapping functions for matrices
— Function: int mpgsl_map_mr__mr (mpgsl_r_r_fun_t * f, void * r, const void * a)
— Function: int mpgsl_map_mr__mr_n (mpgsl_r_rn_fun_t * f, void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_map_mr__mr_o (mpgsl_r_ro_fun_t * f, void * r, const void * a, void * custom_context)
— Function: int mpgsl_map_mr__mr_no (mpgsl_r_rno_fun_t * f, void * r, const void * a, mp_rnd_t n, void * custom_context)

Apply an operation to a single matrix operand.

— Function: int mpgsl_map_mr__mr_mr (mpgsl_r_rr_fun_t * f, void * r, const void * a, const void * b)
— Function: int mpgsl_map_mr__mr_mr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * a, const void * b, mp_rnd_t n)
— Function: int mpgsl_map_mr__mr_mr_o (mpgsl_r_rro_fun_t * f, void * r, const void * a, const void * b, void * custom_context)
— Function: int mpgsl_map_mr__mr_mr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * a, const void * b, mp_rnd_t n, void * custom_context)

Apply an operation to two matrix operands.

— Function: int mpgsl_map_mr__mr_mr_mr (mpgsl_r_rrr_fun_t * f, void * r, const void * a, const void * b, const void * c)
— Function: int mpgsl_map_mr__mr_mr_mr_n (mpgsl_r_rrrn_fun_t * f, void * r, const void * a, const void * b, const void * c, mp_rnd_t n)
— Function: int mpgsl_map_mr__mr_mr_mr_o (mpgsl_r_rrro_fun_t * f, void * r, const void * a, const void * b, const void * c, void * custom_context)
— Function: int mpgsl_map_mr__mr_mr_mr_no (mpgsl_r_rrrno_fun_t * f, void * r, const void * a, const void * b, const void * c, mp_rnd_t n, void * custom_context)

Apply an operation to three matrix operands.

Prototypes
— Function Signature: int mpgsl_mr__mr_fun_t (void * r, const void * a)
— Function Signature: int mpgsl_mr__mr_n_fun_t (void * r, const void * a, mp_rnd_t n)
— Function Signature: int mpgsl_mr__mr_o_fun_t (void * r, const void * a, mp_rnd_t n, void * custom_context)
— Function Signature: int mpgsl_mr__mr_no_fun_t (void * r, const void * a, mp_rnd_t n, mp_rnd_t n, void * custom_context)

Perform an operation on a single matrix operand. This is for functions like mpgsl_matrix_round() and mpgsl_matrix_sin().

— Function Signature: int mpgsl_mr__mr_mr_fun_t (void * r, const void * a, const void * b)
— Function Signature: int mpgsl_mr__mr_mr_n_fun_t (void * r, const void * a, const void * b, mp_rnd_t n)
— Function Signature: int mpgsl_mr__mr_mr_o_fun_t (void * r, const void * a, const void * b, void * custom_context)
— Function Signature: int mpgsl_mr__mr_mr_no_fun_t (void * r, const void * a, const void * b, mp_rnd_t n, void * custom_context)

Perform an operation on two matrix operands. This is for functions like mpgsl_matrix_add().

— Function Signature: int mpgsl_mr__mr_mr_mr_fun_t (void * r, const void * a, const void * b, const void * c)
— Function Signature: int mpgsl_mr__mr_mr_mr_n_fun_t (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)
— Function Signature: int mpgsl_mr__mr_mr_mr_o_fun_t (void * r, const void * a, const void * b, const void * c, void * custom_context)
— Function Signature: int mpgsl_mr__mr_mr_mr_no_fun_t (void * r, const void * a, const void * b, const void * c, mp_rnd_t n, void * custom_context)

Perform an operation on three matrix operands. This is for functions like mpgsl_matrix_fma().


Next: , Previous: compound matrix map matrix, Up: compound matrix map
5.5.14.3 Mapping functions for matrices and scalars
— Function: int mpgsl_map_mr__mr_lr (mpgsl_r_rr_fun_t * f, void * r, const void * a, mpfr_srcptr lambda)
— Function: int mpgsl_map_mr__mr_lr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * a, mpfr_srcptr lambda, mp_rnd_t n)
— Function: int mpgsl_map_mr__mr_lr_o (mpgsl_r_rro_fun_t * f, void * r, const void * a, mpfr_srcptr lambda, void * custom_context)
— Function: int mpgsl_map_mr__mr_lr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * a, mpfr_t lambda, mp_rnd_t n, void * custom_context)

Apply an operation to a matrix operand and a scalar operand.

Prototypes
— Function Signature: int mpgsl_mr__mr_lr_fun_t (void * r, const void * a, mpfr_srcptr lambda)
— Function Signature: int mpgsl_mr__mr_lr_n_fun_t (void * r, const void * a, mpfr_srcptr lambda, mp_rnd_t n)
— Function Signature: int mpgsl_mr__mr_lr_o_fun_t (void * r, const void * a, mpfr_srcptr lambda, void * custom_context)
— Function Signature: int mpgsl_mr__mr_lr_no_fun_t (void * r, const void * a, mpfr_srcptr lambda, mp_rnd_t n, void * custom_context)

Perform an operation on a matrix operand and a scalar operand. This is for functions like mpgsl_matrix_add_scalar().


Previous: compound matrix map scalar, Up: compound matrix map
5.5.14.4 Mapping functions for matrices and vectors
— Function: int mpgsl_map_mr__mr_vr (mpgsl_r_rr_fun_t * f, void * r, const void * m, const void * v)
— Function: int mpgsl_map_mr__mr_vr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * m, const void * v, mp_rnd_t n)
— Function: int mpgsl_map_mr__mr_vr_o (mpgsl_r_rro_fun_t * f, void * r, const void * m, const void * v, void * custom_context)
— Function: int mpgsl_map_mr__mr_vr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * m, const void * v, mp_rnd_t n, void * custom_context)

Apply an operation a matrix operand and a scalar operand, yielding a matrix result. The two matrices must have the same dimensions, and either the number of rows or the number of columns must be 1. The number of elements in the sigle row or column must be equal to the number of elements in the vector.

— Function: int mpgsl_map_vr__mr_vr (mpgsl_r_rr_fun_t * f, void * r, const void * m, const void * v)
— Function: int mpgsl_map_vr__mr_vr_n (mpgsl_r_rrn_fun_t * f, void * r, const void * m, const void * v, mp_rnd_t n)
— Function: int mpgsl_map_vr__mr_vr_o (mpgsl_r_rro_fun_t * f, void * r, const void * m, const void * v, void * custom_context)
— Function: int mpgsl_map_vr__mr_vr_no (mpgsl_r_rrno_fun_t * f, void * r, const void * m, const void * v, mp_rnd_t n, void * custom_context)

Apply an operation a matrix operand and a scalar operand, yielding a vector result. The two vectors must have the same size; either the number of rows or the number of columns in the matrix must be 1. The number of elements in the sigle row or column must be equal to the number of elements in the vectors.

Prototypes
— Function: int mpgsl_mr__mr_vr_fun_t (void * r, const void * a, const void * b)
— Function: int mpgsl_mr__mr_vr_n_fun_t (void * r, const void * a, const void * b, mp_rnd_t n)
— Function: int mpgsl_mr__mr_vr_o_fun_t (void * r, const void * a, const void * b, void * custom_context)
— Function: int mpgsl_mr__mr_vr_no_fun_t (void * r, const void * a, const void * b, mp_rnd_t n, void * custom_context)

Perform an operation on a matrix and a vector, yielding a matrix as result.

— Function: int mpgsl_vr__mr_vr_fun_t (void * r, const void * a, const void * b)
— Function: int mpgsl_vr__mr_vr_n_fun_t (void * r, const void * a, const void * b, mp_rnd_t n)
— Function: int mpgsl_vr__mr_vr_o_fun_t (void * r, const void * a, const void * b, void * custom_context)
— Function: int mpgsl_vr__mr_vr_no_fun_t (void * r, const void * a, const void * b, mp_rnd_t n, void * custom_context)

Perform an operation on a matrix and a vector, yielding a vector as result.


Next: , Previous: compound matrix map, Up: compound matrix

5.5.15 Operations on matrices

All the functions described in this section interpret their void* arguments as pointers to vector or matrix views and return a GSL error code. The n argument is used as rounding mode for the operation. For all the functions:


Next: , Up: compound matrix operations
5.5.15.1 Arithmetics between matrices

Basic Arithmetic Functions.

— Function: int mpgsl_matrix_add (void * r, const void * a, const void * b, mp_rnd_t n)

Add the elements from b to the elements from a.

— Function: int mpgsl_matrix_sub (void * r, const void * a, const void * b, mp_rnd_t n)

Subtract the elements from b to the elements from a.

— Function: int mpgsl_matrix_mul (void * r, const void * a, const void * b, mp_rnd_t n)

Multiply the elements from b to the elements from a.

— Function: int mpgsl_matrix_div (void * r, const void * a, const void * b, mp_rnd_t n)

Divide the elements from a by the elements from b.

— Function: int mpgsl_matrix_fma (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)

Fused multiply and add. Multiply a and b element by element, then add the elements from c.

— Function: int mpgsl_matrix_fms (void * r, const void * a, const void * b, const void * c, mp_rnd_t n)

Fused multiply and subtract. Multiply a and b element by element, then subtract the elements from c.

— Function: int mpgsl_matrix_agm (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the arithmetic–geometric mean between the elements from a and the elements from b.

— Function: int mpgsl_matrix_hypot (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the Euclidean norm of the elements from a and the elements from b. This is the square root of the sum between the square of a and the square of b.

— Function: int mpgsl_matrix_sqr (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_sqrt (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_rec_sqrt (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_cbrt (void * r, const void * a, mp_rnd_t n)

Compute the square, square root, reciprocal square root and cubic square root.

— Function: int mpgsl_matrix_root (void * r, const void * a, unsigned long int k, mp_rnd_t rnd)

Compute the k-th root.

— Function: int mpgsl_matrix_pow (void * r, const void * a, const void * b, mp_rnd_t n)

Compute the power of the elements from a raised to the elements of b.

— Function: int mpgsl_matrix_neg (void * r, const void * a, mp_rnd_t n)

Compute the opposite.

— Function: int mpgsl_matrix_abs (void * r, const void * a, mp_rnd_t n)

Compute the absolute value.


Next: , Previous: compound matrix arithmetics matrix, Up: compound matrix operations
5.5.15.2 Arithmetics between a matrix and a scalar
— Function: int mpgsl_matrix_add_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Add the scalar x to all the elements from a.

— Function: int mpgsl_matrix_sub_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Subtract the scalar x from all the elements from a.

— Function: int mpgsl_matrix_mul_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Multiply the scalar x to all the elements from a.

— Function: int mpgsl_matrix_div_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Divide all the elements from a by the scalar x.

— Function: int mpgsl_matrix_pow_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Compute the power of the elements from a raised to the scalar x.

— Function: int mpgsl_matrix_swapdiv_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Divide the scalar x by all the elements from a.

— Function: int mpgsl_matrix_swappow_scalar (void * r, const void * a, mpfr_srcptr x, mp_rnd_t n)

Compute the power of x raised to the elements from a.


Next: , Previous: compound matrix arithmetics scalar, Up: compound matrix operations
5.5.15.3 Algebra between matrices
— Function: int mpgsl_matrix_dot (void * r, const void * a, const void * b, mp_rnd_t n)

Perform a rows by columns product between a and b, and store the result in r. n is used as rounding mode for the multiplication and addition operations.

— Function: int mpgsl_matrix_dot_rr (void * r, const void * a, const void * b, mp_rnd_t n)

Perform a rows by rows product between a and b, and store the result in r. n is used as rounding mode for the multiplication and addition operations.

This is like a rows by columns product between a and the transposed of b.

— Function: int mpgsl_matrix_dot_cc (void * r, const void * a, const void * b, mp_rnd_t n)

Perform a columns by columns product between a and b, and store the result in r. n is used as rounding mode for the multiplication and addition operations.

This is like a rows by columns product between the transposed of a and b.

— Function: int mpgsl_matrix_dot_cr (void * r, const void * a, const void * b, mp_rnd_t n)

Perform a columns by rows product between a and b, and store the result in r. n is used as rounding mode for the multiplication and addition operations.

This is like a columns by rows product between the transposed of a and the transposed of b.


Next: , Previous: compound matrix algebra matrix, Up: compound matrix operations
5.5.15.4 Logarithms and exponentials

Special Functions.

— Function: int mpgsl_matrix_log (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_log2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_log10 (void * r, const void * a, mp_rnd_t n)

Compute the natural logarithm, the base 2 logarithm and the base 10 logarithm.

— Function: int mpgsl_matrix_log1p (void * r, const void * a, mp_rnd_t n)

Compute the sum between 1 and the elements from a, then apply the natural logarithm to the result.

— Function: int mpgsl_matrix_exp (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_exp2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_exp10 (void * r, const void * a, mp_rnd_t n)

Compute the exponential, the power of 2 and the power of 10.

— Function: int mpgsl_matrix_expm1 (void * r, const void * a, mp_rnd_t n)

Compute the exponential of the elements from a, then subtract 1 from the result.


Next: , Previous: compound matrix func log, Up: compound matrix operations
5.5.15.5 Trigonometric functions

Special Functions.

— Function: int mpgsl_matrix_sin (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_cos (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_tan (void * r, const void * a, mp_rnd_t n)

Compute the sine, cosine and tangent functions.

— Function: int mpgsl_matrix_sec (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_csc (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_cot (void * r, const void * a, mp_rnd_t n)

Compute the secant, cosecant and cotangent functions.

— Function: int mpgsl_matrix_asin (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_acos (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_atan (void * r, const void * a, mp_rnd_t n)

Compute the arc sine, arc cosine and arc tangent functions.

— Function: int mpgsl_matrix_atan2 (void * r, const void * y, const void * x, mp_rnd_t n)

Compute the second form of arc tangent of the operands. Special Functions, for details on the expression.


Next: , Previous: compound matrix func trigo, Up: compound matrix operations
5.5.15.6 Hyperbolic functions

Special Functions.

— Function: int mpgsl_matrix_sinh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_cosh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_tanh (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic sine, cosine and tangent functions.

— Function: int mpgsl_matrix_sech (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_csch (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_coth (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic secant, cosecant and cotangent functions.

— Function: int mpgsl_matrix_asinh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_acosh (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_atanh (void * r, const void * a, mp_rnd_t n)

Compute the hyperbolic arc sine, arc cosine and arc tangent functions.


Next: , Previous: compound matrix func hyper, Up: compound matrix operations
5.5.15.7 Special functions

Special Functions.

— Function: int mpgsl_matrix_eint (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_li2 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_gamma (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_zeta (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_erf (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_erfc (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_j0 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_j1 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_y0 (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_y1 (void * r, const void * a, mp_rnd_t n)

Compute miscellaneous special functions implemented by MPFR. Special Functions, for details.

— Function: int mpgsl_matrix_jn (void * r, long k, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_yn (void * r, long k, const void * a, mp_rnd_t n)

Compute the first and second kinds of Bessel functions of order k.


Previous: compound matrix func special, Up: compound matrix operations
5.5.15.8 Integer related functions

Integer Related Functions

— Function: int mpgsl_matrix_rint (void * r, const void * a, mp_rnd_t n)

Round to the nearest representable integer in the given rounding mode.

— Function: int mpgsl_matrix_rint_ceil (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_rint_floor (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_rint_round (void * r, const void * a, mp_rnd_t n)
— Function: int mpgsl_matrix_rint_trunc (void * r, const void * a, mp_rnd_t n)

Round in various directions. See the original MPFR documentation for details.

— Function: int mpgsl_matrix_ceil (void * r, const void * a)
— Function: int mpgsl_matrix_floor (void * r, const void * a)
— Function: int mpgsl_matrix_round (void * r, const void * a)
— Function: int mpgsl_matrix_trunc (void * r, const void * a)

Round in various directions. See the original MPFR documentation for details.

— Function: int mpgsl_matrix_frac (void * r, const void * a, mp_rnd_t n)

Extract the fractional part.


Previous: compound matrix operations, Up: compound matrix

5.5.16 Input and output functions

The following functions are simple wrappers for mpfr_fprintf(); they are meant moslty to produce debugging messages. Formatted Output Functions. No fancy formatting is supported by this API.

The void * arguments are interpreted as pointers to matrix views. The return values are the numbers of characters output to the channel.

— Function: int mpgsl_matrix_fprintf (FILE * file, const char * format, const void * v)

Print the elements from v to file using format for each number. Usage example:

          mpgsl_matrix     matrix;
          
          mpgsl_matrix_fprintf(stderr, "%Rf", &matrix);
— Function: int mpgsl_matrix_printf (const char * format, const void * v)

Print the elements from v to stdout using format.

— Function: int mpgsl_matrix_print (const void * v)

Print the elements from v to stdout using %Rf as format.

— Function: int mpgsl_matrix_print_stderr (const void * v)

Print the elements from v to stderr using %Rf as format.


Next: , Previous: compound, Up: Top

6 One–dimensional root–finding

This chapter describes routines for finding roots of arbitrary one-dimensional functions. The library provides low level components for a variety of iterative solvers and convergence tests. These can be combined by the user to achieve the desired solution, with full access to the intermediate steps of the iteration. Each class of methods uses the same framework, so that we can switch between solvers at runtime without needing to recompile our program. Each instance of a solver keeps track of its own state, allowing the solvers to be used in multi–threaded programs.


Next: , Up: one root

6.1 Overview

One–dimensional root finding algorithms can be divided into two classes: root bracketing and root polishing. Algorithms which proceed by bracketing a root are guaranteed to converge. Bracketing algorithms begin with a bounded region known to contain a root. The size of this bounded region is reduced, iteratively, until it encloses the root to a desired tolerance. This provides a rigorous error estimate for the location of the root.

The technique of root polishing attempts to improve an initial guess to the root. These algorithms converge only if started “close enough” to a root, and sacrifice a rigorous error bound for speed. By approximating the behavior of a function in the vicinity of a root they attempt to find a higher order improvement of an initial guess. When the behavior of the function is compatible with the algorithm and a good initial guess is available a polishing algorithm can provide rapid convergence.

In GSL both types of algorithm are available in similar frameworks. The user provides a high–level driver for the algorithms, and the library provides the individual functions necessary for each of the steps. There are three main phases of the iteration. The steps are:

  1. Initialize solver state, s, for algorithm T.
  2. Update s using the iteration T.
  3. Test s for convergence, and repeat iteration if necessary.

The state for bracketing solvers is held in a mpgsl_root_fsolver struct. The updating procedure uses only function evaluations (not derivatives). The state for root polishing solvers is held in a mpgsl_root_fdfsolver struct. The updates require both the function and its derivative (hence the name fdf) to be supplied by the user.


Next: , Previous: one root overview, Up: one root

6.2 Caveats

Note that root finding functions can only search for one root at a time. When there are several roots in the search area, the first root to be found will be returned; however it is difficult to predict which of the roots this will be. In most cases, no error will be reported if you try to find a root in an area where there is more than one.

Care must be taken when a function may have a multiple root (such as f(x) = (x-x_0)^2 or f(x) = (x-x_0)^3). It is not possible to use root–bracketing algorithms on even–multiplicity roots. For these algorithms the initial interval must contain a zero-crossing, where the function is negative at one end of the interval and positive at the other end. Roots with even–multiplicity do not cross zero, but only touch it instantaneously. Algorithms based on root bracketing will still work for odd–multiplicity roots (e.g. cubic, quintic, ...). Root polishing algorithms generally work with higher multiplicity roots, but at a reduced rate of convergence. In these cases the Steffenson algorithm can be used to accelerate the convergence of multiple roots.

While it is not absolutely required that f have a root within the search region, numerical root finding functions should not be used haphazardly to check for the existence of roots. There are better ways to do this. Because it is easy to create situations where numerical root finders can fail, it is a bad idea to throw a root finder at a function you do not know much about. In general it is best to examine the function visually by plotting before searching for a root.


Next: , Previous: one root caveats, Up: one root

6.3 Initializing the Solver

Function only solvers
— Function: mpgsl_root_fsolver * mpgsl_root_fsolver_alloc (const mpgsl_root_fsolver_type * T)

This function returns a pointer to a newly allocated instance of a solver of type T. For example, the following code creates an instance of a bisection solver:

          const mpgsl_root_fsolver_type * T
            = mpgsl_root_fsolver_bisection;
          mpgsl_root_fsolver * s
            = mpgsl_root_fsolver_alloc (T);

If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of GSL_ENOMEM.

— Function: void mpgsl_root_fsolver_free (mpgsl_root_fsolver * s)

Free all the memory associated with the solver s.

— Function: int mpgsl_root_fsolver_set (mpgsl_root_fsolver * s, mpgsl_function * f, mpfr_t x_lower, mpfr_t x_upper)

This function initializes, or reinitializes, an existing solver s to use the function f and the initial search interval [x_lower, x_upper].

— Function: const char * mpgsl_root_fsolver_name (const mpgsl_root_fsolver * s)

Return a pointer to the name of the solver. Example:

          printf("s is a '%s' solver\n",
                 mpgsl_root_fsolver_name(s));

would print something like s is a 'bisection' solver.

Function and derivative solvers
— Function: mpgsl_root_fdfsolver * mpgsl_root_fdfsolver_alloc (const mpgsl_root_fdfsolver_type * T)

This function returns a pointer to a newly allocated instance of a derivative-based solver of type T. For example, the following code creates an instance of a Newton-Raphson solver:

          const mpgsl_root_fdfsolver_type * T
            = mpgsl_root_fdfsolver_newton;
          mpgsl_root_fdfsolver * s
            = mpgsl_root_fdfsolver_alloc (T);

If there is insufficient memory to create the solver then the function returns a null pointer and the error handler is invoked with an error code of GSL_ENOMEM.

— Function: void mpgsl_root_fdfsolver_free (mpgsl_root_fdfsolver * s)

Free all the memory associated with the solver s.

— Function: int mpgsl_root_fdfsolver_set (mpgsl_root_fdfsolver * s, mpgsl_function_fdf * fdf, mpfr_t root)

This function initializes, or reinitializes, an existing solver s to use the function and derivative fdf and the initial guess root.

— Function: const char * mpgsl_root_fsolver_name (const mpgsl_root_fsolver * s)

Return a pointer to the name of the solver. Example:

          printf ("s is a '%s' solver\n",
                  mpgsl_root_fsolver_name (s));

would print something like s is a 'bisection' solver.


Next: , Previous: one root init, Up: one root

6.4 Providing the function to solve

We must provide a continuous function of one variable for the root finders to operate on, and, sometimes, its first derivative. In order to allow for general parameters the functions are defined by the data types described in this section.


Next: , Up: one root func

6.4.1 Providing only the function

— Data Type: mpgsl_function

This data type defines a general math function with parameters. Public fields description follows.

int (* function) (mpfr_t result, mpfr_t x, void * params)
Pointer to a C function that should store in the (already initialised) number result, the value of the math function computed in x, using the optional custom parameters referenced by params. If successful the return value must be GSL_SUCCESS, else a different code must be returned.
void * params
A pointer to the parameters of the function.

— Macro: MPGSL_FN_EVAL (mpgsl_function * F, mpfr_t result, mpfr_t x)

Apply the math function described by F to the value x and store the result in result. If successful returns GSL_SUCCESS, else a different code is returned.

Example

Here is an example for the general quadratic function: with a = 3, b = 2, c = 1. The following code defines an mpgsl_function F which we could pass to a root finder, then ivokes it and prints the result.

     #include <mpgsl.h>
     
     typedef struct {
       mpfr_t a, b, c;
     } parameters_tag_t;
     typedef parameters_tag_t *      parameters_t;
     
     static void
     parameters_init (parameters_t params)
     {
       mpfr_init(params->a);
       mpfr_init(params->b);
       mpfr_init(params->c);
     
       mpfr_set_d(params->a, 3.0, GMP_RNDN);
       mpfr_set_d(params->b, 2.0, GMP_RNDN);
       mpfr_set_d(params->c, 1.0, GMP_RNDN);
     }
     static void
     parameters_final (parameters_t params)
     {
       mpfr_clear(params->a);
       mpfr_clear(params->b);
       mpfr_clear(params->c);
     }
     
     static int
     function (mpfr_t result, mpfr_t x, void * p)
     {
       mpfr_t        tmp1, tmp2;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp1);
       mpfr_init(tmp2);
       {
         /* (a * x + b) * x + c */
         mpfr_mul(tmp1, params->a, x, GMP_RNDN);
         mpfr_add(tmp2, tmp1, params->b, GMP_RNDN);
         mpfr_mul(tmp1, tmp2, x, GMP_RNDN);
         mpfr_add(result, tmp1, params->c, GMP_RNDN);
       }
       mpfr_clear(tmp1);
       mpfr_clear(tmp2);
       return GSL_SUCCESS;
     }
     
     int
     main (int argc MPGSL_UNUSED, char * argv [])
     {
       parameters_tag_t      params;
       mpfr_t                result, x;
       mpgsl_function        F = {
         .function = function, .params = &params
       };
       double                y = 2.4;
     
     
       parameters_init(&params);
       mpfr_init(result);
       mpfr_init(x);
       {
         mpfr_set_d(x, 2.4, GMP_RNDN);
         MPGSL_FN_EVAL(&F, result, x);
         mpfr_fprintf(stderr,
           "%s: the function value is: %20RNf, should be close to: %f\n",
           argv[0], result, (3.0 * y * y + 2.0 * y + 1.0));
       }
       mpfr_clear(x);
       mpfr_clear(result);
       parameters_final(&params);
       exit(EXIT_SUCCESS);
     }


Previous: one root func only, Up: one root func

6.4.2 Providing the function and its derivative

— Data Type: mpgsl_function_fdf

This data type defines a general math function with parameters and its first derivative. Public fields description follows.

int (* f) (mpfr_t y, mpfr_t x, void * params)
Pointer to a C function that should store in the (already initialised) number y, the value of the math function computed in x, using the optional custom parameters referenced by params. If successful the return value must be GSL_SUCCESS, else a different code must be returned.
int (* df) (mpfr_t dy, mpfr_t x, void * params)
Pointer to a C function that should store in the (already initialised) number dy, the value of the math derivative computed in x, using the optional custom parameters referenced by params. If successful the return value must be GSL_SUCCESS, else a different code must be returned.
int (* fdf) (mpfr_t x, void * params, mpfr_t y, mpfr_t dy)
Pointer to a C function that should store in the (already initialised) numbers y and dy, the values of the math function and its derivative (respectively) computed in x, using the optional custom parameters referenced by params. If successful the return value must be GSL_SUCCESS, else a different code must be returned.
void * params
Pointer to the parameters of the function.

— Macro: MPGSL_FN_FDF_EVAL_F (mpgsl_function_fdf * FDF, mpfr_t y, mpfr_t x)

Apply the math function described by FDF to the value x and store the result in y. If successful returns GSL_SUCCESS, else a different code is returned.

— Macro: MPGSL_FN_FDF_EVAL_DF (mpgsl_function_fdf * FDF, mpfr_t dy, mpfr_t x)

Apply the math derivative described by FDF to the value x and store the result in dy. If successful returns GSL_SUCCESS, else a different code is returned.

— Macro: MPGSL_FN_FDF_EVAL_F_DF (mpgsl_function_fdf * FDF, mpfr_t x, mpfr_t y, mpfr_t dy)

Apply the math function and its derivative described by FDF to the value x and store the result in y and dy. If successful returns GSL_SUCCESS, else a different code is returned.

Example

Here is an example for the general quadratic function, with a = 3, b = 2, c = 1, and its derivative: The following code defines an mpgsl_function_fdf FDF which we could pass to a root finder, then ivokes it and prints the result.

     #include <mpgsl.h>
     
     typedef struct {
       mpfr_t a, b, c;
     } parameters_tag_t;
     typedef parameters_tag_t *      parameters_t;
     
     static void
     parameters_init (parameters_t params)
     {
       mpfr_init(params->a);
       mpfr_init(params->b);
       mpfr_init(params->c);
     
       mpfr_set_d(params->a, 3.0, GMP_RNDN);
       mpfr_set_d(params->b, 2.0, GMP_RNDN);
       mpfr_set_d(params->c, 1.0, GMP_RNDN);
     }
     static void
     parameters_final (parameters_t params)
     {
       mpfr_clear(params->a);
       mpfr_clear(params->b);
       mpfr_clear(params->c);
     }
     
     static int
     function (mpfr_t y, mpfr_t x, void * p)
     {
       mpfr_t        tmp1, tmp2;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp1);
       mpfr_init(tmp2);
       {
         /* (a * x + b) * x + c */
         mpfr_mul(tmp1, params->a, x, GMP_RNDN);
         mpfr_add(tmp2, tmp1, params->b, GMP_RNDN);
         mpfr_mul(tmp1, tmp2, x, GMP_RNDN);
         mpfr_add(y, tmp1, params->c, GMP_RNDN);
       }
       mpfr_clear(tmp1);
       mpfr_clear(tmp2);
       return GSL_SUCCESS;
     }
     static int
     derivative (mpfr_t dy, mpfr_t x, void * p)
     {
       mpfr_t        tmp1, tmp2;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp1);
       mpfr_init(tmp2);
       {
         /* 2 * a * x + b */
         mpfr_mul_d(tmp1, params->a, 2.0, GMP_RNDN);
         mpfr_mul(tmp2, tmp1, x, GMP_RNDN);
         mpfr_add(dy, tmp2, params->b, GMP_RNDN);
       }
       mpfr_clear(tmp1);
       mpfr_clear(tmp2);
       return GSL_SUCCESS;
     }
     static void
     function_and_derivative (mpfr_t x , void * p, mpfr_t y, mpfr_t dy)
     {
       int   e;
       e = function(y, x, p);
       if (GSL_SUCCESS != e) return e;
       return derivative(dy, x, p);
     }
     
     int
     main (int argc MPGSL_UNUSED, char * argv [])
     {
       parameters_tag_t      params;
       mpfr_t                dy, y, x;
       mpgsl_function_fdf    FDF = {
         .f          = function,
         .df         = derivative,
         .fdf        = function_and_derivative,
         .params     = &params
       };
       double                xd = 2.4;
     
     
       parameters_init(&params);
       mpfr_init(dy);
       mpfr_init(y);
       mpfr_init(x);
       {
         mpfr_set_d(x, 2.4, GMP_RNDN);
         MPGSL_FN_FDF_EVAL_F_DF(&FDF, x, y, dy);
         mpfr_fprintf(stderr,
                      "%s: y = %20RNf, dy = %20RNf, should be %f and %f\n",
                      argv[0], y, dy,
                      (3.0 * xd * xd + 2.0 * xd + 1.0),
                      (2.0 * 3.0 * xd + 2.0));
         MPGSL_FN_FDF_EVAL_F(&FDF, y, x);
         mpfr_fprintf(stderr, "%s: y  = %20RNf\n", argv[0], y);
         MPGSL_FN_FDF_EVAL_DF(&FDF, dy, x);
         mpfr_fprintf(stderr, "%s: dy = %20RNf\n", argv[0], dy);
       }
       mpfr_clear(x);
       mpfr_clear(y);
       mpfr_clear(dy);
       parameters_final(&params);
       exit(EXIT_SUCCESS);
     }

Here is another example for the function: and its derivative: in which we take advantage of the math function and derivative expressions to save some computation of intermediate results (see the C function function_and_derivative()).

     #include <mpgsl.h>
     
     typedef struct {
       mpfr_t        ell;
     } parameters_tag_t;
     typedef parameters_tag_t *      parameters_t;
     
     static void
     parameters_init (parameters_t params)
     {
       mpfr_init(params->ell);
       mpfr_set_d(params->ell, 2.0, GMP_RNDN);
     }
     static void
     parameters_final (parameters_t params)
     {
       mpfr_clear(params->ell);
     }
     
     static int
     function (mpfr_t y, mpfr_t x, void * p)
     {
       mpfr_t        tmp;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp);
       {
         /* exp(ell * x) */
         mpfr_mul(tmp, params->ell, x, GMP_RNDN);
         mpfr_exp(y, tmp, GMP_RNDN);
       }
       mpfr_clear(tmp);
       return GSL_SUCCESS;
     }
     static int
     derivative (mpfr_t dy, mpfr_t x, void * p)
     {
       mpfr_t        tmp1, tmp2;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp1);
       mpfr_init(tmp2);
       {
         /* ell * exp(ell * x) */
         mpfr_mul(tmp1, params->ell, x, GMP_RNDN);
         mpfr_exp(tmp2, tmp1, GMP_RNDN);
         mpfr_mul(dy, params->ell, tmp2, GMP_RNDN);
       }
       mpfr_clear(tmp1);
       mpfr_clear(tmp2);
       return GSL_SUCCESS;
     }
     static int
     function_and_derivative (mpfr_t x , void * p, mpfr_t y, mpfr_t dy)
     {
       mpfr_t        tmp;
       parameters_t  params = (parameters_t)p;
     
       mpfr_init(tmp);
       {
         /*
           y  = exp(ell * x)
           dy = ell * exp(ell * x) = ell * y
         */
         mpfr_mul(tmp, params->ell, x, GMP_RNDN);
         mpfr_exp(y, tmp, GMP_RNDN);
         mpfr_mul(dy, params->ell, y, GMP_RNDN);
       }
       mpfr_clear(tmp);
       return GSL_SUCCESS;
     }
     
     int
     main (int argc MPGSL_UNUSED, char * argv [])
     {
       parameters_tag_t      params;
       mpfr_t                dy, y, x;
       mpgsl_function_fdf    FDF = {
         .f          = function,
         .df         = derivative,
         .fdf        = function_and_derivative,
         .params     = &params
       };
       double                xd = 2.4;
     
     
       parameters_init(&params);
       mpfr_init(dy);
       mpfr_init(y);
       mpfr_init(x);
       {
         mpfr_set_d(x, 2.4, GMP_RNDN);
         MPGSL_FN_FDF_EVAL_F_DF(&FDF, x, y, dy);
         mpfr_fprintf(stderr,
                      "%s: y = %20RNf, dy = %20RNf, should be %f and %f\n",
                      argv[0], y, dy,
                      (exp(2.0 * xd)),
                      (2.0 * exp(2.0 * xd)));
         MPGSL_FN_FDF_EVAL_F(&FDF, y, x);
         mpfr_fprintf(stderr, "%s: y  = %20RNf\n", argv[0], y);
         MPGSL_FN_FDF_EVAL_DF(&FDF, dy, x);
         mpfr_fprintf(stderr, "%s: dy = %20RNf\n", argv[0], dy);
       }
       mpfr_clear(x);
       mpfr_clear(y);
       mpfr_clear(dy);
       parameters_final(&params);
       exit(EXIT_SUCCESS);
     }


Next: , Previous: one root func, Up: one root

6.5 Search Bounds and Guesses

We provide either search bounds or an initial guess; this section explains how search bounds and guesses work and how function arguments control them.

A guess is simply an x value which is iterated until it is within the desired precision of a root. It takes the form of an mpfr_t.

Search bounds are the endpoints of an interval which is iterated until the length of the interval is smaller than the requested precision. The interval is defined by two values, the lower limit and the upper limit. Whether the endpoints are intended to be included in the interval or not depends on the context in which the interval is used.


Next: , Previous: one root bounds, Up: one root

6.6 Root Finding Iteration

The following functions drive the iteration of each algorithm. Each function performs one iteration to update the state of any solver of the corresponding type. The same functions work for all solvers so that different methods can be substituted at runtime without modifications to the code.

— Function: int mpgsl_root_fsolver_iterate (mpgsl_root_fsolver * s)
— Function: int mpgsl_root_fdfsolver_iterate (mpgsl_root_fdfsolver * s)

These functions perform a single iteration of the solver s. If the iteration encounters an unexpected problem then an error code will be returned:

GSL_EBADFUNC
The iteration encountered a singular point where the function or its derivative evaluated to Inf or NaN.
GSL_EZERODIV
The derivative of the function vanished at the iteration point, preventing the algorithm from continuing without a division by zero.

The solver maintains a current best estimate of the root at all times. The bracketing solvers also keep track of the current best interval bounding the root. This information can be accessed with the following auxiliary functions.

— Function: mpfr_ptr mpgsl_root_fsolver_root (const mpgsl_root_fsolver * s)
— Function: mpfr_ptr mpgsl_root_fdfsolver_root (const mpgsl_root_fdfsolver * s)

These functions return the current estimate of the root for the solver s. The referenced number is part of solver's state, so it must not be modified.

— Function: mpfr_ptr mpgsl_root_fsolver_x_lower (const mpgsl_root_fsolver * s)
— Function: mpfr_ptr mpgsl_root_fsolver_x_upper (const mpgsl_root_fsolver * s)

These functions return the current bracketing interval for the solver s. The referenced numbers are part of solver's state, so they must not be modified.


Next: , Previous: one root iteration, Up: one root

6.7 Search Stopping Parameters

A root finding procedure should stop when one of the following conditions is true:

The handling of these conditions is under user control. The functions below allow the user to test the precision of the current result in several standard ways.

— Function: int mpgsl_root_test_interval (mpfr_t x_lower, mpfr_t x_upper, mpfr_t epsabs, mpfr_t epsrel)

This function tests for the convergence of the interval [x_lower, x_upper] with absolute error epsabs and relative error epsrel. The test returns GSL_SUCCESS if the following condition is achieved: when the interval x = [a,b] does not include the origin. If the interval includes the origin then \min(|a|,|b|) is replaced by zero (which is the minimum value of |x| over the interval). This ensures that the relative error is accurately estimated for roots close to the origin.

This condition on the interval also implies that any estimate of the root r in the interval satisfies the same condition with respect to the true root r^*: assuming that the true root r^* is contained within the interval.

— Function: int mpgsl_root_test_delta (mpfr_t x1, mpfr_t x0, mpfr_t epsabs, mpfr_t epsrel)

Test for the convergence of the sequence ..., x0, x1 with absolute error epsabs and relative error epsrel. The test returns GSL_SUCCESS if the following condition is achieved: and returns GSL_CONTINUE otherwise.

— Function: int mpgsl_root_test_residual (mpfr_t f, mpfr_t epsabs)

Test the residual value f against the absolute error bound epsabs. The test returns GSL_SUCCESS if the following condition is achieved: and returns GSL_CONTINUE otherwise. This criterion is suitable for situations where the precise location of the root, x, is unimportant provided a value can be found where the residual, |f(x)|, is small enough.


Next: , Previous: one root stop, Up: one root

6.8 Root Bracketing Algorithms

The root bracketing algorithms described in this section require an initial interval which is guaranteed to contain a root—if a and b are the endpoints of the interval then f(a) must differ in sign from f(b). This ensures that the function crosses zero at least once in the interval. If a valid initial interval is used then these algorithm cannot fail, provided the function is well-behaved.

Note that a bracketing algorithm cannot find roots of even degree, since these do not cross the x-axis.

— Solver: mpgsl_root_fsolver_bisection

The bisection algorithm is the simplest method of bracketing the roots of a function. It is the slowest algorithm provided by the library, with linear convergence.

On each iteration, the interval is bisected and the value of the function at the midpoint is calculated. The sign of this value is used to determine which half of the interval does not contain a root. That half is discarded to give a new, smaller interval containing the root. This procedure can be continued indefinitely until the interval is sufficiently small.

At any time the current estimate of the root is taken as the midpoint of the interval.

— Solver: mpgsl_root_fsolver_falsepos

The false position algorithm is a method of finding roots based on linear interpolation. Its convergence is linear, but it is usually faster than bisection.

On each iteration a line is drawn between the endpoints (a,f(a)) and (b,f(b)) and the point where this line crosses the x-axis taken as a “midpoint”. The value of the function at this point is calculated and its sign is used to determine which side of the interval does not contain a root. That side is discarded to give a new, smaller interval containing the root. This procedure can be continued indefinitely until the interval is sufficiently small.

The best estimate of the root is taken from the linear interpolation of the interval on the current iteration.

— Solver: mpgsl_root_fsolver_brent

The Brent-Dekker method (referred to here as Brent's method) combines an interpolation strategy with the bisection algorithm. This produces a fast algorithm which is still robust.

On each iteration Brent's method approximates the function using an interpolating curve. On the first iteration this is a linear interpolation of the two endpoints. For subsequent iterations the algorithm uses an inverse quadratic fit to the last three points, for higher accuracy. The intercept of the interpolating curve with the x-axis is taken as a guess for the root. If it lies within the bounds of the current interval then the interpolating point is accepted, and used to generate a smaller interval. If the interpolating point is not accepted then the algorithm falls back to an ordinary bisection step.

The best estimate of the root is taken from the most recent interpolation or bisection.


Next: , Previous: one root bracketing, Up: one root

6.9 Root Finding Algorithms using Derivatives

The root polishing algorithms described in this section require an initial guess for the location of the root. There is no absolute guarantee of convergence; the function must be suitable for this technique and the initial guess must be sufficiently close to the root for it to work. When these conditions are satisfied then convergence is quadratic.

These algorithms make use of both the function and its derivative.

— Derivative Solver: mpgsl_root_fdfsolver_newton

Newton's Method is the standard root–polishing algorithm. The algorithm begins with an initial guess for the location of the root. On each iteration, a line tangent to the function f is drawn at that position. The point where this line crosses the x-axis becomes the new guess. The iteration is defined by the following sequence: Newton's method converges quadratically for single roots, and linearly for multiple roots.

— Derivative Solver: mpgsl_root_fdfsolver_secant

The secant method is a simplified version of Newton's method which does not require the computation of the derivative on every step.

On its first iteration the algorithm begins with Newton's method, using the derivative to compute a first step: Subsequent iterations avoid the evaluation of the derivative by replacing it with a numerical estimate, the slope of the line through the previous two points: When the derivative does not change significantly in the vicinity of the root the secant method gives a useful saving. Asymptotically the secant method is faster than Newton's method whenever the cost of evaluating the derivative is more than 0.44 times the cost of evaluating the function itself. As with all methods of computing a numerical derivative the estimate can suffer from cancellation errors if the separation of the points becomes too small.

On single roots, the method has a convergence of order (1 + \sqrt 5)/2 (approximately 1.62). It converges linearly for multiple roots.

— Derivative Solver: mpgsl_root_fdfsolver_steffenson

The Steffenson Method provides the fastest convergence of all the routines. It combines the basic Newton algorithm with an Aitken “delta–squared” acceleration. If the Newton iterates are x_i then the acceleration procedure generates a new sequence R_i: which converges faster than the original sequence under reasonable conditions. The new sequence requires three terms before it can produce its first value so the method returns accelerated values on the second and subsequent iterations. On the first iteration it returns the ordinary Newton estimate. The Newton iterate is also returned if the denominator of the acceleration term ever becomes zero.

As with all acceleration procedures this method can become unstable if the function is not well–behaved.


Next: , Previous: one root polishing, Up: one root

6.10 Examples

Root bracketing algorithms

Root polishing algorithms


Next: , Up: one root examples

6.10.1 Bisection algorithm with interval stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fsolver * solver;
       mpfr_t                x_lower, x_upper;
       mpfr_t                epsabs, epsrel;
       int                   e, iteration_count = 0;
       mpgsl_function        F = { .function = sine_function,
                                   .params   = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tbisection algorithm,\n\
     \tinterval stop criterion.\n");
     
       solver = mpgsl_root_fsolver_alloc(mpgsl_root_fsolver_bisection);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(x_lower);
       mpfr_init(x_upper);
       mpfr_init(epsabs);
       mpfr_init(epsrel);
       {
         mpfr_set_d(x_lower, -1.0, GMP_RNDN);
         mpfr_set_d(x_upper, +0.5, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
         mpfr_set_d(epsrel,  1e-3, GMP_RNDN);
     
         e = mpgsl_root_fsolver_set(solver, &F, x_lower, x_upper);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         mpfr_printf("iteration %2d: [%-26Rf, %-26Rf]\n",
                     iteration_count++,
                     mpgsl_root_fsolver_x_lower(solver),
                     mpgsl_root_fsolver_x_upper(solver));
     
         do {
           e = mpgsl_root_fsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           mpfr_printf("iteration %2d: [%-26Rf, %-26Rf]\n",
                       iteration_count++,
                       mpgsl_root_fsolver_x_lower(solver),
                       mpgsl_root_fsolver_x_upper(solver));
     
           e = mpgsl_root_test_interval(mpgsl_root_fsolver_x_lower(solver),
                                        mpgsl_root_fsolver_x_upper(solver),
                                        epsabs, epsrel);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fsolver_root(solver));
       }
      end:
       mpfr_clear(epsrel);
       mpfr_clear(epsabs);
       mpfr_clear(x_upper);
       mpfr_clear(x_lower);
       mpgsl_root_fsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Next: , Previous: one root examples bisection, Up: one root examples

6.10.2 Falsepos algorithm with delta stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fsolver *  solver;
       mpfr_t                x_lower, x_upper;
       mpfr_t                epsabs, epsrel, x1;
       int                   e, iteration_count = 0;
       mpgsl_function        F = { .function = sine_function,
                                   .params   = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tfalsepos algorithm,\n\
     \tdelta stop criterion.\n");
     
       solver = mpgsl_root_fsolver_alloc(mpgsl_root_fsolver_falsepos);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(x_lower);
       mpfr_init(x_upper);
       mpfr_init(epsabs);
       mpfr_init(epsrel);
       mpfr_init(x1);
       {
         mpfr_set_d(x_lower, -1.0, GMP_RNDN);
         mpfr_set_d(x_upper, +0.5, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
         mpfr_set_d(epsrel,  1e-3, GMP_RNDN);
     
         mpfr_set(x1, x_lower, GMP_RNDN);
     
         e = mpgsl_root_fsolver_set(solver, &F, x_lower, x_upper);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         mpfr_printf("iteration %2d: x = %-28Rf\n",
                     iteration_count++, mpgsl_root_fsolver_root(solver));
     
         do {
           e = mpgsl_root_fsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           mpfr_printf("iteration %2d: x = %-28Rf\n",
                       iteration_count++, mpgsl_root_fsolver_root(solver));
     
           e = mpgsl_root_test_delta(x1, mpgsl_root_fsolver_root(solver),
                                     epsabs, epsrel);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             mpfr_set(x1, mpgsl_root_fsolver_root(solver), GMP_RNDN);
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fsolver_root(solver));
       }
      end:
       mpfr_clear(epsrel);
       mpfr_clear(epsabs);
       mpfr_clear(x_upper);
       mpfr_clear(x_lower);
       mpgsl_root_fsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Next: , Previous: one root examples falsepos, Up: one root examples

6.10.3 Brent algorithm with residual stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fsolver *  solver;
       mpfr_t                x_lower, x_upper;
       mpfr_t                epsabs, residual;
       int                   e, iteration_count = 0;
       mpgsl_function        F = { .function = sine_function,
                                   .params   = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tbrent algorithm,\n\
     \tresidual stop criterion.\n");
     
       solver = mpgsl_root_fsolver_alloc(mpgsl_root_fsolver_brent);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(x_lower);
       mpfr_init(x_upper);
       mpfr_init(epsabs);
       mpfr_init(residual);
       {
         mpfr_set_d(x_lower, -1.0, GMP_RNDN);
         mpfr_set_d(x_upper, +0.5, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
     
         e = mpgsl_root_fsolver_set(solver, &F, x_lower, x_upper);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         MPGSL_FN_EVAL(&F, residual, mpgsl_root_fsolver_root(solver));
         mpfr_printf("iteration %2d: x = %-30Rf, f(x) = %-30Rf\n",
                     iteration_count++,
                     mpgsl_root_fsolver_root(solver),
                     residual);
     
         do {
           e = mpgsl_root_fsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           MPGSL_FN_EVAL(&F, residual, mpgsl_root_fsolver_root(solver));
           mpfr_printf("iteration %2d: x = %-30Rf, f(x) = %-30Rf\n",
                       iteration_count++,
                       mpgsl_root_fsolver_root(solver),
                       residual);
     
           e = mpgsl_root_test_residual(residual, epsabs);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fsolver_root(solver));
       }
      end:
       mpfr_clear(residual);
       mpfr_clear(epsabs);
       mpfr_clear(x_upper);
       mpfr_clear(x_lower);
       mpgsl_root_fsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Next: , Previous: one root examples brent, Up: one root examples

6.10.4 Newton algorithm with delta stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     cosine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_cos(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     sine_and_cosine_function (mpfr_t x, void * p MPGSL_UNUSED,
                               mpfr_t y, mpfr_t dy)
     {
       mpfr_sin(y,  x, GMP_RNDN);
       mpfr_cos(dy, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fdfsolver *solver;
       mpfr_t                epsabs, epsrel, x1;
       int                   e, iteration_count = 0;
       mpgsl_function_fdf    F = { .f      = sine_function,
                                   .df     = cosine_function,
                                   .fdf    = sine_and_cosine_function,
                                   .params = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tnewton algorithm,\n\
     \tdelta stop criterion.\n");
     
       solver = mpgsl_root_fdfsolver_alloc(mpgsl_root_fdfsolver_newton);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(epsabs);
       mpfr_init(epsrel);
       mpfr_init(x1);
       {
         mpfr_set_d(x1,     -1.0, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
         mpfr_set_d(epsrel,  1e-6, GMP_RNDN);
     
         e = mpgsl_root_fdfsolver_set(solver, &F, x1);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         mpfr_printf("iteration %2d: x = %-32Rf\n",
                     iteration_count++,
                     mpgsl_root_fdfsolver_root(solver));
     
         do {
           e = mpgsl_root_fdfsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           mpfr_printf("iteration %2d: x = %-32Rf\n",
                       iteration_count++,
                       mpgsl_root_fdfsolver_root(solver));
     
           e = mpgsl_root_test_delta(x1, mpgsl_root_fdfsolver_root(solver),
                                     epsabs, epsrel);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             mpfr_set(x1, mpgsl_root_fdfsolver_root(solver), GMP_RNDN);
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fdfsolver_root(solver));
       }
      end:
       mpfr_clear(epsabs);
       mpfr_clear(epsrel);
       mpfr_clear(x1);
       mpgsl_root_fdfsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Next: , Previous: one root examples newton, Up: one root examples

6.10.5 Secant algorithm with delta stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     cosine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_cos(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     sine_and_cosine_function (mpfr_t x, void * p MPGSL_UNUSED,
                               mpfr_t y, mpfr_t dy)
     {
       mpfr_sin(y,  x, GMP_RNDN);
       mpfr_cos(dy, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fdfsolver *solver;
       mpfr_t                epsabs, epsrel, x1;
       int                   e, iteration_count = 0;
       mpgsl_function_fdf    F = { .f      = sine_function,
                                   .df     = cosine_function,
                                   .fdf    = sine_and_cosine_function,
                                   .params = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tsecant algorithm,\n\
     \tdelta stop criterion.\n");
     
       solver = mpgsl_root_fdfsolver_alloc(mpgsl_root_fdfsolver_secant);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(epsabs);
       mpfr_init(epsrel);
       mpfr_init(x1);
       {
         mpfr_set_d(x1,     -1.0, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
         mpfr_set_d(epsrel,  1e-6, GMP_RNDN);
     
         e = mpgsl_root_fdfsolver_set(solver, &F, x1);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         mpfr_printf("iteration %2d: x = %-32Rf\n",
                     iteration_count++,
                     mpgsl_root_fdfsolver_root(solver));
     
         do {
           e = mpgsl_root_fdfsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           mpfr_printf("iteration %2d: x = %-32Rf\n",
                       iteration_count++,
                       mpgsl_root_fdfsolver_root(solver));
     
           e = mpgsl_root_test_delta(x1, mpgsl_root_fdfsolver_root(solver),
                                     epsabs, epsrel);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             mpfr_set(x1, mpgsl_root_fdfsolver_root(solver), GMP_RNDN);
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fdfsolver_root(solver));
       }
      end:
       mpfr_clear(epsabs);
       mpfr_clear(epsrel);
       mpfr_clear(x1);
       mpgsl_root_fdfsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Previous: one root examples secant, Up: one root examples

6.10.6 Steffenson algorithm with residual stop criterion

     #include <stdlib.h>
     #include <mpgsl.h>
     
     static int
     sine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_sin(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     cosine_function (mpfr_t y, mpfr_t x, void * p MPGSL_UNUSED)
     {
       mpfr_cos(y, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     static int
     sine_and_cosine_function (mpfr_t x, void * p MPGSL_UNUSED,
                               mpfr_t y, mpfr_t dy)
     {
       mpfr_sin(y,  x, GMP_RNDN);
       mpfr_cos(dy, x, GMP_RNDN);
       return GSL_SUCCESS;
     }
     
     int
     main (void)
     {
       mpgsl_root_fdfsolver *solver;
       mpfr_t                guess, epsabs, residual;
       int                   e, iteration_count = 0;
       mpgsl_function_fdf    F = { .f      = sine_function,
                                   .df     = cosine_function,
                                   .fdf    = sine_and_cosine_function,
                                   .params = NULL };
     
       printf("*** One-dimensional root finding:\n\
     \tsteffenson algorithm,\n\
     \tresidual stop criterion.\n");
     
       solver = mpgsl_root_fdfsolver_alloc(mpgsl_root_fdfsolver_steffenson);
       if (NULL == solver) {
         perror("error initialising solver");
         exit(EXIT_FAILURE);
       }
     
       mpfr_init(guess);
       mpfr_init(epsabs);
       mpfr_init(residual);
       {
         mpfr_set_d(guess, -1.0, GMP_RNDN);
         mpfr_set_d(epsabs,  1e-6, GMP_RNDN);
     
         e = mpgsl_root_fdfsolver_set(solver, &F, guess);
         if (GSL_SUCCESS != e) {
           fprintf(stderr, "error setting up: %s\n", gsl_strerror(e));
           goto end;
         }
     
         MPGSL_FN_FDF_EVAL_F(&F, residual,
                             mpgsl_root_fdfsolver_root(solver));
         mpfr_printf("iteration %2d: x = %-32Rf, f(x) = %-32Rf\n",
                     iteration_count++,
                     mpgsl_root_fdfsolver_root(solver),
                     residual);
     
         do {
           e = mpgsl_root_fdfsolver_iterate(solver);
           if (GSL_SUCCESS != e) {
             fprintf(stderr, "error iterating: %s\n", gsl_strerror(e));
             goto end;
           }
     
           MPGSL_FN_FDF_EVAL_F(&F, residual,
                               mpgsl_root_fdfsolver_root(solver));
           mpfr_printf("iteration %2d: x = %-32Rf, f(x) = %-32Rf\n",
                       iteration_count++,
                       mpgsl_root_fdfsolver_root(solver),
                       residual);
     
           e = mpgsl_root_test_residual(residual, epsabs);
           switch (e) {
           case GSL_SUCCESS:
             goto solved;
           case GSL_CONTINUE:
             break;
           default:
             fprintf(stderr, "error testing: %s\n", gsl_strerror(e));
             goto end;
           }
         } while (GSL_CONTINUE == e);
       solved:
         mpfr_printf("result = %Rf\n", mpgsl_root_fdfsolver_root(solver));
       }
      end:
       mpfr_clear(residual);
       mpfr_clear(epsabs);
       mpfr_clear(guess);
       mpgsl_root_fdfsolver_free(solver);
       exit(EXIT_SUCCESS);
     }


Previous: one root examples, Up: one root

6.11 References and Further Reading

For information on the Brent-Dekker algorithm see the following two papers:


Next: , Previous: one root, Up: Top

Appendix A How to interpret the library version

NOTE This section is not specific to this package; it is an explanation of library interface version numbering that the author uses for all its C libraries.

This package is distributing one or more libraries; this package has a version number and each library has an interface version number, package and interface versions are independent.

A.1 Package version

The package version tracks the development of the source code: every time the source code is modified the package version is incremented. There are three numbers: major, minor, patch level.

major
Incremented before big rewritings or additions.
minor
Incremented to track incremental development of the code.
patch level
For stable versions: incremented for bug fixes and very small additions that were forgotten in the last minor update.

The package has also a letter that indicates the status of the release: a for experimental or “alpha”, b for testing or “beta”, . for stable.

A full version number looks like this:

1.2a3
major 1, minor 2, patch level 3, status alpha;
10.23b34
major 10, minor 23, patch level 34, status beta;
1.2.3
major 1, minor 2, patch level 3, status stable.

The reason for the letter to be between the minor number and the patch level is that it is more important:

summary:

     1.1.x < 1.2a0 < 1.2a1 < ... < 1.2b0 < 1.2b1 < ... < 1.2.0 < 1.2.1

A.2 Interface version

The interface version changes only when the public interface of a library (public function prototypes and data types) is modified. There are two numbers: major and minor.

major
Establishes a set of: functions, data types, code behaviour that will not change until the next major number update. Code written for the first release of the library with a major interface number (1.0), will work unchanged with all the libraries with the same major interface number (1.1, 1.2, ...); if it does not: it is a library bug, please signal this to the author using the bug tracking system.

The only exception to this rule is the wrong behaviour caused by bugs; code that relies on buggy behaviour of a package version may not work with subsequent package versions with the same major interface number.

minor
Establishes a set of: functions, data types, code behaviour that will not change until the next major number update. This set is added to the set of the previous interface version.

The library file installed on our system has name composed with the interface number; example: if the package has version 1.2.3 and the interface has version 1.1, the library file is called libmy1.1.so.

Even if the interface number has not changed, and the library file has the same name: when an installed shared library is updated to a new package version, programs and libraries depending on it will have to be recompiled, because the internals of the library file will have changed.


Next: , Previous: library version, Up: Top

Appendix B GNU General Public License

Version 3, 29 June 2007
     Copyright © 2007 Free Software Foundation, Inc. http://fsf.org/
     
     Everyone is permitted to copy and distribute verbatim copies of this
     license document, but changing it is not allowed.

Preamble

The GNU General Public License is a free, copyleft license for software and other kinds of works.

The licenses for most software and other practical works are designed to take away your freedom to share and change the works. By contrast, the GNU General Public License is intended to guarantee your freedom to share and change all versions of a program—to make sure it remains free software for all its users. We, the Free Software Foundation, use the GNU General Public License for most of our software; it applies also to any other work released this way by its authors. You can apply it to your programs, too.

When we speak of free software, we are referring to freedom, not price. Our General Public Licenses are designed to make sure that you have the freedom to distribute copies of free software (and charge for them if you wish), that you receive source code or can get it if you want it, that you can change the software or use pieces of it in new free programs, and that you know you can do these things.

To protect your rights, we need to prevent others from denying you these rights or asking you to surrender the rights. Therefore, you have certain responsibilities if you distribute copies of the software, or if you modify it: responsibilities to respect the freedom of others.

For example, if you distribute copies of such a program, whether gratis or for a fee, you must pass on to the recipients the same freedoms that you received. You must make sure that they, too, receive or can get the source code. And you must show them these terms so they know their rights.

Developers that use the GNU GPL protect your rights with two steps: (1) assert copyright on the software, and (2) offer you this License giving you legal permission to copy, distribute and/or modify it.

For the developers' and authors' protection, the GPL clearly explains that there is no warranty for this free software. For both users' and authors' sake, the GPL requires that modified versions be marked as changed, so that their problems will not be attributed erroneously to authors of previous versions.

Some devices are designed to deny users access to install or run modified versions of the software inside them, although the manufacturer can do so. This is fundamentally incompatible with the aim of protecting users' freedom to change the software. The systematic pattern of such abuse occurs in the area of products for individuals to use, which is precisely where it is most unacceptable. Therefore, we have designed this version of the GPL to prohibit the practice for those products. If such problems arise substantially in other domains, we stand ready to extend this provision to those domains in future versions of the GPL, as needed to protect the freedom of users.

Finally, every program is threatened constantly by software patents. States should not allow patents to restrict development and use of software on general-purpose computers, but in those that do, we wish to avoid the special danger that patents applied to a free program could make it effectively proprietary. To prevent this, the GPL assures that patents cannot be used to render the program non-free.

The precise terms and conditions for copying, distribution and modification follow.

TERMS AND CONDITIONS

  1. Definitions.

    “This License” refers to version 3 of the GNU General Public License.

    “Copyright” also means copyright-like laws that apply to other kinds of works, such as semiconductor masks.

    “The Program” refers to any copyrightable work licensed under this License. Each licensee is addressed as “you”. “Licensees” and “recipients” may be individuals or organizations.

    To “modify” a work means to copy from or adapt all or part of the work in a fashion requiring copyright permission, other than the making of an exact copy. The resulting work is called a “modified version” of the earlier work or a work “based on” the earlier work.

    A “covered work” means either the unmodified Program or a work based on the Program.

    To “propagate” a work means to do anything with it that, without permission, would make you directly or secondarily liable for infringement under applicable copyright law, except executing it on a computer or modifying a private copy. Propagation includes copying, distribution (with or without modification), making available to the public, and in some countries other activities as well.

    To “convey” a work means any kind of propagation that enables other parties to make or receive copies. Mere interaction with a user through a computer network, with no transfer of a copy, is not conveying.

    An interactive user interface displays “Appropriate Legal Notices” to the extent that it includes a convenient and prominently visible feature that (1) displays an appropriate copyright notice, and (2) tells the user that there is no warranty for the work (except to the extent that warranties are provided), that licensees may convey the work under this License, and how to view a copy of this License. If the interface presents a list of user commands or options, such as a menu, a prominent item in the list meets this criterion.

  2. Source Code.

    The “source code” for a work means the preferred form of the work for making modifications to it. “Object code” means any non-source form of a work.

    A “Standard Interface” means an interface that either is an official standard defined by a recognized standards body, or, in the case of interfaces specified for a particular programming language, one that is widely used among developers working in that language.

    The “System Libraries” of an executable work include anything, other than the work as a whole, that (a) is included in the normal form of packaging a Major Component, but which is not part of that Major Component, and (b) serves only to enable use of the work with that Major Component, or to implement a Standard Interface for which an implementation is available to the public in source code form. A “Major Component”, in this context, means a major essential component (kernel, window system, and so on) of the specific operating system (if any) on which the executable work runs, or a compiler used to produce the work, or an object code interpreter used to run it.

    The “Corresponding Source” for a work in object code form means all the source code needed to generate, install, and (for an executable work) run the object code and to modify the work, including scripts to control those activities. However, it does not include the work's System Libraries, or general-purpose tools or generally available free programs which are used unmodified in performing those activities but which are not part of the work. For example, Corresponding Source includes interface definition files associated with source files for the work, and the source code for shared libraries and dynamically linked subprograms that the work is specifically designed to require, such as by intimate data communication or control flow between those subprograms and other parts of the work.

    The Corresponding Source need not include anything that users can regenerate automatically from other parts of the Corresponding Source.

    The Corresponding Source for a work in source code form is that same work.

  3. Basic Permissions.

    All rights granted under this License are granted for the term of copyright on the Program, and are irrevocable provided the stated conditions are met. This License explicitly affirms your unlimited permission to run the unmodified Program. The output from running a covered work is covered by this License only if the output, given its content, constitutes a covered work. This License acknowledges your rights of fair use or other equivalent, as provided by copyright law.

    You may make, run and propagate covered works that you do not convey, without conditions so long as your license otherwise remains in force. You may convey covered works to others for the sole purpose of having them make modifications exclusively for you, or provide you with facilities for running those works, provided that you comply with the terms of this License in conveying all material for which you do not control copyright. Those thus making or running the covered works for you must do so exclusively on your behalf, under your direction and control, on terms that prohibit them from making any copies of your copyrighted material outside their relationship with you.

    Conveying under any other circumstances is permitted solely under the conditions stated below. Sublicensing is not allowed; section 10 makes it unnecessary.

  4. Protecting Users' Legal Rights From Anti-Circumvention Law.

    No covered work shall be deemed part of an effective technological measure under any applicable law fulfilling obligations under article 11 of the WIPO copyright treaty adopted on 20 December 1996, or similar laws prohibiting or restricting circumvention of such measures.

    When you convey a covered work, you waive any legal power to forbid circumvention of technological measures to the extent such circumvention is effected by exercising rights under this License with respect to the covered work, and you disclaim any intention to limit operation or modification of the work as a means of enforcing, against the work's users, your or third parties' legal rights to forbid circumvention of technological measures.

  5. Conveying Verbatim Copies.

    You may convey verbatim copies of the Program's source code as you receive it, in any medium, provided that you conspicuously and appropriately publish on each copy an appropriate copyright notice; keep intact all notices stating that this License and any non-permissive terms added in accord with section 7 apply to the code; keep intact all notices of the absence of any warranty; and give all recipients a copy of this License along with the Program.

    You may charge any price or no price for each copy that you convey, and you may offer support or warranty protection for a fee.

  6. Conveying Modified Source Versions.

    You may convey a work based on the Program, or the modifications to produce it from the Program, in the form of source code under the terms of section 4, provided that you also meet all of these conditions:

    1. The work must carry prominent notices stating that you modified it, and giving a relevant date.
    2. The work must carry prominent notices stating that it is released under this License and any conditions added under section 7. This requirement modifies the requirement in section 4 to “keep intact all notices”.
    3. You must license the entire work, as a whole, under this License to anyone who comes into possession of a copy. This License will therefore apply, along with any applicable section 7 additional terms, to the whole of the work, and all its parts, regardless of how they are packaged. This License gives no permission to license the work in any other way, but it does not invalidate such permission if you have separately received it.
    4. If the work has interactive user interfaces, each must display Appropriate Legal Notices; however, if the Program has interactive interfaces that do not display Appropriate Legal Notices, your work need not make them do so.

    A compilation of a covered work with other separate and independent works, which are not by their nature extensions of the covered work, and which are not combined with it such as to form a larger program, in or on a volume of a storage or distribution medium, is called an “aggregate” if the compilation and its resulting copyright are not used to limit the access or legal rights of the compilation's users beyond what the individual works permit. Inclusion of a covered work in an aggregate does not cause this License to apply to the other parts of the aggregate.

  7. Conveying Non-Source Forms.

    You may convey a covered work in object code form under the terms of sections 4 and 5, provided that you also convey the machine-readable Corresponding Source under the terms of this License, in one of these ways:

    1. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by the Corresponding Source fixed on a durable physical medium customarily used for software interchange.
    2. Convey the object code in, or embodied in, a physical product (including a physical distribution medium), accompanied by a written offer, valid for at least three years and valid for as long as you offer spare parts or customer support for that product model, to give anyone who possesses the object code either (1) a copy of the Corresponding Source for all the software in the product that is covered by this License, on a durable physical medium customarily used for software interchange, for a price no more than your reasonable cost of physically performing this conveying of source, or (2) access to copy the Corresponding Source from a network server at no charge.
    3. Convey individual copies of the object code with a copy of the written offer to provide the Corresponding Source. This alternative is allowed only occasionally and noncommercially, and only if you received the object code with such an offer, in accord with subsection 6b.
    4. Convey the object code by offering access from a designated place (gratis or for a charge), and offer equivalent access to the Corresponding Source in the same way through the same place at no further charge. You need not require recipients to copy the Corresponding Source along with the object code. If the place to copy the object code is a network server, the Corresponding Source may be on a different server (operated by you or a third party) that supports equivalent copying facilities, provided you maintain clear directions next to the object code saying where to find the Corresponding Source. Regardless of what server hosts the Corresponding Source, you remain obligated to ensure that it is available for as long as needed to satisfy these requirements.
    5. Convey the object code using peer-to-peer transmission, provided you inform other peers where the object code and Corresponding Source of the work are being offered to the general public at no charge under subsection 6d.

    A separable portion of the object code, whose source code is excluded from the Corresponding Source as a System Library, need not be included in conveying the object code work.

    A “User Product” is either (1) a “consumer product”, which means any tangible personal property which is normally used for personal, family, or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, “normally used” refers to a typical or common use of that class of product, regardless of the status of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial commercial, industrial or non-consumer uses, unless such uses represent the only significant mode of use of the product.

    “Installation Information” for a User Product means any methods, procedures, authorization keys, or other information required to install and execute modified versions of a covered work in that User Product from a modified version of its Corresponding Source. The information must suffice to ensure that the continued functioning of the modified object code is in no case prevented or interfered with solely because modification has been made.

    If you convey an object code work under this section in, or with, or specifically for use in, a User Product, and the conveying occurs as part of a transaction in which the right of possession and use of the User Product is transferred to the recipient in perpetuity or for a fixed term (regardless of how the transaction is characterized), the Corresponding Source conveyed under this section must be accompanied by the Installation Information. But this requirement does not apply if neither you nor any third party retains the ability to install modified object code on the User Product (for example, the work has been installed in ROM).

    The requirement to provide Installation Information does not include a requirement to continue to provide support service, warranty, or updates for a work that has been modified or installed by the recipient, or for the User Product in which it has been modified or installed. Access to a network may be denied when the modification itself materially and adversely affects the operation of the network or violates the rules and protocols for communication across the network.

    Corresponding Source conveyed, and Installation Information provided, in accord with this section must be in a format that is publicly documented (and with an implementation available to the public in source code form), and must require no special password or key for unpacking, reading or copying.

  8. Additional Terms.

    “Additional permissions” are terms that supplement the terms of this License by making exceptions from one or more of its conditions. Additional permissions that are applicable to the entire Program shall be treated as though they were included in this License, to the extent that they are valid under applicable law. If additional permissions apply only to part of the Program, that part may be used separately under those permissions, but the entire Program remains governed by this License without regard to the additional permissions.

    When you convey a copy of a covered work, you may at your option remove any additional permissions from that copy, or from any part of it. (Additional permissions may be written to require their own removal in certain cases when you modify the work.) You may place additional permissions on material, added by you to a covered work, for which you have or can give appropriate copyright permission.

    Notwithstanding any other provision of this License, for material you add to a covered work, you may (if authorized by the copyright holders of that material) supplement the terms of this License with terms:

    1. Disclaiming warranty or limiting liability differently from the terms of sections 15 and 16 of this License; or
    2. Requiring preservation of specified reasonable legal notices or author attributions in that material or in the Appropriate Legal Notices displayed by works containing it; or
    3. Prohibiting misrepresentation of the origin of that material, or requiring that modified versions of such material be marked in reasonable ways as different from the original version; or
    4. Limiting the use for publicity purposes of names of licensors or authors of the material; or
    5. Declining to grant rights under trademark law for use of some trade names, trademarks, or service marks; or
    6. Requiring indemnification of licensors and authors of that material by anyone who conveys the material (or modified versions of it) with contractual assumptions of liability to the recipient, for any liability that these contractual assumptions directly impose on those licensors and authors.

    All other non-permissive additional terms are considered “further restrictions” within the meaning of section 10. If the Program as you received it, or any part of it, contains a notice stating that it is governed by this License along with a term that is a further restriction, you may remove that term. If a license document contains a further restriction but permits relicensing or conveying under this License, you may add to a covered work material governed by the terms of that license document, provided that the further restriction does not survive such relicensing or conveying.

    If you add terms to a covered work in accord with this section, you must place, in the relevant source files, a statement of the additional terms that apply to those files, or a notice indicating where to find the applicable terms.

    Additional terms, permissive or non-permissive, may be stated in the form of a separately written license, or stated as exceptions; the above requirements apply either way.

  9. Termination.

    You may not propagate or modify a covered work except as expressly provided under this License. Any attempt otherwise to propagate or modify it is void, and will automatically terminate your rights under this License (including any patent licenses granted under the third paragraph of section 11).

    However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation.

    Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice.

    Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, you do not qualify to receive new licenses for the same material under section 10.

  10. Acceptance Not Required for Having Copies.

    You are not required to accept this License in order to receive or run a copy of the Program. Ancillary propagation of a covered work occurring solely as a consequence of using peer-to-peer transmission to receive a copy likewise does not require acceptance. However, nothing other than this License grants you permission to propagate or modify any covered work. These actions infringe copyright if you do not accept this License. Therefore, by modifying or propagating a covered work, you indicate your acceptance of this License to do so.

  11. Automatic Licensing of Downstream Recipients.

    Each time you convey a covered work, the recipient automatically receives a license from the original licensors, to run, modify and propagate that work, subject to this License. You are not responsible for enforcing compliance by third parties with this License.

    An “entity transaction” is a transaction transferring control of an organization, or substantially all assets of one, or subdividing an organization, or merging organizations. If propagation of a covered work results from an entity transaction, each party to that transaction who receives a copy of the work also receives whatever licenses to the work the party's predecessor in interest had or could give under the previous paragraph, plus a right to possession of the Corresponding Source of the work from the predecessor in interest, if the predecessor has it or can get it with reasonable efforts.

    You may not impose any further restrictions on the exercise of the rights granted or affirmed under this License. For example, you may not impose a license fee, royalty, or other charge for exercise of rights granted under this License, and you may not initiate litigation (including a cross-claim or counterclaim in a lawsuit) alleging that any patent claim is infringed by making, using, selling, offering for sale, or importing the Program or any portion of it.

  12. Patents.

    A “contributor” is a copyright holder who authorizes use under this License of the Program or a work on which the Program is based. The work thus licensed is called the contributor's “contributor version”.

    A contributor's “essential patent claims” are all patent claims owned or controlled by the contributor, whether already acquired or hereafter acquired, that would be infringed by some manner, permitted by this License, of making, using, or selling its contributor version, but do not include claims that would be infringed only as a consequence of further modification of the contributor version. For purposes of this definition, “control” includes the right to grant patent sublicenses in a manner consistent with the requirements of this License.

    Each contributor grants you a non-exclusive, worldwide, royalty-free patent license under the contributor's essential patent claims, to make, use, sell, offer for sale, import and otherwise run, modify and propagate the contents of its contributor version.

    In the following three paragraphs, a “patent license” is any express agreement or commitment, however denominated, not to enforce a patent (such as an express permission to practice a patent or covenant not to sue for patent infringement). To “grant” such a patent license to a party means to make such an agreement or commitment not to enforce a patent against the party.

    If you convey a covered work, knowingly relying on a patent license, and the Corresponding Source of the work is not available for anyone to copy, free of charge and under the terms of this License, through a publicly available network server or other readily accessible means, then you must either (1) cause the Corresponding Source to be so available, or (2) arrange to deprive yourself of the benefit of the patent license for this particular work, or (3) arrange, in a manner consistent with the requirements of this License, to extend the patent license to downstream recipients. “Knowingly relying” means you have actual knowledge that, but for the patent license, your conveying the covered work in a country, or your recipient's use of the covered work in a country, would infringe one or more identifiable patents in that country that you have reason to believe are valid.

    If, pursuant to or in connection with a single transaction or arrangement, you convey, or propagate by procuring conveyance of, a covered work, and grant a patent license to some of the parties receiving the covered work authorizing them to use, propagate, modify or convey a specific copy of the covered work, then the patent license you grant is automatically extended to all recipients of the covered work and works based on it.

    A patent license is “discriminatory” if it does not include within the scope of its coverage, prohibits the exercise of, or is conditioned on the non-exercise of one or more of the rights that are specifically granted under this License. You may not convey a covered work if you are a party to an arrangement with a third party that is in the business of distributing software, under which you make payment to the third party based on the extent of your activity of conveying the work, and under which the third party grants, to any of the parties who would receive the covered work from you, a discriminatory patent license (a) in connection with copies of the covered work conveyed by you (or copies made from those copies), or (b) primarily for and in connection with specific products or compilations that contain the covered work, unless you entered into that arrangement, or that patent license was granted, prior to 28 March 2007.

    Nothing in this License shall be construed as excluding or limiting any implied license or other defenses to infringement that may otherwise be available to you under applicable patent law.

  13. No Surrender of Others' Freedom.

    If conditions are imposed on you (whether by court order, agreement or otherwise) that contradict the conditions of this License, they do not excuse you from the conditions of this License. If you cannot convey a covered work so as to satisfy simultaneously your obligations under this License and any other pertinent obligations, then as a consequence you may not convey it at all. For example, if you agree to terms that obligate you to collect a royalty for further conveying from those to whom you convey the Program, the only way you could satisfy both those terms and this License would be to refrain entirely from conveying the Program.

  14. Use with the GNU Affero General Public License.

    Notwithstanding any other provision of this License, you have permission to link or combine any covered work with a work licensed under version 3 of the GNU Affero General Public License into a single combined work, and to convey the resulting work. The terms of this License will continue to apply to the part which is the covered work, but the special requirements of the GNU Affero General Public License, section 13, concerning interaction through a network will apply to the combination as such.

  15. Revised Versions of this License.

    The Free Software Foundation may publish revised and/or new versions of the GNU General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns.

    Each version is given a distinguishing version number. If the Program specifies that a certain numbered version of the GNU General Public License “or any later version” applies to it, you have the option of following the terms and conditions either of that numbered version or of any later version published by the Free Software Foundation. If the Program does not specify a version number of the GNU General Public License, you may choose any version ever published by the Free Software Foundation.

    If the Program specifies that a proxy can decide which future versions of the GNU General Public License can be used, that proxy's public statement of acceptance of a version permanently authorizes you to choose that version for the Program.

    Later license versions may give you additional or different permissions. However, no additional obligations are imposed on any author or copyright holder as a result of your choosing to follow a later version.

  16. Disclaimer of Warranty.

    THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM “AS IS” WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.

  17. Limitation of Liability.

    IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

  18. Interpretation of Sections 15 and 16.

    If the disclaimer of warranty and limitation of liability provided above cannot be given local legal effect according to their terms, reviewing courts shall apply local law that most closely approximates an absolute waiver of all civil liability in connection with the Program, unless a warranty or assumption of liability accompanies a copy of the Program in return for a fee.

END OF TERMS AND CONDITIONS

How to Apply These Terms to Your New Programs

If you develop a new program, and you want it to be of the greatest possible use to the public, the best way to achieve this is to make it free software which everyone can redistribute and change under these terms.

To do so, attach the following notices to the program. It is safest to attach them to the start of each source file to most effectively state the exclusion of warranty; and each file should have at least the “copyright” line and a pointer to where the full notice is found.

     one line to give the program's name and a brief idea of what it does.
     Copyright (C) year name of author
     
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or (at
     your option) any later version.
     
     This program is distributed in the hope that it will be useful, but
     WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     General Public License for more details.
     
     You should have received a copy of the GNU General Public License
     along with this program.  If not, see http://www.gnu.org/licenses/.

Also add information on how to contact you by electronic and paper mail.

If the program does terminal interaction, make it output a short notice like this when it starts in an interactive mode:

     program Copyright (C) year name of author
     This program comes with ABSOLUTELY NO WARRANTY; for details type ‘show w’.
     This is free software, and you are welcome to redistribute it
     under certain conditions; type ‘show c’ for details.

The hypothetical commands ‘show w’ and ‘show c’ should show the appropriate parts of the General Public License. Of course, your program's commands might be different; for a GUI interface, you would use an “about box”.

You should also get your employer (if you work as a programmer) or school, if any, to sign a “copyright disclaimer” for the program, if necessary. For more information on this, and how to apply and follow the GNU GPL, see http://www.gnu.org/licenses/.

The GNU General Public License does not permit incorporating your program into proprietary programs. If your program is a subroutine library, you may consider it more useful to permit linking proprietary applications with the library. If this is what you want to do, use the GNU Lesser General Public License instead of this License. But first, please read http://www.gnu.org/philosophy/why-not-lgpl.html.


Next: , Previous: Package License, Up: Top

Appendix C GNU Free Documentation License

Version 1.3, 3 November 2008
     Copyright © 2000, 2001, 2002, 2007, 2008 Free Software Foundation, Inc.
     http://fsf.org/
     
     Everyone is permitted to copy and distribute verbatim copies
     of this license document, but changing it is not allowed.
  1. PREAMBLE

    The purpose of this License is to make a manual, textbook, or other functional and useful document free in the sense of freedom: to assure everyone the effective freedom to copy and redistribute it, with or without modifying it, either commercially or noncommercially. Secondarily, this License preserves for the author and publisher a way to get credit for their work, while not being considered responsible for modifications made by others.

    This License is a kind of “copyleft”, which means that derivative works of the document must themselves be free in the same sense. It complements the GNU General Public License, which is a copyleft license designed for free software.

    We have designed this License in order to use it for manuals for free software, because free software needs free documentation: a free program should come with manuals providing the same freedoms that the software does. But this License is not limited to software manuals; it can be used for any textual work, regardless of subject matter or whether it is published as a printed book. We recommend this License principally for works whose purpose is instruction or reference.

  2. APPLICABILITY AND DEFINITIONS

    This License applies to any manual or other work, in any medium, that contains a notice placed by the copyright holder saying it can be distributed under the terms of this License. Such a notice grants a world-wide, royalty-free license, unlimited in duration, to use that work under the conditions stated herein. The “Document”, below, refers to any such manual or work. Any member of the public is a licensee, and is addressed as “you”. You accept the license if you copy, modify or distribute the work in a way requiring permission under copyright law.

    A “Modified Version” of the Document means any work containing the Document or a portion of it, either copied verbatim, or with modifications and/or translated into another language.

    A “Secondary Section” is a named appendix or a front-matter section of the Document that deals exclusively with the relationship of the publishers or authors of the Document to the Document's overall subject (or to related matters) and contains nothing that could fall directly within that overall subject. (Thus, if the Document is in part a textbook of mathematics, a Secondary Section may not explain any mathematics.) The relationship could be a matter of historical connection with the subject or with related matters, or of legal, commercial, philosophical, ethical or political position regarding them.

    The “Invariant Sections” are certain Secondary Sections whose titles are designated, as being those of Invariant Sections, in the notice that says that the Document is released under this License. If a section does not fit the above definition of Secondary then it is not allowed to be designated as Invariant. The Document may contain zero Invariant Sections. If the Document does not identify any Invariant Sections then there are none.

    The “Cover Texts” are certain short passages of text that are listed, as Front-Cover Texts or Back-Cover Texts, in the notice that says that the Document is released under this License. A Front-Cover Text may be at most 5 words, and a Back-Cover Text may be at most 25 words.

    A “Transparent” copy of the Document means a machine-readable copy, represented in a format whose specification is available to the general public, that is suitable for revising the document straightforwardly with generic text editors or (for images composed of pixels) generic paint programs or (for drawings) some widely available drawing editor, and that is suitable for input to text formatters or for automatic translation to a variety of formats suitable for input to text formatters. A copy made in an otherwise Transparent file format whose markup, or absence of markup, has been arranged to thwart or discourage subsequent modification by readers is not Transparent. An image format is not Transparent if used for any substantial amount of text. A copy that is not “Transparent” is called “Opaque”.

    Examples of suitable formats for Transparent copies include plain ascii without markup, Texinfo input format, LaTeX input format, SGML or XML using a publicly available DTD, and standard-conforming simple HTML, PostScript or PDF designed for human modification. Examples of transparent image formats include PNG, XCF and JPG. Opaque formats include proprietary formats that can be read and edited only by proprietary word processors, SGML or XML for which the DTD and/or processing tools are not generally available, and the machine-generated HTML, PostScript or PDF produced by some word processors for output purposes only.

    The “Title Page” means, for a printed book, the title page itself, plus such following pages as are needed to hold, legibly, the material this License requires to appear in the title page. For works in formats which do not have any title page as such, “Title Page” means the text near the most prominent appearance of the work's title, preceding the beginning of the body of the text.

    The “publisher” means any person or entity that distributes copies of the Document to the public.

    A section “Entitled XYZ” means a named subunit of the Document whose title either is precisely XYZ or contains XYZ in parentheses following text that translates XYZ in another language. (Here XYZ stands for a specific section name mentioned below, such as “Acknowledgements”, “Dedications”, “Endorsements”, or “History”.) To “Preserve the Title” of such a section when you modify the Document means that it remains a section “Entitled XYZ” according to this definition.

    The Document may include Warranty Disclaimers next to the notice which states that this License applies to the Document. These Warranty Disclaimers are considered to be included by reference in this License, but only as regards disclaiming warranties: any other implication that these Warranty Disclaimers may have is void and has no effect on the meaning of this License.

  3. VERBATIM COPYING

    You may copy and distribute the Document in any medium, either commercially or noncommercially, provided that this License, the copyright notices, and the license notice saying this License applies to the Document are reproduced in all copies, and that you add no other conditions whatsoever to those of this License. You may not use technical measures to obstruct or control the reading or further copying of the copies you make or distribute. However, you may accept compensation in exchange for copies. If you distribute a large enough number of copies you must also follow the conditions in section 3.

    You may also lend copies, under the same conditions stated above, and you may publicly display copies.

  4. COPYING IN QUANTITY

    If you publish printed copies (or copies in media that commonly have printed covers) of the Document, numbering more than 100, and the Document's license notice requires Cover Texts, you must enclose the copies in covers that carry, clearly and legibly, all these Cover Texts: Front-Cover Texts on the front cover, and Back-Cover Texts on the back cover. Both covers must also clearly and legibly identify you as the publisher of these copies. The front cover must present the full title with all words of the title equally prominent and visible. You may add other material on the covers in addition. Copying with changes limited to the covers, as long as they preserve the title of the Document and satisfy these conditions, can be treated as verbatim copying in other respects.

    If the required texts for either cover are too voluminous to fit legibly, you should put the first ones listed (as many as fit reasonably) on the actual cover, and continue the rest onto adjacent pages.

    If you publish or distribute Opaque copies of the Document numbering more than 100, you must either include a machine-readable Transparent copy along with each Opaque copy, or state in or with each Opaque copy a computer-network location from which the general network-using public has access to download using public-standard network protocols a complete Transparent copy of the Document, free of added material. If you use the latter option, you must take reasonably prudent steps, when you begin distribution of Opaque copies in quantity, to ensure that this Transparent copy will remain thus accessible at the stated location until at least one year after the last time you distribute an Opaque copy (directly or through your agents or retailers) of that edition to the public.

    It is requested, but not required, that you contact the authors of the Document well before redistributing any large number of copies, to give them a chance to provide you with an updated version of the Document.

  5. MODIFICATIONS

    You may copy and distribute a Modified Version of the Document under the conditions of sections 2 and 3 above, provided that you release the Modified Version under precisely this License, with the Modified Version filling the role of the Document, thus licensing distribution and modification of the Modified Version to whoever possesses a copy of it. In addition, you must do these things in the Modified Version:

    1. Use in the Title Page (and on the covers, if any) a title distinct from that of the Document, and from those of previous versions (which should, if there were any, be listed in the History section of the Document). You may use the same title as a previous version if the original publisher of that version gives permission.
    2. List on the Title Page, as authors, one or more persons or entities responsible for authorship of the modifications in the Modified Version, together with at least five of the principal authors of the Document (all of its principal authors, if it has fewer than five), unless they release you from this requirement.
    3. State on the Title page the name of the publisher of the Modified Version, as the publisher.
    4. Preserve all the copyright notices of the Document.
    5. Add an appropriate copyright notice for your modifications adjacent to the other copyright notices.
    6. Include, immediately after the copyright notices, a license notice giving the public permission to use the Modified Version under the terms of this License, in the form shown in the Addendum below.
    7. Preserve in that license notice the full lists of Invariant Sections and required Cover Texts given in the Document's license notice.
    8. Include an unaltered copy of this License.
    9. Preserve the section Entitled “History”, Preserve its Title, and add to it an item stating at least the title, year, new authors, and publisher of the Modified Version as given on the Title Page. If there is no section Entitled “History” in the Document, create one stating the title, year, authors, and publisher of the Document as given on its Title Page, then add an item describing the Modified Version as stated in the previous sentence.
    10. Preserve the network location, if any, given in the Document for public access to a Transparent copy of the Document, and likewise the network locations given in the Document for previous versions it was based on. These may be placed in the “History” section. You may omit a network location for a work that was published at least four years before the Document itself, or if the original publisher of the version it refers to gives permission.
    11. For any section Entitled “Acknowledgements” or “Dedications”, Preserve the Title of the section, and preserve in the section all the substance and tone of each of the contributor acknowledgements and/or dedications given therein.
    12. Preserve all the Invariant Sections of the Document, unaltered in their text and in their titles. Section numbers or the equivalent are not considered part of the section titles.
    13. Delete any section Entitled “Endorsements”. Such a section may not be included in the Modified Version.
    14. Do not retitle any existing section to be Entitled “Endorsements” or to conflict in title with any Invariant Section.
    15. Preserve any Warranty Disclaimers.

    If the Modified Version includes new front-matter sections or appendices that qualify as Secondary Sections and contain no material copied from the Document, you may at your option designate some or all of these sections as invariant. To do this, add their titles to the list of Invariant Sections in the Modified Version's license notice. These titles must be distinct from any other section titles.

    You may add a section Entitled “Endorsements”, provided it contains nothing but endorsements of your Modified Version by various parties—for example, statements of peer review or that the text has been approved by an organization as the authoritative definition of a standard.

    You may add a passage of up to five words as a Front-Cover Text, and a passage of up to 25 words as a Back-Cover Text, to the end of the list of Cover Texts in the Modified Version. Only one passage of Front-Cover Text and one of Back-Cover Text may be added by (or through arrangements made by) any one entity. If the Document already includes a cover text for the same cover, previously added by you or by arrangement made by the same entity you are acting on behalf of, you may not add another; but you may replace the old one, on explicit permission from the previous publisher that added the old one.

    The author(s) and publisher(s) of the Document do not by this License give permission to use their names for publicity for or to assert or imply endorsement of any Modified Version.

  6. COMBINING DOCUMENTS

    You may combine the Document with other documents released under this License, under the terms defined in section 4 above for modified versions, provided that you include in the combination all of the Invariant Sections of all of the original documents, unmodified, and list them all as Invariant Sections of your combined work in its license notice, and that you preserve all their Warranty Disclaimers.

    The combined work need only contain one copy of this License, and multiple identical Invariant Sections may be replaced with a single copy. If there are multiple Invariant Sections with the same name but different contents, make the title of each such section unique by adding at the end of it, in parentheses, the name of the original author or publisher of that section if known, or else a unique number. Make the same adjustment to the section titles in the list of Invariant Sections in the license notice of the combined work.

    In the combination, you must combine any sections Entitled “History” in the various original documents, forming one section Entitled “History”; likewise combine any sections Entitled “Acknowledgements”, and any sections Entitled “Dedications”. You must delete all sections Entitled “Endorsements.”

  7. COLLECTIONS OF DOCUMENTS

    You may make a collection consisting of the Document and other documents released under this License, and replace the individual copies of this License in the various documents with a single copy that is included in the collection, provided that you follow the rules of this License for verbatim copying of each of the documents in all other respects.

    You may extract a single document from such a collection, and distribute it individually under this License, provided you insert a copy of this License into the extracted document, and follow this License in all other respects regarding verbatim copying of that document.

  8. AGGREGATION WITH INDEPENDENT WORKS

    A compilation of the Document or its derivatives with other separate and independent documents or works, in or on a volume of a storage or distribution medium, is called an “aggregate” if the copyright resulting from the compilation is not used to limit the legal rights of the compilation's users beyond what the individual works permit. When the Document is included in an aggregate, this License does not apply to the other works in the aggregate which are not themselves derivative works of the Document.

    If the Cover Text requirement of section 3 is applicable to these copies of the Document, then if the Document is less than one half of the entire aggregate, the Document's Cover Texts may be placed on covers that bracket the Document within the aggregate, or the electronic equivalent of covers if the Document is in electronic form. Otherwise they must appear on printed covers that bracket the whole aggregate.

  9. TRANSLATION

    Translation is considered a kind of modification, so you may distribute translations of the Document under the terms of section 4. Replacing Invariant Sections with translations requires special permission from their copyright holders, but you may include translations of some or all Invariant Sections in addition to the original versions of these Invariant Sections. You may include a translation of this License, and all the license notices in the Document, and any Warranty Disclaimers, provided that you also include the original English version of this License and the original versions of those notices and disclaimers. In case of a disagreement between the translation and the original version of this License or a notice or disclaimer, the original version will prevail.

    If a section in the Document is Entitled “Acknowledgements”, “Dedications”, or “History”, the requirement (section 4) to Preserve its Title (section 1) will typically require changing the actual title.

  10. TERMINATION

    You may not copy, modify, sublicense, or distribute the Document except as expressly provided under this License. Any attempt otherwise to copy, modify, sublicense, or distribute it is void, and will automatically terminate your rights under this License.

    However, if you cease all violation of this License, then your license from a particular copyright holder is reinstated (a) provisionally, unless and until the copyright holder explicitly and finally terminates your license, and (b) permanently, if the copyright holder fails to notify you of the violation by some reasonable means prior to 60 days after the cessation.

    Moreover, your license from a particular copyright holder is reinstated permanently if the copyright holder notifies you of the violation by some reasonable means, this is the first time you have received notice of violation of this License (for any work) from that copyright holder, and you cure the violation prior to 30 days after your receipt of the notice.

    Termination of your rights under this section does not terminate the licenses of parties who have received copies or rights from you under this License. If your rights have been terminated and not permanently reinstated, receipt of a copy of some or all of the same material does not give you any rights to use it.

  11. FUTURE REVISIONS OF THIS LICENSE

    The Free Software Foundation may publish new, revised versions of the GNU Free Documentation License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns. See http://www.gnu.org/copyleft/.

    Each version of the License is given a distinguishing version number. If the Document specifies that a particular numbered version of this License “or any later version” applies to it, you have the option of following the terms and conditions either of that specified version or of any later version that has been published (not as a draft) by the Free Software Foundation. If the Document does not specify a version number of this License, you may choose any version ever published (not as a draft) by the Free Software Foundation. If the Document specifies that a proxy can decide which future versions of this License can be used, that proxy's public statement of acceptance of a version permanently authorizes you to choose that version for the Document.

  12. RELICENSING

    “Massive Multiauthor Collaboration Site” (or “MMC Site”) means any World Wide Web server that publishes copyrightable works and also provides prominent facilities for anybody to edit those works. A public wiki that anybody can edit is an example of such a server. A “Massive Multiauthor Collaboration” (or “MMC”) contained in the site means any set of copyrightable works thus published on the MMC site.

    “CC-BY-SA” means the Creative Commons Attribution-Share Alike 3.0 license published by Creative Commons Corporation, a not-for-profit corporation with a principal place of business in San Francisco, California, as well as future copyleft versions of that license published by that same organization.

    “Incorporate” means to publish or republish a Document, in whole or in part, as part of another Document.

    An MMC is “eligible for relicensing” if it is licensed under this License, and if all works that were first published under this License somewhere other than this MMC, and subsequently incorporated in whole or in part into the MMC, (1) had no cover texts or invariant sections, and (2) were thus incorporated prior to November 1, 2008.

    The operator of an MMC Site may republish an MMC contained in the site under CC-BY-SA on the same site at any time before August 1, 2009, provided the MMC is eligible for relicensing.

ADDENDUM: How to use this License for your documents

To use this License in a document you have written, include a copy of the License in the document and put the following copyright and license notices just after the title page:

       Copyright (C)  year  your name.
       Permission is granted to copy, distribute and/or modify this document
       under the terms of the GNU Free Documentation License, Version 1.3
       or any later version published by the Free Software Foundation;
       with no Invariant Sections, no Front-Cover Texts, and no Back-Cover
       Texts.  A copy of the license is included in the section entitled ``GNU
       Free Documentation License''.

If you have Invariant Sections, Front-Cover Texts and Back-Cover Texts, replace the “with...Texts.” line with this:

         with the Invariant Sections being list their titles, with
         the Front-Cover Texts being list, and with the Back-Cover Texts
         being list.

If you have Invariant Sections without Cover Texts, or some other combination of the three, merge those two alternatives to suit the situation.

If your document contains nontrivial examples of program code, we recommend releasing these examples in parallel under your choice of free software license, such as the GNU General Public License, to permit their use in free software.


Next: , Previous: Documentation License, Up: Top

Appendix D Bibliography and references


Next: , Previous: references, Up: Top

Appendix E An entry for each concept


Next: , Previous: concept index, Up: Top

Appendix F An entry for each function.


Next: , Previous: function index, Up: Top

Appendix G An entry for each variable.


Previous: variable index, Up: Top

Appendix H An entry for each type.

Table of Contents