This file documents the GNU Scientific Library (GSL), a collection of numerical routines for scientific computing. It corresponds to release 1.8 of the library. Please report any errors in this manual to bug-gsl@gnu.org.
More information about GSL can be found at the project homepage, http://www.gnu.org/software/gsl/.
Printed copies of this manual can be purchased from Network Theory Ltd at http://www.network-theory.co.uk/gsl/manual/. The money raised from sales of the manual helps support the development of GSL.
Copyright © 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 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.2 or any later version published by the Free Software Foundation; with the Invariant Sections being “GNU General Public License” and “Free Software Needs Free Documentation”, the Front-Cover text being “A GNU Manual”, and with the Back-Cover Text being (a) (see below). A copy of the license is included in the section entitled “GNU Free Documentation License”.
(a) The Back-Cover Text is: “You have freedom to copy and modify this GNU Manual, like GNU software.”
The GNU Scientific Library (GSL) is a collection of routines for numerical computing. The routines have been written from scratch in C, and present a modern Applications Programming Interface (API) for C programmers, allowing wrappers to be written for very high level languages. The source code is distributed under the GNU General Public License.
The library covers a wide range of topics in numerical computing. Routines are available for the following areas,
| Complex Numbers | Roots of Polynomials
| |
| Special Functions | Vectors and Matrices
| |
| Permutations | Combinations
| |
| Sorting | BLAS Support
| |
| Linear Algebra | CBLAS Library
| |
| Fast Fourier Transforms | Eigensystems
| |
| Random Numbers | Quadrature
| |
| Random Distributions | Quasi-Random Sequences
| |
| Histograms | Statistics
| |
| Monte Carlo Integration | N-Tuples
| |
| Differential Equations | Simulated Annealing
| |
| Numerical Differentiation | Interpolation
| |
| Series Acceleration | Chebyshev Approximations
| |
| Root-Finding | Discrete Hankel Transforms
| |
| Least-Squares Fitting | Minimization
| |
| IEEE Floating-Point | Physical Constants
| |
| Wavelets
|
The use of these routines is described in this manual. Each chapter provides detailed definitions of the functions, followed by example programs and references to the articles on which the algorithms are based.
Where possible the routines have been based on reliable public-domain packages such as FFTPACK and QUADPACK, which the developers of GSL have reimplemented in C with modern coding conventions.
The subroutines in the GNU Scientific Library are “free software”; this means that everyone is free to use them, and to redistribute them in other free programs. The library is not in the public domain; it is copyrighted and there are conditions on its distribution. These conditions are designed to permit everything that a good cooperating citizen would want to do. What is not allowed is to try to prevent others from further sharing any version of the software that they might get from you.
Specifically, we want to make sure that you have the right to share copies of programs that you are given which use the GNU Scientific Library, that you receive their source code or else can get it if you want it, that you can change these programs or use pieces of them in new free programs, and that you know you can do these things.
To make sure that everyone has such rights, we have to forbid you to deprive anyone else of these rights. For example, if you distribute copies of any code which uses the GNU Scientific Library, you must give the recipients all the rights that you have received. You must make sure that they, too, receive or can get the source code, both to the library and the code which uses it. And you must tell them their rights. This means that the library should not be redistributed in proprietary programs.
Also, for our own protection, we must make certain that everyone finds out that there is no warranty for the GNU Scientific Library. If these programs are modified by someone else and passed on, we want their recipients to know that what they have is not what we distributed, so that any problems introduced by others will not reflect on our reputation.
The precise conditions for the distribution of software related to the GNU Scientific Library are found in the GNU General Public License (see GNU General Public License). Further information about this license is available from the GNU Project webpage Frequently Asked Questions about the GNU GPL,
The Free Software Foundation also operates a license consulting service for commercial users (contact details available from http://www.fsf.org/).
The source code for the library can be obtained in different ways, by copying it from a friend, purchasing it on cdrom or downloading it from the internet. A list of public ftp servers which carry the source code can be found on the GNU website,
The preferred platform for the library is a GNU system, which allows it to take advantage of additional features in the GNU C compiler and GNU C library. However, the library is fully portable and should compile on most systems with a C compiler. Precompiled versions of the library can be purchased from commercial redistributors listed on the website above.
Announcements of new releases, updates and other relevant events are
made on the info-gsl@gnu.org mailing list. To subscribe to this
low-volume list, send an email of the following form:
To: info-gsl-request@gnu.org
Subject: subscribe
You will receive a response asking you to reply in order to confirm your subscription.
The software described in this manual has no warranty, it is provided “as is”. It is your responsibility to validate the behavior of the routines and their accuracy using the source code provided, or to purchase support and warranties from commercial redistributors. Consult the GNU General Public license for further details (see GNU General Public License).
A list of known bugs can be found in the BUGS file included in the GSL distribution. Details of compilation problems can be found in the INSTALL file.
If you find a bug which is not listed in these files, please report it to bug-gsl@gnu.org.
All bug reports should include:
It is useful if you can check whether the same problem occurs when the library is compiled without optimization. Thank you.
Any errors or omissions in this manual can also be reported to the same address.
Additional information, including online copies of this manual, links to related projects, and mailing list archives are available from the website mentioned above.
Any questions about the use and installation of the library can be asked
on the mailing list help-gsl@gnu.org. To subscribe to this
list, send an email of the following form:
To: help-gsl-request@gnu.org
Subject: subscribe
This mailing list can be used to ask questions not covered by this manual, and to contact the developers of the library.
If you would like to refer to the GNU Scientific Library in a journal article, the recommended way is to cite this reference manual, e.g. M. Galassi et al, GNU Scientific Library Reference Manual (2nd Ed.), ISBN 0954161734.
If you want to give a url, use “http://www.gnu.org/software/gsl/”.
This manual contains many examples which can be typed at the keyboard. A command entered at the terminal is shown like this,
$ command
The first character on the line is the terminal prompt, and should not be typed. The dollar sign `$' is used as the standard prompt in this manual, although some systems may use a different character.
The examples assume the use of the GNU operating system. There may be
minor differences in the output on other systems. The commands for
setting environment variables use the Bourne shell syntax of the
standard GNU shell (bash).
This chapter describes how to compile programs that use GSL, and introduces its conventions.
The following short program demonstrates the use of the library by computing the value of the Bessel function J_0(x) for x=5,
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
int
main (void)
{
double x = 5.0;
double y = gsl_sf_bessel_J0 (x);
printf ("J0(%g) = %.18e\n", x, y);
return 0;
}
The output is shown below, and should be correct to double-precision accuracy,
J0(5) = -1.775967713143382920e-01
The steps needed to compile this program are described in the following sections.
The library header files are installed in their own gsl directory. You should write any preprocessor include statements with a gsl/ directory prefix thus,
#include <gsl/gsl_math.h>
If the directory is not installed on the standard search path of your
compiler you will also need to provide its location to the preprocessor
as a command line flag. The default location of the gsl
directory is /usr/local/include/gsl. A typical compilation
command for a source file example.c with the GNU C compiler
gcc is,
$ gcc -Wall -I/usr/local/include -c example.c
This results in an object file example.o. The default
include path for gcc searches /usr/local/include automatically so
the -I option can actually be omitted when GSL is installed
in its default location.
The library is installed as a single file, libgsl.a. A shared version of the library libgsl.so is also installed on systems that support shared libraries. The default location of these files is /usr/local/lib. If this directory is not on the standard search path of your linker you will also need to provide its location as a command line flag.
To link against the library you need to specify both the main library and a supporting cblas library, which provides standard basic linear algebra subroutines. A suitable cblas implementation is provided in the library libgslcblas.a if your system does not provide one. The following example shows how to link an application with the library,
$ gcc -L/usr/local/lib example.o -lgsl -lgslcblas -lm
The default library path for gcc searches /usr/local/lib
automatically so the -L option can be omitted when GSL is
installed in its default location.
The following command line shows how you would link the same application with an alternative cblas library called libcblas,
$ gcc example.o -lgsl -lcblas -lm
For the best performance an optimized platform-specific cblas
library should be used for -lcblas. The library must conform to
the cblas standard. The atlas package provides a portable
high-performance blas library with a cblas interface. It is
free software and should be installed for any work requiring fast vector
and matrix operations. The following command line will link with the
atlas library and its cblas interface,
$ gcc example.o -lgsl -lcblas -latlas -lm
For more information see BLAS Support.
To run a program linked with the shared version of the library the operating system must be able to locate the corresponding .so file at runtime. If the library cannot be found, the following error will occur:
$ ./a.out
./a.out: error while loading shared libraries:
libgsl.so.0: cannot open shared object file: No such
file or directory
To avoid this error, define the shell variable LD_LIBRARY_PATH to
include the directory where the library is installed.
For example, in the Bourne shell (/bin/sh or /bin/bash),
the library search path can be set with the following commands:
$ LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
$ export LD_LIBRARY_PATH
$ ./example
In the C-shell (/bin/csh or /bin/tcsh) the equivalent
command is,
% setenv LD_LIBRARY_PATH /usr/local/lib:$LD_LIBRARY_PATH
The standard prompt for the C-shell in the example above is the percent character `%', and should not be typed as part of the command.
To save retyping these commands each session they should be placed in an individual or system-wide login file.
To compile a statically linked version of the program, use the
-static flag in gcc,
$ gcc -static example.o -lgsl -lgslcblas -lm
The library is written in ANSI C and is intended to conform to the ANSI C standard (C89). It should be portable to any system with a working ANSI C compiler.
The library does not rely on any non-ANSI extensions in the interface it exports to the user. Programs you write using GSL can be ANSI compliant. Extensions which can be used in a way compatible with pure ANSI C are supported, however, via conditional compilation. This allows the library to take advantage of compiler extensions on those platforms which support them.
When an ANSI C feature is known to be broken on a particular system the library will exclude any related functions at compile-time. This should make it impossible to link a program that would use these functions and give incorrect results.
To avoid namespace conflicts all exported function names and variables
have the prefix gsl_, while exported macros have the prefix
GSL_.
The inline keyword is not part of the original ANSI C standard
(C89) and the library does not export any inline function definitions by
default. However, the library provides optional inline versions of
performance-critical functions by conditional compilation. The inline
versions of these functions can be included by defining the macro
HAVE_INLINE when compiling an application,
$ gcc -Wall -c -DHAVE_INLINE example.c
If you use autoconf this macro can be defined automatically. If
you do not define the macro HAVE_INLINE then the slower
non-inlined versions of the functions will be used instead.
Note that the actual usage of the inline keyword is extern
inline, which eliminates unnecessary function definitions in gcc.
If the form extern inline causes problems with other compilers a
stricter autoconf test can be used, see Autoconf Macros.
The extended numerical type long double is part of the ANSI C
standard and should be available in every modern compiler. However, the
precision of long double is platform dependent, and this should
be considered when using it. The IEEE standard only specifies the
minimum precision of extended precision numbers, while the precision of
double is the same on all platforms.
In some system libraries the stdio.h formatted input/output
functions printf and scanf are not implemented correctly
for long double. Undefined or incorrect results are avoided by
testing these functions during the configure stage of library
compilation and eliminating certain GSL functions which depend on them
if necessary. The corresponding line in the configure output
looks like this,
checking whether printf works with long double... no
Consequently when long double formatted input/output does not
work on a given system it should be impossible to link a program which
uses GSL functions dependent on this.
If it is necessary to work on a system which does not support formatted
long double input/output then the options are to use binary
formats or to convert long double results into double for
reading and writing.
To help in writing portable applications GSL provides some implementations of functions that are found in other libraries, such as the BSD math library. You can write your application to use the native versions of these functions, and substitute the GSL versions via a preprocessor macro if they are unavailable on another platform.
For example, after determining whether the BSD function hypot is
available you can include the following macro definitions in a file
config.h with your application,
/* Substitute gsl_hypot for missing system hypot */
#ifndef HAVE_HYPOT
#define hypot gsl_hypot
#endif
The application source files can then use the include command
#include <config.h> to replace each occurrence of hypot by
gsl_hypot when hypot is not available. This substitution
can be made automatically if you use autoconf, see Autoconf Macros.
In most circumstances the best strategy is to use the native versions of these functions when available, and fall back to GSL versions otherwise, since this allows your application to take advantage of any platform-specific optimizations in the system library. This is the strategy used within GSL itself.
The main implementation of some functions in the library will not be optimal on all architectures. For example, there are several ways to compute a Gaussian random variate and their relative speeds are platform-dependent. In cases like this the library provides alternative implementations of these functions with the same interface. If you write your application using calls to the standard implementation you can select an alternative version later via a preprocessor definition. It is also possible to introduce your own optimized functions this way while retaining portability. The following lines demonstrate the use of a platform-dependent choice of methods for sampling from the Gaussian distribution,
#ifdef SPARC
#define gsl_ran_gaussian gsl_ran_gaussian_ratio_method
#endif
#ifdef INTEL
#define gsl_ran_gaussian my_gaussian
#endif
These lines would be placed in the configuration header file config.h of the application, which should then be included by all the source files. Note that the alternative implementations will not produce bit-for-bit identical results, and in the case of random number distributions will produce an entirely different stream of random variates.
Many functions in the library are defined for different numeric types.
This feature is implemented by varying the name of the function with a
type-related modifier—a primitive form of C++ templates. The
modifier is inserted into the function name after the initial module
prefix. The following table shows the function names defined for all
the numeric types of an imaginary module gsl_foo with function
fn,
gsl_foo_fn double
gsl_foo_long_double_fn long double
gsl_foo_float_fn float
gsl_foo_long_fn long
gsl_foo_ulong_fn unsigned long
gsl_foo_int_fn int
gsl_foo_uint_fn unsigned int
gsl_foo_short_fn short
gsl_foo_ushort_fn unsigned short
gsl_foo_char_fn char
gsl_foo_uchar_fn unsigned char
The normal numeric precision double is considered the default and
does not require a suffix. For example, the function
gsl_stats_mean computes the mean of double precision numbers,
while the function gsl_stats_int_mean computes the mean of
integers.
A corresponding scheme is used for library defined types, such as
gsl_vector and gsl_matrix. In this case the modifier is
appended to the type name. For example, if a module defines a new
type-dependent struct or typedef gsl_foo it is modified for other
types in the following way,
gsl_foo double
gsl_foo_long_double long double
gsl_foo_float float
gsl_foo_long long
gsl_foo_ulong unsigned long
gsl_foo_int int
gsl_foo_uint unsigned int
gsl_foo_short short
gsl_foo_ushort unsigned short
gsl_foo_char char
gsl_foo_uchar unsigned char
When a module contains type-dependent definitions the library provides individual header files for each type. The filenames are modified as shown in the below. For convenience the default header includes the definitions for all the types. To include only the double precision header file, or any other specific type, use its individual filename.
#include <gsl/gsl_foo.h> All types
#include <gsl/gsl_foo_double.h> double
#include <gsl/gsl_foo_long_double.h> long double
#include <gsl/gsl_foo_float.h> float
#include <gsl/gsl_foo_long.h> long
#include <gsl/gsl_foo_ulong.h> unsigned long
#include <gsl/gsl_foo_int.h> int
#include <gsl/gsl_foo_uint.h> unsigned int
#include <gsl/gsl_foo_short.h> short
#include <gsl/gsl_foo_ushort.h> unsigned short
#include <gsl/gsl_foo_char.h> char
#include <gsl/gsl_foo_uchar.h> unsigned char
The library header files automatically define functions to have
extern "C" linkage when included in C++ programs. This allows
the functions to be called directly from C++.
To use C++ exception handling within user-defined functions passed to
the library as parameters, the library must be built with the
additional CFLAGS compilation option -fexceptions.
The library assumes that arrays, vectors and matrices passed as
modifiable arguments are not aliased and do not overlap with each other.
This removes the need for the library to handle overlapping memory
regions as a special case, and allows additional optimizations to be
used. If overlapping memory regions are passed as modifiable arguments
then the results of such functions will be undefined. If the arguments
will not be modified (for example, if a function prototype declares them
as const arguments) then overlapping or aliased memory regions
can be safely used.
The library can be used in multi-threaded programs. All the functions
are thread-safe, in the sense that they do not use static variables.
Memory is always associated with objects and not with functions. For
functions which use workspace objects as temporary storage the
workspaces should be allocated on a per-thread basis. For functions
which use table objects as read-only memory the tables can be used
by multiple threads simultaneously. Table arguments are always declared
const in function prototypes, to indicate that they may be
safely accessed by different threads.
There are a small number of static global variables which are used to control the overall behavior of the library (e.g. whether to use range-checking, the function to call on fatal error, etc). These variables are set directly by the user, so they should be initialized once at program startup and not modified by different threads.
From time to time, it may be necessary for the definitions of some
functions to be altered or removed from the library. In these
circumstances the functions will first be declared deprecated and
then removed from subsequent versions of the library. Functions that
are deprecated can be disabled in the current release by setting the
preprocessor definition GSL_DISABLE_DEPRECATED. This allows
existing code to be tested for forwards compatibility.
Where possible the routines in the library have been written to avoid
dependencies between modules and files. This should make it possible to
extract individual functions for use in your own applications, without
needing to have the whole library installed. You may need to define
certain macros such as GSL_ERROR and remove some #include
statements in order to compile the files as standalone units. Reuse of
the library code in this way is encouraged, subject to the terms of the
GNU General Public License.
This chapter describes the way that GSL functions report and handle errors. By examining the status information returned by every function you can determine whether it succeeded or failed, and if it failed you can find out what the precise cause of failure was. You can also define your own error handling functions to modify the default behavior of the library.
The functions described in this section are declared in the header file gsl_errno.h.
The library follows the thread-safe error reporting conventions of the
posix Threads library. Functions return a non-zero error code to
indicate an error and 0 to indicate success.
int status = gsl_function (...)
if (status) { /* an error occurred */
.....
/* status value specifies the type of error */
}
The routines report an error whenever they cannot perform the task requested of them. For example, a root-finding function would return a non-zero error code if could not converge to the requested accuracy, or exceeded a limit on the number of iterations. Situations like this are a normal occurrence when using any mathematical library and you should check the return status of the functions that you call.
Whenever a routine reports an error the return value specifies the type
of error. The return value is analogous to the value of the variable
errno in the C library. The caller can examine the return code
and decide what action to take, including ignoring the error if it is
not considered serious.
In addition to reporting errors by return codes the library also has an
error handler function gsl_error. This function is called by
other library functions when they report an error, just before they
return to the caller. The default behavior of the error handler is to
print a message and abort the program,
gsl: file.c:67: ERROR: invalid argument supplied by user
Default GSL error handler invoked.
Aborted
The purpose of the gsl_error handler is to provide a function
where a breakpoint can be set that will catch library errors when
running under the debugger. It is not intended for use in production
programs, which should handle any errors using the return codes.
The error code numbers returned by library functions are defined in the
file gsl_errno.h. They all have the prefix GSL_ and
expand to non-zero constant integer values. Many of the error codes use
the same base name as the corresponding error code in the C library. Here are
some of the most common error codes,
Domain error; used by mathematical functions when an argument value does not fall into the domain over which the function is defined (like EDOM in the C library)
Range error; used by mathematical functions when the result value is not representable because of overflow or underflow (like ERANGE in the C library)
No memory available. The system cannot allocate more virtual memory because its capacity is full (like ENOMEM in the C library). This error is reported when a GSL routine encounters problems when trying to allocate memory with
malloc.
Invalid argument. This is used to indicate various kinds of problems with passing the wrong argument to a library function (like EINVAL in the C library).
The error codes can be converted into an error message using the
function gsl_strerror.
This function returns a pointer to a string describing the error code gsl_errno. For example,
printf ("error: %s\n", gsl_strerror (status));would print an error message like
error: output range errorfor a status value ofGSL_ERANGE.
The default behavior of the GSL error handler is to print a short
message and call abort(). When this default is in use programs
will stop with a core-dump whenever a library routine reports an error.
This is intended as a fail-safe default for programs which do not check
the return status of library routines (we don't encourage you to write
programs this way).
If you turn off the default error handler it is your responsibility to check the return values of routines and handle them yourself. You can also customize the error behavior by providing a new error handler. For example, an alternative error handler could log all errors to a file, ignore certain error conditions (such as underflows), or start the debugger and attach it to the current process when an error occurs.
All GSL error handlers have the type gsl_error_handler_t, which is
defined in gsl_errno.h,
This is the type of GSL error handler functions. An error handler will be passed four arguments which specify the reason for the error (a string), the name of the source file in which it occurred (also a string), the line number in that file (an integer) and the error number (an integer). The source file and line number are set at compile time using the
__FILE__and__LINE__directives in the preprocessor. An error handler function returns typevoid. Error handler functions should be defined like this,void handler (const char * reason, const char * file, int line, int gsl_errno)
To request the use of your own error handler you need to call the
function gsl_set_error_handler which is also declared in
gsl_errno.h,
This function sets a new error handler, new_handler, for the GSL library routines. The previous handler is returned (so that you can restore it later). Note that the pointer to a user defined error handler function is stored in a static variable, so there can be only one error handler per program. This function should be not be used in multi-threaded programs except to set up a program-wide error handler from a master thread. The following example shows how to set and restore a new error handler,
/* save original handler, install new handler */ old_handler = gsl_set_error_handler (&my_handler); /* code uses new handler */ ..... /* restore original handler */ gsl_set_error_handler (old_handler);To use the default behavior (
aborton error) set the error handler toNULL,old_handler = gsl_set_error_handler (NULL);
This function turns off the error handler by defining an error handler which does nothing. This will cause the program to continue after any error, so the return values from any library routines must be checked. This is the recommended behavior for production programs. The previous handler is returned (so that you can restore it later).
The error behavior can be changed for specific applications by
recompiling the library with a customized definition of the
GSL_ERROR macro in the file gsl_errno.h.
If you are writing numerical functions in a program which also uses GSL code you may find it convenient to adopt the same error reporting conventions as in the library.
To report an error you need to call the function gsl_error with a
string describing the error and then return an appropriate error code
from gsl_errno.h, or a special value, such as NaN. For
convenience the file gsl_errno.h defines two macros which carry
out these steps:
This macro reports an error using the GSL conventions and returns a status value of
gsl_errno. It expands to the following code fragment,gsl_error (reason, __FILE__, __LINE__, gsl_errno); return gsl_errno;The macro definition in gsl_errno.h actually wraps the code in a
do { ... } while (0)block to prevent possible parsing problems.
Here is an example of how the macro could be used to report that a
routine did not achieve a requested tolerance. To report the error the
routine needs to return the error code GSL_ETOL.
if (residual > tolerance)
{
GSL_ERROR("residual exceeds tolerance", GSL_ETOL);
}
This macro is the same as
GSL_ERRORbut returns a user-defined value of value instead of an error code. It can be used for mathematical functions that return a floating point value.
The following example shows how to return a NaN at a mathematical
singularity using the GSL_ERROR_VAL macro,
if (x == 0)
{
GSL_ERROR_VAL("argument lies on singularity",
GSL_ERANGE, GSL_NAN);
}
Here is an example of some code which checks the return value of a function where an error might be reported,
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
...
int status;
size_t n = 37;
gsl_set_error_handler_off();
status = gsl_fft_complex_radix2_forward (data, n);
if (status) {
if (status == GSL_EINVAL) {
fprintf (stderr, "invalid argument, n=%d\n", n);
} else {
fprintf (stderr, "failed, gsl_errno=%d\n",
status);
}
exit (-1);
}
...
The function gsl_fft_complex_radix2 only accepts integer lengths
which are a power of two. If the variable n is not a power of
two then the call to the library function will return GSL_EINVAL,
indicating that the length argument is invalid. The function call to
gsl_set_error_handler_off() stops the default error handler from
aborting the program. The else clause catches any other possible
errors.
This chapter describes basic mathematical functions. Some of these functions are present in system libraries, but the alternative versions given here can be used as a substitute when the system functions are not available.
The functions and macros described in this chapter are defined in the header file gsl_math.h.
The library ensures that the standard bsd mathematical constants are defined. For reference, here is a list of the constants:
M_EM_LOG2EM_LOG10EM_SQRT2M_SQRT1_2M_SQRT3M_PIM_PI_2M_PI_4M_SQRTPIM_2_SQRTPIM_1_PIM_2_PIM_LN10M_LN2M_LNPIM_EULERThis macro contains the IEEE representation of positive infinity, +\infty. It is computed from the expression
+1.0/0.0.
This macro contains the IEEE representation of negative infinity, -\infty. It is computed from the expression
-1.0/0.0.
This macro contains the IEEE representation of the Not-a-Number symbol,
NaN. It is computed from the ratio0.0/0.0.
This function returns +1 if x is positive infinity, -1 if x is negative infinity and 0 otherwise.
This function returns 1 if x is a real number, and 0 if it is infinite or not-a-number.
The following routines provide portable implementations of functions
found in the BSD math library. When native versions are not available
the functions described here can be used instead. The substitution can
be made automatically if you use autoconf to compile your
application (see Portability functions).
This function computes the value of \log(1+x) in a way that is accurate for small x. It provides an alternative to the BSD math function
log1p(x).
This function computes the value of \exp(x)-1 in a way that is accurate for small x. It provides an alternative to the BSD math function
expm1(x).
This function computes the value of \sqrt{x^2 + y^2} in a way that avoids overflow. It provides an alternative to the BSD math function
hypot(x,y).
This function computes the value of \arccosh(x). It provides an alternative to the standard math function
acosh(x).
This function computes the value of \arcsinh(x). It provides an alternative to the standard math function
asinh(x).
This function computes the value of \arctanh(x). It provides an alternative to the standard math function
atanh(x).
This function computes the value of x * 2^e. It provides an alternative to the standard math function
ldexp(x,e).
This function splits the number x into its normalized fraction f and exponent e, such that x = f * 2^e and 0.5 <= f < 1. The function returns f and stores the exponent in e. If x is zero, both f and e are set to zero. This function provides an alternative to the standard math function
frexp(x, e).
A common complaint about the standard C library is its lack of a function for calculating (small) integer powers. GSL provides a simple functions to fill this gap. For reasons of efficiency, these functions do not check for overflow or underflow conditions.
This routine computes the power x^n for integer n. The power is computed efficiently—for example, x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications. A version of this function which also computes the numerical error in the result is available as
gsl_sf_pow_int_e.
These functions can be used to compute small integer powers x^2, x^3, etc. efficiently. The functions will be inlined when possible so that use of these functions should be as efficient as explicitly writing the corresponding product expression.
#include <gsl/gsl_math.h>
double y = gsl_pow_4 (3.141) /* compute 3.141**4 */
This macro returns the sign of x. It is defined as
((x) >= 0 ? 1 : -1). Note that with this definition the sign of zero is positive (regardless of its ieee sign bit).
This macro evaluates to 1 if n is odd and 0 if n is even. The argument n must be of integer type.
This macro is the opposite of
GSL_IS_ODD(n). It evaluates to 1 if n is even and 0 if n is odd. The argument n must be of integer type.
This macro returns the maximum of a and b. It is defined as
((a) > (b) ? (a):(b)).
This macro returns the minimum of a and b. It is defined as
((a) < (b) ? (a):(b)).
This function returns the maximum of the double precision numbers a and b using an inline function. The use of a function allows for type checking of the arguments as an extra safety feature. On platforms where inline functions are not available the macro
GSL_MAXwill be automatically substituted.
This function returns the minimum of the double precision numbers a and b using an inline function. The use of a function allows for type checking of the arguments as an extra safety feature. On platforms where inline functions are not available the macro
GSL_MINwill be automatically substituted.
These functions return the maximum or minimum of the integers a and b using an inline function. On platforms where inline functions are not available the macros
GSL_MAXorGSL_MINwill be automatically substituted.
These functions return the maximum or minimum of the long doubles a and b using an inline function. On platforms where inline functions are not available the macros
GSL_MAXorGSL_MINwill be automatically substituted.
It is sometimes useful to be able to compare two floating point numbers approximately, to allow for rounding and truncation errors. The following function implements the approximate floating-point comparison algorithm proposed by D.E. Knuth in Section 4.2.2 of Seminumerical Algorithms (3rd edition).
This function determines whether x and y are approximately equal to a relative accuracy epsilon.
The relative accuracy is measured using an interval of size 2 \delta, where \delta = 2^k \epsilon and k is the maximum base-2 exponent of x and y as computed by the function
frexp().If x and y lie within this interval, they are considered approximately equal and the function returns 0. Otherwise if x < y, the function returns -1, or if x > y, the function returns +1.
The implementation is based on the package
fcmpby T.C. Belding.
The functions described in this chapter provide support for complex numbers. The algorithms take care to avoid unnecessary intermediate underflows and overflows, allowing the functions to be evaluated over as much of the complex plane as possible.
For multiple-valued functions the branch cuts have been chosen to follow the conventions of Abramowitz and Stegun in the Handbook of Mathematical Functions. The functions return principal values which are the same as those in GNU Calc, which in turn are the same as those in Common Lisp, The Language (Second Edition)1 and the HP-28/48 series of calculators.
The complex types are defined in the header file gsl_complex.h, while the corresponding complex functions and arithmetic operations are defined in gsl_complex_math.h.
Complex numbers are represented using the type gsl_complex. The
internal representation of this type may vary across platforms and
should not be accessed directly. The functions and macros described
below allow complex numbers to be manipulated in a portable way.
For reference, the default form of the gsl_complex type is
given by the following struct,
typedef struct
{
double dat[2];
} gsl_complex;
The real and imaginary part are stored in contiguous elements of a two
element array. This eliminates any padding between the real and
imaginary parts, dat[0] and dat[1], allowing the struct to
be mapped correctly onto packed complex arrays.
This function uses the rectangular cartesian components (x,y) to return the complex number z = x + i y.
This function returns the complex number z = r \exp(i \theta) = r (\cos(\theta) + i \sin(\theta)) from the polar representation (r,theta).
These macros return the real and imaginary parts of the complex number z.
This macro uses the cartesian components (x,y) to set the real and imaginary parts of the complex number pointed to by zp. For example,
GSL_SET_COMPLEX(&z, 3, 4)sets z to be 3 + 4i.
These macros allow the real and imaginary parts of the complex number pointed to by zp to be set independently.
This function returns the argument of the complex number z, \arg(z), where -\pi < \arg(z) <= \pi.
This function returns the magnitude of the complex number z, |z|.
This function returns the squared magnitude of the complex number z, |z|^2.
This function returns the natural logarithm of the magnitude of the complex number z, \log|z|. It allows an accurate evaluation of \log|z| when |z| is close to one. The direct evaluation of
log(gsl_complex_abs(z))would lead to a loss of precision in this case.
This function returns the sum of the complex numbers a and b, z=a+b.
This function returns the difference of the complex numbers a and b, z=a-b.
This function returns the product of the complex numbers a and b, z=ab.
This function returns the quotient of the complex numbers a and b, z=a/b.
This function returns the sum of the complex number a and the real number x, z=a+x.
This function returns the difference of the complex number a and the real number x, z=a-x.
This function returns the product of the complex number a and the real number x, z=ax.
This function returns the quotient of the complex number a and the real number x, z=a/x.
This function returns the sum of the complex number a and the imaginary number iy, z=a+iy.
This function returns the difference of the complex number a and the imaginary number iy, z=a-iy.
This function returns the product of the complex number a and the imaginary number iy, z=a*(iy).
This function returns the quotient of the complex number a and the imaginary number iy, z=a/(iy).
This function returns the complex conjugate of the complex number z, z^* = x - i y.
This function returns the inverse, or reciprocal, of the complex number z, 1/z = (x - i y)/(x^2 + y^2).
This function returns the negative of the complex number z, -z = (-x) + i(-y).
This function returns the square root of the complex number z, \sqrt z. The branch cut is the negative real axis. The result always lies in the right half of the complex plane.
This function returns the complex square root of the real number x, where x may be negative.
The function returns the complex number z raised to the complex power a, z^a. This is computed as \exp(\log(z)*a) using complex logarithms and complex exponentials.
This function returns the complex number z raised to the real power x, z^x.
This function returns the complex exponential of the complex number z, \exp(z).
This function returns the complex natural logarithm (base e) of the complex number z, \log(z). The branch cut is the negative real axis.
This function returns the complex base-10 logarithm of the complex number z, \log_10 (z).
This function returns the complex base-b logarithm of the complex number z, \log_b(z). This quantity is computed as the ratio \log(z)/\log(b).
This function returns the complex sine of the complex number z, \sin(z) = (\exp(iz) - \exp(-iz))/(2i).
This function returns the complex cosine of the complex number z, \cos(z) = (\exp(iz) + \exp(-iz))/2.
This function returns the complex tangent of the complex number z, \tan(z) = \sin(z)/\cos(z).
This function returns the complex secant of the complex number z, \sec(z) = 1/\cos(z).
This function returns the complex cosecant of the complex number z, \csc(z) = 1/\sin(z).
This function returns the complex cotangent of the complex number z, \cot(z) = 1/\tan(z).
This function returns the complex arcsine of the complex number z, \arcsin(z). The branch cuts are on the real axis, less than -1 and greater than 1.
This function returns the complex arcsine of the real number z, \arcsin(z). For z between -1 and 1, the function returns a real value in the range [-\pi/2,\pi/2]. For z less than -1 the result has a real part of -\pi/2 and a positive imaginary part. For z greater than 1 the result has a real part of \pi/2 and a negative imaginary part.
This function returns the complex arccosine of the complex number z, \arccos(z). The branch cuts are on the real axis, less than -1 and greater than 1.
This function returns the complex arccosine of the real number z, \arccos(z). For z between -1 and 1, the function returns a real value in the range [0,\pi]. For z less than -1 the result has a real part of \pi and a negative imaginary part. For z greater than 1 the result is purely imaginary and positive.
This function returns the complex arctangent of the complex number z, \arctan(z). The branch cuts are on the imaginary axis, below -i and above i.
This function returns the complex arcsecant of the complex number z, \arcsec(z) = \arccos(1/z).
This function returns the complex arcsecant of the real number z, \arcsec(z) = \arccos(1/z).
This function returns the complex arccosecant of the complex number z, \arccsc(z) = \arcsin(1/z).
This function returns the complex arccosecant of the real number z, \arccsc(z) = \arcsin(1/z).
This function returns the complex arccotangent of the complex number z, \arccot(z) = \arctan(1/z).
This function returns the complex hyperbolic sine of the complex number z, \sinh(z) = (\exp(z) - \exp(-z))/2.
This function returns the complex hyperbolic cosine of the complex number z, \cosh(z) = (\exp(z) + \exp(-z))/2.
This function returns the complex hyperbolic tangent of the complex number z, \tanh(z) = \sinh(z)/\cosh(z).
This function returns the complex hyperbolic secant of the complex number z, \sech(z) = 1/\cosh(z).
This function returns the complex hyperbolic cosecant of the complex number z, \csch(z) = 1/\sinh(z).
This function returns the complex hyperbolic cotangent of the complex number z, \coth(z) = 1/\tanh(z).
This function returns the complex hyperbolic arcsine of the complex number z, \arcsinh(z). The branch cuts are on the imaginary axis, below -i and above i.
This function returns the complex hyperbolic arccosine of the complex number z, \arccosh(z). The branch cut is on the real axis, less than 1.
This function returns the complex hyperbolic arccosine of the real number z, \arccosh(z).
This function returns the complex hyperbolic arctangent of the complex number z, \arctanh(z). The branch cuts are on the real axis, less than -1 and greater than 1.
This function returns the complex hyperbolic arctangent of the real number z, \arctanh(z).
This function returns the complex hyperbolic arcsecant of the complex number z, \arcsech(z) = \arccosh(1/z).
This function returns the complex hyperbolic arccosecant of the complex number z, \arccsch(z) = \arcsin(1/z).
This function returns the complex hyperbolic arccotangent of the complex number z, \arccoth(z) = \arctanh(1/z).
The implementations of the elementary and trigonometric functions are based on the following papers,
The general formulas and details of branch cuts can be found in the following books,
This chapter describes functions for evaluating and solving polynomials.
There are routines for finding real and complex roots of quadratic and
cubic equations using analytic methods. An iterative polynomial solver
is also available for finding the roots of general polynomials with real
coefficients (of any order). The functions are declared in the header
file gsl_poly.h.
This function evaluates the polynomial c[0] + c[1] x + c[2] x^2 + \dots + c[len-1] x^{len-1} using Horner's method for stability. The function is inlined when possible.
The functions described here manipulate polynomials stored in Newton's divided-difference representation. The use of divided-differences is described in Abramowitz & Stegun sections 25.1.4 and 25.2.26.
This function computes a divided-difference representation of the interpolating polynomial for the points (xa, ya) stored in the arrays xa and ya of length size. On output the divided-differences of (xa,ya) are stored in the array dd, also of length size.
This function evaluates the polynomial stored in divided-difference form in the arrays dd and xa of length size at the point x.
This function converts the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation is supplied in the arrays dd and xa of length size. On output the Taylor coefficients of the polynomial expanded about the point xp are stored in the array c also of length size. A workspace of length size must be provided in the array w.
This function finds the real roots of the quadratic equation,
a x^2 + b x + c = 0The number of real roots (either zero, one or two) is returned, and their locations are stored in x0 and x1. If no real roots are found then x0 and x1 are not modified. If one real root is found (i.e. if a=0) then it is stored in x0. When two real roots are found they are stored in x0 and x1 in ascending order. The case of coincident roots is not considered special. For example (x-1)^2=0 will have two roots, which happen to have exactly equal values.
The number of roots found depends on the sign of the discriminant b^2 - 4 a c. This will be subject to rounding and cancellation errors when computed in double precision, and will also be subject to errors if the coefficients of the polynomial are inexact. These errors may cause a discrete change in the number of roots. However, for polynomials with small integer coefficients the discriminant can always be computed exactly.
This function finds the complex roots of the quadratic equation,
a z^2 + b z + c = 0The number of complex roots is returned (either one or two) and the locations of the roots are stored in z0 and z1. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components. If only one real root is found (i.e. if a=0) then it is stored in z0.
This function finds the real roots of the cubic equation,
x^3 + a x^2 + b x + c = 0with a leading coefficient of unity. The number of real roots (either one or three) is returned, and their locations are stored in x0, x1 and x2. If one real root is found then only x0 is modified. When three real roots are found they are stored in x0, x1 and x2 in ascending order. The case of coincident roots is not considered special. For example, the equation (x-1)^3=0 will have three roots with exactly equal values.
This function finds the complex roots of the cubic equation,
z^3 + a z^2 + b z + c = 0The number of complex roots is returned (always three) and the locations of the roots are stored in z0, z1 and z2. The roots are returned in ascending order, sorted first by their real components and then by their imaginary components.
The roots of polynomial equations cannot be found analytically beyond the special cases of the quadratic, cubic and quartic equation. The algorithm described in this section uses an iterative method to find the approximate locations of roots of higher order polynomials.
This function allocates space for a
gsl_poly_complex_workspacestruct and a workspace suitable for solving a polynomial with n coefficients using the routinegsl_poly_complex_solve.The function returns a pointer to the newly allocated
gsl_poly_complex_workspaceif no errors were detected, and a null pointer in the case of error.
This function frees all the memory associated with the workspace w.
This function computes the roots of the general polynomial P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using balanced-QR reduction of the companion matrix. The parameter n specifies the length of the coefficient array. The coefficient of the highest order term must be non-zero. The function requires a workspace w of the appropriate size. The n-1 roots are returned in the packed complex array z of length 2(n-1), alternating real and imaginary parts.
The function returns
GSL_SUCCESSif all the roots are found andGSL_EFAILEDif the QR reduction does not converge. Note that due to finite precision, roots of higher multiplicity are returned as a cluster of simple roots with reduced accuracy. The solution of polynomials with higher-order roots requires specialized algorithms that take the multiplicity structure into account (see e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical Software, Volume 30, Issue 2 (2004), pp 218–236).
To demonstrate the use of the general polynomial solver we will take the polynomial P(x) = x^5 - 1 which has the following roots,
1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}
The following program will find these roots.
#include <stdio.h>
#include <gsl/gsl_poly.h>
int
main (void)
{
int i;
/* coefficients of P(x) = -1 + x^5 */
double a[6] = { -1, 0, 0, 0, 0, 1 };
double z[10];
gsl_poly_complex_workspace * w
= gsl_poly_complex_workspace_alloc (6);
gsl_poly_complex_solve (a, 6, w, z);
gsl_poly_complex_workspace_free (w);
for (i = 0; i < 5; i++)
{
printf ("z%d = %+.18f %+.18f\n",
i, z[2*i], z[2*i+1]);
}
return 0;
}
The output of the program is,
$ ./a.out
z0 = -0.809016994374947451 +0.587785252292473137
z1 = -0.809016994374947451 -0.587785252292473137
z2 = +0.309016994374947451 +0.951056516295153642
z3 = +0.309016994374947451 -0.951056516295153642
z4 = +1.000000000000000000 +0.000000000000000000
which agrees with the analytic result, z_n = \exp(2 \pi n i/5).
The balanced-QR method and its error analysis are described in the following papers,
The formulas for divided differences are given in Abramowitz and Stegun,
This chapter describes the GSL special function library. The library includes routines for calculating the values of Airy functions, Bessel functions, Clausen functions, Coulomb wave functions, Coupling coefficients, the Dawson function, Debye functions, Dilogarithms, Elliptic integrals, Jacobi elliptic functions, Error functions, Exponential integrals, Fermi-Dirac functions, Gamma functions, Gegenbauer functions, Hypergeometric functions, Laguerre functions, Legendre functions and Spherical Harmonics, the Psi (Digamma) Function, Synchrotron functions, Transport functions, Trigonometric functions and Zeta functions. Each routine also computes an estimate of the numerical error in the calculated value of the function.
The functions in this chapter are declared in individual header files, such as gsl_sf_airy.h, gsl_sf_bessel.h, etc. The complete set of header files can be included using the file gsl_sf.h.
The special functions are available in two calling conventions, a natural form which returns the numerical value of the function and an error-handling form which returns an error code. The two types of function provide alternative ways of accessing the same underlying code.
The natural form returns only the value of the function and can be used directly in mathematical expressions. For example, the following function call will compute the value of the Bessel function J_0(x),
double y = gsl_sf_bessel_J0 (x);
There is no way to access an error code or to estimate the error using this method. To allow access to this information the alternative error-handling form stores the value and error in a modifiable argument,
gsl_sf_result result;
int status = gsl_sf_bessel_J0_e (x, &result);
The error-handling functions have the suffix _e. The returned
status value indicates error conditions such as overflow, underflow or
loss of precision. If there are no errors the error-handling functions
return GSL_SUCCESS.
The error handling form of the special functions always calculate an error estimate along with the value of the result. Therefore, structures are provided for amalgamating a value and error estimate. These structures are declared in the header file gsl_sf_result.h.
The gsl_sf_result struct contains value and error fields.
typedef struct
{
double val;
double err;
} gsl_sf_result;
The field val contains the value and the field err contains an estimate of the absolute error in the value.
In some cases, an overflow or underflow can be detected and handled by a
function. In this case, it may be possible to return a scaling exponent
as well as an error/value pair in order to save the result from
exceeding the dynamic range of the built-in types. The
gsl_sf_result_e10 struct contains value and error fields as well
as an exponent field such that the actual result is obtained as
result * 10^(e10).
typedef struct
{
double val;
double err;
int e10;
} gsl_sf_result_e10;
The goal of the library is to achieve double precision accuracy wherever
possible. However the cost of evaluating some special functions to
double precision can be significant, particularly where very high order
terms are required. In these cases a mode argument allows the
accuracy of the function to be reduced in order to improve performance.
The following precision levels are available for the mode argument,
GSL_PREC_DOUBLEGSL_PREC_SINGLEGSL_PREC_APPROXThe approximate mode provides the fastest evaluation at the lowest accuracy.
The Airy functions Ai(x) and Bi(x) are defined by the integral representations,
Ai(x) = (1/\pi) \int_0^\infty \cos((1/3) t^3 + xt) dt
Bi(x) = (1/\pi) \int_0^\infty (e^(-(1/3) t^3) + \sin((1/3) t^3 + xt)) dt
For further information see Abramowitz & Stegun, Section 10.4. The Airy functions are defined in the header file gsl_sf_airy.h.
These routines compute the Airy function Ai(x) with an accuracy specified by mode.
These routines compute the Airy function Bi(x) with an accuracy specified by mode.
These routines compute a scaled version of the Airy function S_A(x) Ai(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3) x^(3/2)), and is 1 for x<0.
These routines compute a scaled version of the Airy function S_B(x) Bi(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)), and is 1 for x<0.
These routines compute the Airy function derivative Ai'(x) with an accuracy specified by mode.
These routines compute the Airy function derivative Bi'(x) with an accuracy specified by mode.
These routines compute the scaled Airy function derivative S_A(x) Ai'(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3) x^(3/2)), and is 1 for x<0.
These routines compute the scaled Airy function derivative S_B(x) Bi'(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)), and is 1 for x<0.
These routines compute the location of the s-th zero of the Airy function Ai(x).
These routines compute the location of the s-th zero of the Airy function Bi(x).
These routines compute the location of the s-th zero of the Airy function derivative Ai'(x).
These routines compute the location of the s-th zero of the Airy function derivative Bi'(x).
The routines described in this section compute the Cylindrical Bessel functions J_n(x), Y_n(x), Modified cylindrical Bessel functions I_n(x), K_n(x), Spherical Bessel functions j_l(x), y_l(x), and Modified Spherical Bessel functions i_l(x), k_l(x). For more information see Abramowitz & Stegun, Chapters 9 and 10. The Bessel functions are defined in the header file gsl_sf_bessel.h.
These routines compute the regular cylindrical Bessel function of zeroth order, J_0(x).
These routines compute the regular cylindrical Bessel function of first order, J_1(x).
These routines compute the regular cylindrical Bessel function of order n, J_n(x).
This routine computes the values of the regular cylindrical Bessel functions J_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the irregular cylindrical Bessel function of zeroth order, Y_0(x), for x>0.
These routines compute the irregular cylindrical Bessel function of first order, Y_1(x), for x>0.
These routines compute the irregular cylindrical Bessel function of order n, Y_n(x), for x>0.
This routine computes the values of the irregular cylindrical Bessel functions Y_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The domain of the function is x>0. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the regular modified cylindrical Bessel function of zeroth order, I_0(x).
These routines compute the regular modified cylindrical Bessel function of first order, I_1(x).
These routines compute the regular modified cylindrical Bessel function of order n, I_n(x).
This routine computes the values of the regular modified cylindrical Bessel functions I_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The start of the range nmin must be positive or zero. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the scaled regular modified cylindrical Bessel function of zeroth order \exp(-|x|) I_0(x).
These routines compute the scaled regular modified cylindrical Bessel function of first order \exp(-|x|) I_1(x).
These routines compute the scaled regular modified cylindrical Bessel function of order n, \exp(-|x|) I_n(x)
This routine computes the values of the scaled regular cylindrical Bessel functions \exp(-|x|) I_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The start of the range nmin must be positive or zero. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the irregular modified cylindrical Bessel function of zeroth order, K_0(x), for x > 0.
These routines compute the irregular modified cylindrical Bessel function of first order, K_1(x), for x > 0.
These routines compute the irregular modified cylindrical Bessel function of order n, K_n(x), for x > 0.
This routine computes the values of the irregular modified cylindrical Bessel functions K_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The start of the range nmin must be positive or zero. The domain of the function is x>0. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the scaled irregular modified cylindrical Bessel function of zeroth order \exp(x) K_0(x) for x>0.
These routines compute the scaled irregular modified cylindrical Bessel function of first order \exp(x) K_1(x) for x>0.
These routines compute the scaled irregular modified cylindrical Bessel function of order n, \exp(x) K_n(x), for x>0.
This routine computes the values of the scaled irregular cylindrical Bessel functions \exp(x) K_n(x) for n from nmin to nmax inclusive, storing the results in the array result_array. The start of the range nmin must be positive or zero. The domain of the function is x>0. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the regular spherical Bessel function of zeroth order, j_0(x) = \sin(x)/x.
These routines compute the regular spherical Bessel function of first order, j_1(x) = (\sin(x)/x - \cos(x))/x.
These routines compute the regular spherical Bessel function of second order, j_2(x) = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x.
These routines compute the regular spherical Bessel function of order l, j_l(x), for l >= 0 and x >= 0.
This routine computes the values of the regular spherical Bessel functions j_l(x) for l from 0 to lmax inclusive for lmax >= 0 and x >= 0, storing the results in the array result_array. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
This routine uses Steed's method to compute the values of the regular spherical Bessel functions j_l(x) for l from 0 to lmax inclusive for lmax >= 0 and x >= 0, storing the results in the array result_array. The Steed/Barnett algorithm is described in Comp. Phys. Comm. 21, 297 (1981). Steed's method is more stable than the recurrence used in the other functions but is also slower.
These routines compute the irregular spherical Bessel function of zeroth order, y_0(x) = -\cos(x)/x.
These routines compute the irregular spherical Bessel function of first order, y_1(x) = -(\cos(x)/x + \sin(x))/x.
These routines compute the irregular spherical Bessel function of second order, y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x).
These routines compute the irregular spherical Bessel function of order l, y_l(x), for l >= 0.
This routine computes the values of the irregular spherical Bessel functions y_l(x) for l from 0 to lmax inclusive for lmax >= 0, storing the results in the array result_array. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
The regular modified spherical Bessel functions i_l(x) are related to the modified Bessel functions of fractional order, i_l(x) = \sqrt{\pi/(2x)} I_{l+1/2}(x)
These routines compute the scaled regular modified spherical Bessel function of zeroth order, \exp(-|x|) i_0(x).
These routines compute the scaled regular modified spherical Bessel function of first order, \exp(-|x|) i_1(x).
These routines compute the scaled regular modified spherical Bessel function of second order, \exp(-|x|) i_2(x)
These routines compute the scaled regular modified spherical Bessel function of order l, \exp(-|x|) i_l(x)
This routine computes the values of the scaled regular modified cylindrical Bessel functions \exp(-|x|) i_l(x) for l from 0 to lmax inclusive for lmax >= 0, storing the results in the array result_array. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
The irregular modified spherical Bessel functions k_l(x) are related to the irregular modified Bessel functions of fractional order, k_l(x) = \sqrt{\pi/(2x)} K_{l+1/2}(x).
These routines compute the scaled irregular modified spherical Bessel function of zeroth order, \exp(x) k_0(x), for x>0.
These routines compute the scaled irregular modified spherical Bessel function of first order, \exp(x) k_1(x), for x>0.
These routines compute the scaled irregular modified spherical Bessel function of second order, \exp(x) k_2(x), for x>0.
These routines compute the scaled irregular modified spherical Bessel function of order l, \exp(x) k_l(x), for x>0.
This routine computes the values of the scaled irregular modified spherical Bessel functions \exp(x) k_l(x) for l from 0 to lmax inclusive for lmax >= 0 and x>0, storing the results in the array result_array. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values.
These routines compute the regular cylindrical Bessel function of fractional order \nu, J_\nu(x).
This function computes the regular cylindrical Bessel function of fractional order \nu, J_\nu(x), evaluated at a series of x values. The array v of length size contains the x values. They are assumed to be strictly ordered and positive. The array is over-written with the values of J_\nu(x_i).
These routines compute the irregular cylindrical Bessel function of fractional order \nu, Y_\nu(x).
These routines compute the regular modified Bessel function of fractional order \nu, I_\nu(x) for x>0, \nu>0.
These routines compute the scaled regular modified Bessel function of fractional order \nu, \exp(-|x|)I_\nu(x) for x>0, \nu>0.
These routines compute the irregular modified Bessel function of fractional order \nu, K_\nu(x) for x>0, \nu>0.
These routines compute the logarithm of the irregular modified Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0, \nu>0.
These routines compute the scaled irregular modified Bessel function of fractional order \nu, \exp(+|x|) K_\nu(x) for x>0, \nu>0.
These routines compute the location of the s-th positive zero of the Bessel function J_0(x).
These routines compute the location of the s-th positive zero of the Bessel function J_1(x).
These routines compute the location of the s-th positive zero of the Bessel function J_\nu(x). The current implementation does not support negative values of nu.
The Clausen function is defined by the following integral,
Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2))
It is related to the dilogarithm by Cl_2(\theta) = \Im Li_2(\exp(i\theta)). The Clausen functions are declared in the header file gsl_sf_clausen.h.
These routines compute the Clausen integral Cl_2(x).
The prototypes of the Coulomb functions are declared in the header file gsl_sf_coulomb.h. Both bound state and scattering solutions are available.
These routines compute the lowest-order normalized hydrogenic bound state radial wavefunction R_1 := 2Z \sqrt{Z} \exp(-Z r).
These routines compute the n-th normalized hydrogenic bound state radial wavefunction,
R_n := 2 (Z^{3/2}/n^2) \sqrt{(n-l-1)!/(n+l)!} \exp(-Z r/n) (2Zr/n)^l L^{2l+1}_{n-l-1}(2Zr/n).where L^a_b(x) is the generalized Laguerre polynomial (see Laguerre Functions). The normalization is chosen such that the wavefunction \psi is given by \psi(n,l,r) = R_n Y_{lm}.
The Coulomb wave functions F_L(\eta,x), G_L(\eta,x) are
described in Abramowitz & Stegun, Chapter 14. Because there can be a
large dynamic range of values for these functions, overflows are handled
gracefully. If an overflow occurs, GSL_EOVRFLW is signalled and
exponent(s) are returned through the modifiable parameters exp_F,
exp_G. The full solution can be reconstructed from the following
relations,
F_L(eta,x) = fc[k_L] * exp(exp_F)
G_L(eta,x) = gc[k_L] * exp(exp_G)
F_L'(eta,x) = fcp[k_L] * exp(exp_F)
G_L'(eta,x) = gcp[k_L] * exp(exp_G)
This function computes the Coulomb wave functions F_L(\eta,x), G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x), G'_{L-k}(\eta,x) with respect to x. The parameters are restricted to L, L-k > -1/2, x > 0 and integer k. Note that L itself is not restricted to being an integer. The results are stored in the parameters F, G for the function values and Fp, Gp for the derivative values. If an overflow occurs,
GSL_EOVRFLWis returned and scaling exponents are stored in the modifiable parameters exp_F, exp_G.
This function computes the Coulomb wave function F_L(\eta,x) for L = Lmin \dots Lmin + kmax, storing the results in fc_array. In the case of overflow the exponent is stored in F_exponent.
This function computes the functions F_L(\eta,x), G_L(\eta,x) for L = Lmin \dots Lmin + kmax storing the results in fc_array and gc_array. In the case of overflow the exponents are stored in F_exponent and G_exponent.
This function computes the functions F_L(\eta,x), G_L(\eta,x) and their derivatives F'_L(\eta,x), G'_L(\eta,x) for L = Lmin \dots Lmin + kmax storing the results in fc_array, gc_array, fcp_array and gcp_array. In the case of overflow the exponents are stored in F_exponent and G_exponent.
This function computes the Coulomb wave function divided by the argument F_L(\eta, x)/x for L = Lmin \dots Lmin + kmax, storing the results in fc_array. In the case of overflow the exponent is stored in F_exponent. This function reduces to spherical Bessel functions in the limit \eta \to 0.
The Coulomb wave function normalization constant is defined in Abramowitz 14.1.7.
This function computes the Coulomb wave function normalization constant C_L(\eta) for L > -1.
This function computes the Coulomb wave function normalization constant C_L(\eta) for L = Lmin \dots Lmin + kmax, Lmin > -1.
The Wigner 3-j, 6-j and 9-j symbols give the coupling coefficients for combined angular momentum vectors. Since the arguments of the standard coupling coefficient functions are integer or half-integer, the arguments of the following functions are, by convention, integers equal to twice the actual spin value. For information on the 3-j coefficients see Abramowitz & Stegun, Section 27.9. The functions described in this section are declared in the header file gsl_sf_coupling.h.
These routines compute the Wigner 3-j coefficient,
(ja jb jc ma mb mc)where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.
These routines compute the Wigner 6-j coefficient,
{ja jb jc jd je jf}where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.
These routines compute the Wigner 9-j coefficient,
{ja jb jc jd je jf jg jh ji}where the arguments are given in half-integer units, ja = two_ja/2, ma = two_ma/2, etc.
The Dawson integral is defined by \exp(-x^2) \int_0^x dt \exp(t^2). A table of Dawson's integral can be found in Abramowitz & Stegun, Table 7.5. The Dawson functions are declared in the header file gsl_sf_dawson.h.
These routines compute the value of Dawson's integral for x.
The Debye functions D_n(x) are defined by the following integral,
D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1))
For further information see Abramowitz & Stegun, Section 27.1. The Debye functions are declared in the header file gsl_sf_debye.h.
These routines compute the first-order Debye function D_1(x) = (1/x) \int_0^x dt (t/(e^t - 1)).
These routines compute the second-order Debye function D_2(x) = (2/x^2) \int_0^x dt (t^2/(e^t - 1)).
These routines compute the third-order Debye function D_3(x) = (3/x^3) \int_0^x dt (t^3/(e^t - 1)).
These routines compute the fourth-order Debye function D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1)).
These routines compute the fifth-order Debye function D_5(x) = (5/x^5) \int_0^x dt (t^5/(e^t - 1)).
These routines compute the fourth-order Debye function D_6(x) = (6/x^6) \int_0^x dt (t^6/(e^t - 1)).
The functions described in this section are declared in the header file gsl_sf_dilog.h.
These routines compute the dilogarithm for a real argument. In Lewin's notation this is Li_2(x), the real part of the dilogarithm of a real x. It is defined by the integral representation Li_2(x) = - \Re \int_0^x ds \log(1-s) / s. Note that \Im(Li_2(x)) = 0 for x <= 1, and -\pi\log(x) for x > 1.
This function computes the full complex-valued dilogarithm for the complex argument z = r \exp(i \theta). The real and imaginary parts of the result are returned in result_re, result_im.
The following functions allow for the propagation of errors when combining quantities by multiplication. The functions are declared in the header file gsl_sf_elementary.h.
This function multiplies x and y storing the product and its associated error in result.
This function multiplies x and y with associated absolute errors dx and dy. The product xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2) is stored in result.
The functions described in this section are declared in the header file gsl_sf_ellint.h.
The Legendre forms of elliptic integrals F(\phi,k), E(\phi,k) and P(\phi,k,n) are defined by,
F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t)))
E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t)))
P(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))
The complete Legendre forms are denoted by K(k) = F(\pi/2, k) and E(k) = E(\pi/2, k). Further information on the Legendre forms of elliptic integrals can be found in Abramowitz & Stegun, Chapter 17. The notation used here is based on Carlson, Numerische Mathematik 33 (1979) 1 and differs slightly from that used by Abramowitz & Stegun.
The Carlson symmetric forms of elliptical integrals RC(x,y), RD(x,y,z), RF(x,y,z) and RJ(x,y,z,p) are defined by,
RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1)
RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2)
RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2)
RJ(x,y,z,p) = 3/2 \int_0^\infty dt
(t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1)
These routines compute the complete elliptic integral K(k) to the accuracy specified by the mode variable mode.
These routines compute the complete elliptic integral E(k) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral F(\phi,k) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral E(\phi,k) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral P(\phi,k,n) to the accuracy specified by the mode variable mode.
These functions compute the incomplete elliptic integral D(\phi,k,n) which is defined through the Carlson form RD(x,y,z) by the following relation,
D(\phi,k,n) = RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).
These routines compute the incomplete elliptic integral RC(x,y) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral RD(x,y,z) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral RF(x,y,z) to the accuracy specified by the mode variable mode.
These routines compute the incomplete elliptic integral RJ(x,y,z,p) to the accuracy specified by the mode variable mode.
The Jacobian Elliptic functions are defined in Abramowitz & Stegun, Chapter 16. The functions are declared in the header file gsl_sf_elljac.h.
This function computes the Jacobian elliptic functions sn(u|m), cn(u|m), dn(u|m) by descending Landen transformations.
The error function is described in Abramowitz & Stegun, Chapter 7. The functions in this section are declared in the header file gsl_sf_erf.h.
These routines compute the error function erf(x), where erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2).
These routines compute the complementary error function erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2).
These routines compute the logarithm of the complementary error function \log(\erfc(x)).
The probability functions for the Normal or Gaussian distribution are described in Abramowitz & Stegun, Section 26.2.
These routines compute the Gaussian probability density function Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2).
These routines compute the upper tail of the Gaussian probability function Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2).
The hazard function for the normal distribution, also known as the inverse Mill's ratio, is defined as,
h(x) = Z(x)/Q(x) = \sqrt{2/\pi} \exp(-x^2 / 2) / \erfc(x/\sqrt 2)
It decreases rapidly as x approaches -\infty and asymptotes to h(x) \sim x as x approaches +\infty.
These routines compute the hazard function for the normal distribution.
The functions described in this section are declared in the header file gsl_sf_exp.h.
These routines provide an exponential function \exp(x) using GSL semantics and error checking.
This function computes the exponential \exp(x) using the
gsl_sf_result_e10type to return a result with extended range. This function may be useful if the value of \exp(x) would overflow the numeric range ofdouble.
These routines exponentiate x and multiply by the factor y to return the product y \exp(x).
This function computes the product y \exp(x) using the
gsl_sf_result_e10type to return a result with extended numeric range.
These routines compute the quantity \exp(x)-1 using an algorithm that is accurate for small x.
These routines compute the quantity (\exp(x)-1)/x using an algorithm that is accurate for small x. For small x the algorithm is based on the expansion (\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + \dots.
These routines compute the quantity 2(\exp(x)-1-x)/x^2 using an algorithm that is accurate for small x. For small x the algorithm is based on the expansion 2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + \dots.
These routines compute the N-relative exponential, which is the n-th generalization of the functions
gsl_sf_exprelandgsl_sf_exprel2. The N-relative exponential is given by,exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!) = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ... = 1F1 (1,1+N,x)
This function exponentiates x with an associated absolute error dx.
This function exponentiates a quantity x with an associated absolute error dx using the
gsl_sf_result_e10type to return a result with extended range.
This routine computes the product y \exp(x) for the quantities x, y with associated absolute errors dx, dy.
This routine computes the product y \exp(x) for the quantities x, y with associated absolute errors dx, dy using the
gsl_sf_result_e10type to return a result with extended range.
Information on the exponential integrals can be found in Abramowitz & Stegun, Chapter 5. These functions are declared in the header file gsl_sf_expint.h.
These routines compute the exponential integral E_1(x),
E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t.
These routines compute the second-order exponential integral E_2(x),
E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.
These routines compute the exponential integral Ei(x),
Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t)where PV denotes the principal value of the integral.
These routines compute the integral Shi(x) = \int_0^x dt \sinh(t)/t.
These routines compute the integral Chi(x) := \Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t] , where \gamma_E is the Euler constant (available as the macro
M_EULER).
These routines compute the third-order exponential integral Ei_3(x) = \int_0^xdt \exp(-t^3) for x >= 0.
These routines compute the Sine integral Si(x) = \int_0^x dt \sin(t)/t.
These routines compute the Cosine integral Ci(x) = -\int_x^\infty dt \cos(t)/t for x > 0.
These routines compute the Arctangent integral, which is defined as AtanInt(x) = \int_0^x dt \arctan(t)/t.
The functions described in this section are declared in the header file gsl_sf_fermi_dirac.h.
The complete Fermi-Dirac integral F_j(x) is given by,
F_j(x) := (1/r\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1))
These routines compute the complete Fermi-Dirac integral with an index of -1. This integral is given by F_{-1}(x) = e^x / (1 + e^x).
These routines compute the complete Fermi-Dirac integral with an index of 0. This integral is given by F_0(x) = \ln(1 + e^x).
These routines compute the complete Fermi-Dirac integral with an index of 1, F_1(x) = \int_0^\infty dt (t /(\exp(t-x)+1)).
These routines compute the complete Fermi-Dirac integral with an index of 2, F_2(x) = (1/2) \int_0^\infty dt (t^2 /(\exp(t-x)+1)).
These routines compute the complete Fermi-Dirac integral with an integer index of j, F_j(x) = (1/\Gamma(j+1)) \int_0^\infty dt (t^j /(\exp(t-x)+1)).
These routines compute the complete Fermi-Dirac integral F_{-1/2}(x).
These routines compute the complete Fermi-Dirac integral F_{1/2}(x).
These routines compute the complete Fermi-Dirac integral F_{3/2}(x).
The incomplete Fermi-Dirac integral F_j(x,b) is given by,
F_j(x,b) := (1/\Gamma(j+1)) \int_b^\infty dt (t^j / (\Exp(t-x) + 1))
These routines compute the incomplete Fermi-Dirac integral with an index of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x).
The functions described in this section are declared in the header file gsl_sf_gamma.h.
The Gamma function is defined by the following integral,
\Gamma(x) = \int_0^\infty dt t^{x-1} \exp(-t)
It is related to the factorial function by \Gamma(n)=(n-1)! for positive integer n. Further information on the Gamma function can be found in Abramowitz & Stegun, Chapter 6. The functions described in this section are declared in the header file gsl_sf_gamma.h.
These routines compute the Gamma function \Gamma(x), subject to x not being a negative integer. The function is computed using the real Lanczos method. The maximum value of x such that \Gamma(x) is not considered an overflow is given by the macro
GSL_SF_GAMMA_XMAXand is 171.0.
These routines compute the logarithm of the Gamma function, \log(\Gamma(x)), subject to x not a being negative integer. For x<0 the real part of \log(\Gamma(x)) is returned, which is equivalent to \log(|\Gamma(x)|). The function is computed using the real Lanczos method.
This routine computes the sign of the gamma function and the logarithm of its magnitude, subject to x not being a negative integer. The function is computed using the real Lanczos method. The value of the gamma function can be reconstructed using the relation \Gamma(x) = sgn * \exp(resultlg).
These routines compute the regulated Gamma Function \Gamma^*(x) for x > 0. The regulated gamma function is given by,
\Gamma^*(x) = \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x)) = (1 + (1/12x) + ...) for x \to \inftyand is a useful suggestion of Temme.
These routines compute the reciprocal of the gamma function, 1/\Gamma(x) using the real Lanczos method.
This routine computes \log(\Gamma(z)) for complex z=z_r+i z_i and z not a negative integer, using the complex Lanczos method. The returned parameters are lnr = \log|\Gamma(z)| and arg = \arg(\Gamma(z)) in (-\pi,\pi]. Note that the phase part (arg) is not well-determined when |z| is very large, due to inevitable roundoff in restricting to (-\pi,\pi]. This will result in a
GSL_ELOSSerror when it occurs. The absolute value part (lnr), however, never suffers from loss of precision.
Although factorials can be computed from the Gamma function, using the relation n! = \Gamma(n+1) for non-negative integer n, it is usually more efficient to call the functions in this section, particularly for small values of n, whose factorial values are maintained in hardcoded tables.
These routines compute the factorial n!. The factorial is related to the Gamma function by n! = \Gamma(n+1). The maximum value of n such that n! is not considered an overflow is given by the macro
GSL_SF_FACT_NMAXand is 170.
These routines compute the double factorial n!! = n(n-2)(n-4) \dots. The maximum value of n such that n!! is not considered an overflow is given by the macro
GSL_SF_DOUBLEFACT_NMAXand is 297.
These routines compute the logarithm of the factorial of n, \log(n!). The algorithm is faster than computing \ln(\Gamma(n+1)) via
gsl_sf_lngammafor n < 170, but defers for larger n.
These routines compute the logarithm of the double factorial of n, \log(n!!).
These routines compute the combinatorial factor
n choose m= n!/(m!(n-m)!)
These routines compute the logarithm of
n choose m. This is equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).
These routines compute the Taylor coefficient x^n / n! for x >= 0, n >= 0.
These routines compute the Pochhammer symbol (a)_x = \Gamma(a + x)/\Gamma(a), subject to a and a+x not being negative integers. The Pochhammer symbol is also known as the Apell symbol and sometimes written as (a,x).
These routines compute the logarithm of the Pochhammer symbol, \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)) for a > 0, a+x > 0.
These routines compute the sign of the Pochhammer symbol and the logarithm of its magnitude. The computed parameters are result = \log(|(a)_x|) and sgn = \sgn((a)_x) where (a)_x = \Gamma(a + x)/\Gamma(a), subject to a, a+x not being negative integers.
These routines compute the relative Pochhammer symbol ((a)_x - 1)/x where (a)_x = \Gamma(a + x)/\Gamma(a).
These functions compute the unnormalized incomplete Gamma Function \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x >= 0.
These routines compute the normalized incomplete Gamma Function Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0, x >= 0.
These routines compute the complementary normalized incomplete Gamma Function P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt t^{a-1} \exp(-t) for a > 0, x >= 0.
Note that Abramowitz & Stegun call P(a,x) the incomplete gamma function (section 6.5).
These routines compute the Beta Function, B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) for a > 0, b > 0.
These routines compute the logarithm of the Beta Function, \log(B(a,b)) for a > 0, b > 0.
These routines compute the normalized incomplete Beta function B_x(a,b)/B(a,b) where B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt for a > 0, b > 0, and 0 <= x <= 1.
The Gegenbauer polynomials are defined in Abramowitz & Stegun, Chapter 22, where they are known as Ultraspherical polynomials. The functions described in this section are declared in the header file gsl_sf_gegenbauer.h.
These functions evaluate the Gegenbauer polynomials C^{(\lambda)}_n(x) using explicit representations for n =1, 2, 3.
These functions evaluate the Gegenbauer polynomial C^{(\lambda)}_n(x) for a specific value of n, lambda, x subject to \lambda > -1/2, n >= 0.
This function computes an array of Gegenbauer polynomials C^{(\lambda)}_n(x) for n = 0, 1, 2, \dots, nmax, subject to \lambda > -1/2, nmax >= 0.
Hypergeometric functions are described in Abramowitz & Stegun, Chapters 13 and 15. These functions are declared in the header file gsl_sf_hyperg.h.
These routines compute the hypergeometric function 0F1(c,x).
These routines compute the confluent hypergeometric function 1F1(m,n,x) = M(m,n,x) for integer parameters m, n.
These routines compute the confluent hypergeometric function 1F1(a,b,x) = M(a,b,x) for general parameters a, b.
These routines compute the confluent hypergeometric function U(m,n,x) for integer parameters m, n.
This routine computes the confluent hypergeometric function U(m,n,x) for integer parameters m, n using the
gsl_sf_result_e10type to return a result with extended range.
These routines compute the confluent hypergeometric function U(a,b,x).
This routine computes the confluent hypergeometric function U(a,b,x) using the
gsl_sf_result_e10type to return a result with extended range.
These routines compute the Gauss hypergeometric function 2F1(a,b,c,x) for |x| < 1.
If the arguments (a,b,c,x) are too close to a singularity then the function can return the error code
GSL_EMAXITERwhen the series approximation converges too slowly. This occurs in the region of x=1, c - a - b = m for integer m.
These routines compute the Gauss hypergeometric function 2F1(a_R + i a_I, a_R - i a_I, c, x) with complex parameters for |x| < 1. exceptions:
These routines compute the renormalized Gauss hypergeometric function 2F1(a,b,c,x) / \Gamma(c) for |x| < 1.
These routines compute the renormalized Gauss hypergeometric function 2F1(a_R + i a_I, a_R - i a_I, c, x) / \Gamma(c) for |x| < 1.
These routines compute the hypergeometric function 2F0(a,b,x). The series representation is a divergent hypergeometric series. However, for x < 0 we have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)
The generalized Laguerre polynomials are defined in terms of confluent hypergeometric functions as L^a_n(x) = ((a+1)_n / n!) 1F1(-n,a+1,x), and are sometimes referred to as the associated Laguerre polynomials. They are related to the plain Laguerre polynomials L_n(x) by L^0_n(x) = L_n(x) and L^k_n(x) = (-1)^k (d^k/dx^k) L_(n+k)(x). For more information see Abramowitz & Stegun, Chapter 22.
The functions described in this section are declared in the header file gsl_sf_laguerre.h.
These routines evaluate the generalized Laguerre polynomials L^a_1(x), L^a_2(x), L^a_3(x) using explicit representations.
These routines evaluate the generalized Laguerre polynomials L^a_n(x) for a > -1, n >= 0.
Lambert's W functions, W(x), are defined to be solutions of the equation W(x) \exp(W(x)) = x. This function has multiple branches for x < 0; however, it has only two real-valued branches. We define W_0(x) to be the principal branch, where W > -1 for x < 0, and W_{-1}(x) to be the other real branch, where W < -1 for x < 0. The Lambert functions are declared in the header file gsl_sf_lambert.h.
These compute the principal branch of the Lambert W function, W_0(x).
These compute the secondary real-valued branch of the Lambert W function, W_{-1}(x).
The Legendre Functions and Legendre Polynomials are described in Abramowitz & Stegun, Chapter 8. These functions are declared in the header file gsl_sf_legendre.h.
These functions evaluate the Legendre polynomials P_l(x) using explicit representations for l=1, 2, 3.
These functions evaluate the Legendre polynomial P_l(x) for a specific value of l, x subject to l >= 0, |x| <= 1
These functions compute an array of Legendre polynomials P_l(x), and optionally their derivatives dP_l(x)/dx, for l = 0, \dots, lmax, |x| <= 1
These routines compute the Legendre function Q_0(x) for x > -1, x != 1.
These routines compute the Legendre function Q_1(x) for x > -1, x != 1.
These routines compute the Legendre function Q_l(x) for x > -1, x != 1 and l >= 0.
The following functions compute the associated Legendre Polynomials
P_l^m(x). Note that this function grows combinatorially with
l and can overflow for l larger than about 150. There is
no trouble for small m, but overflow occurs when m and
l are both large. Rather than allow overflows, these functions
refuse to calculate P_l^m(x) and return GSL_EOVRFLW when
they can sense that l and m are too big.
If you want to calculate a spherical harmonic, then do not use
these functions. Instead use gsl_sf_legendre_sphPlm() below,
which uses a similar recursion, but with the normalized functions.
These routines compute the associated Legendre polynomial P_l^m(x) for m >= 0, l >= m, |x| <= 1.
These functions compute an array of Legendre polynomials P_l^m(x), and optionally their derivatives dP_l^m(x)/dx, for m >= 0, l = |m|, ..., lmax, |x| <= 1.
These routines compute the normalized associated Legendre polynomial $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$ suitable for use in spherical harmonics. The parameters must satisfy m >= 0, l >= m, |x| <= 1. Theses routines avoid the overflows that occur for the standard normalization of P_l^m(x).
These functions compute an array of normalized associated Legendre functions $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$, and optionally their derivatives, for m >= 0, l = |m|, ..., lmax, |x| <= 1.0
This function returns the size of result_array[] needed for the array versions of P_l^m(x), lmax - m + 1.
The Conical Functions P^\mu_{-(1/2)+i\lambda}(x) and Q^\mu_{-(1/2)+i\lambda} are described in Abramowitz & Stegun, Section 8.12.
These routines compute the irregular Spherical Conical Function P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.
These routines compute the regular Spherical Conical Function P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.
These routines compute the conical function P^0_{-1/2 + i \lambda}(x) for x > -1.
These routines compute the conical function P^1_{-1/2 + i \lambda}(x) for x > -1.
These routines compute the Regular Spherical Conical Function P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.
These routines compute the Regular Cylindrical Conical Function P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1.
The following spherical functions are specializations of Legendre functions which give the regular eigenfunctions of the Laplacian on a 3-dimensional hyperbolic space H3d. Of particular interest is the flat limit, \lambda \to \infty, \eta \to 0, \lambda\eta fixed.
These routines compute the zeroth radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space, L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta)) for \eta >= 0. In the flat limit this takes the form L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta).
These routines compute the first radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space, L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1} \sin(\lambda \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda \cot(\lambda\eta)) for \eta >= 0. In the flat limit this takes the form L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta).
These routines compute the l-th radial eigenfunction of the Laplacian on the 3-dimensional hyperbolic space \eta >= 0, l >= 0. In the flat limit this takes the form L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).
This function computes an array of radial eigenfunctions L^{H3d}_l(\lambda, \eta) for 0 <= l <= lmax.
Information on the properties of the Logarithm function can be found in Abramowitz & Stegun, Chapter 4. The functions described in this section are declared in the header file gsl_sf_log.h.
These routines compute the logarithm of x, \log(x), for x > 0.
These routines compute the logarithm of the magnitude of x, \log(|x|), for x \ne 0.
This routine computes the complex logarithm of z = z_r + i z_i. The results are returned as lnr, theta such that \exp(lnr + i \theta) = z_r + i z_i, where \theta lies in the range [-\pi,\pi].
These routines compute \log(1 + x) for x > -1 using an algorithm that is accurate for small x.
These routines compute \log(1 + x) - x for x > -1 using an algorithm that is accurate for small x.
The following functions are equivalent to the function gsl_pow_int
(see Small integer powers) with an error estimate. These functions are
declared in the header file gsl_sf_pow_int.h.
These routines compute the power x^n for integer n. The power is computed using the minimum number of multiplications. For example, x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications. For reasons of efficiency, these functions do not check for overflow or underflow conditions.
#include <gsl/gsl_sf_pow_int.h>
/* compute 3.0**12 */
double y = gsl_sf_pow_int(3.0, 12);
The polygamma functions of order m are defined by
\psi^{(m)}(x) = (d/dx)^m \psi(x) = (d/dx)^{m+1} \log(\Gamma(x))
where \psi(x) = \Gamma'(x)/\Gamma(x) is known as the digamma function. These functions are declared in the header file gsl_sf_psi.h.
These routines compute the digamma function \psi(n) for positive integer n. The digamma function is also called the Psi function.
These routines compute the digamma function \psi(x) for general x, x \ne 0.
These routines compute the real part of the digamma function on the line 1+i y, \Re[\psi(1 + i y)].
These routines compute the Trigamma function \psi'(n) for positive integer n.
These routines compute the Trigamma function \psi'(x) for general x.
These routines compute the polygamma function \psi^{(m)}(x) for m >= 0, x > 0.
The functions described in this section are declared in the header file gsl_sf_synchrotron.h.
These routines compute the first synchrotron function x \int_x^\infty dt K_{5/3}(t) for x >= 0.
These routines compute the second synchrotron function x K_{2/3}(x) for x >= 0.
The transport functions J(n,x) are defined by the integral representations J(n,x) := \int_0^x dt t^n e^t /(e^t - 1)^2. They are declared in the header file gsl_sf_transport.h.
These routines compute the transport function J(2,x).
These routines compute the transport function J(3,x).
These routines compute the transport function J(4,x).
These routines compute the transport function J(5,x).
The library includes its own trigonometric functions in order to provide consistency across platforms and reliable error estimates. These functions are declared in the header file gsl_sf_trig.h.
These routines compute the hypotenuse function \sqrt{x^2 + y^2} avoiding overflow and underflow.
These routines compute \sinc(x) = \sin(\pi x) / (\pi x) for any value of x.
This function computes the complex sine, \sin(z_r + i z_i) storing the real and imaginary parts in szr, szi.
This function computes the complex cosine, \cos(z_r + i z_i) storing the real and imaginary parts in szr, szi.
This function computes the logarithm of the complex sine, \log(\sin(z_r + i z_i)) storing the real and imaginary parts in szr, szi.
This function converts the polar coordinates (r,theta) to rectilinear coordinates (x,y), x = r\cos(\theta), y = r\sin(\theta).
This function converts the rectilinear coordinates (x,y) to polar coordinates (r,theta), such that x = r\cos(\theta), y = r\sin(\theta). The argument theta lies in the range [-\pi, \pi].
These routines force the angle theta to lie in the range (-\pi,\pi].
These routines force the angle theta to lie in the range [0, 2\pi).
This routine computes the sine of an angle x with an associated absolute error dx, \sin(x \pm dx). Note that this function is provided in the error-handling form only since its purpose is to compute the propagated error.
This routine computes the cosine of an angle x with an associated absolute error dx, \cos(x \pm dx). Note that this function is provided in the error-handling form only since its purpose is to compute the propagated error.
The Riemann zeta function is defined in Abramowitz & Stegun, Section 23.2. The functions described in this section are declared in the header file gsl_sf_zeta.h.
The Riemann zeta function is defined by the infinite sum \zeta(s) = \sum_{k=1}^\infty k^{-s}.
These routines compute the Riemann zeta function \zeta(n) for integer n, n \ne 1.
These routines compute the Riemann zeta function \zeta(s) for arbitrary s, s \ne 1.
For large positive argument, the Riemann zeta function approaches one. In this region the fractional part is interesting, and therefore we need a function to evaluate it explicitly.
These routines compute \zeta(n) - 1 for integer n, n \ne 1.
These routines compute \zeta(s) - 1 for arbitrary s, s \ne 1.
The Hurwitz zeta function is defined by \zeta(s,q) = \sum_0^\infty (k+q)^{-s}.
These routines compute the Hurwitz zeta function \zeta(s,q) for s > 1, q > 0.
The eta function is defined by \eta(s) = (1-2^{1-s}) \zeta(s).
These routines compute the eta function \eta(n) for integer n.
These routines compute the eta function \eta(s) for arbitrary s.
The following example demonstrates the use of the error handling form of the special functions, in this case to compute the Bessel function J_0(5.0),
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
int
main (void)
{
double x = 5.0;
gsl_sf_result result;
double expected = -0.17759677131433830434739701;
int status = gsl_sf_bessel_J0_e (x, &result);
printf ("status = %s\n", gsl_strerror(status));
printf ("J0(5.0) = %.18f\n"
" +/- % .18f\n",
result.val, result.err);
printf ("exact = %.18f\n", expected);
return status;
}
Here are the results of running the program,
$ ./a.out
status = success
J0(5.0) = -0.177596771314338292
+/- 0.000000000000000193
exact = -0.177596771314338292
The next program computes the same quantity using the natural form of the function. In this case the error term result.err and return status are not accessible.
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>
int
main (void)
{
double x = 5.0;
double expected = -0.17759677131433830434739701;
double y = gsl_sf_bessel_J0 (x);
printf ("J0(5.0) = %.18f\n", y);
printf ("exact = %.18f\n", expected);
return 0;
}
The results of the function are the same,
$ ./a.out
J0(5.0) = -0.177596771314338292
exact = -0.177596771314338292
The library follows the conventions of Abramowitz & Stegun where possible,
The following papers contain information on the algorithms used to compute the special functions,
The functions described in this chapter provide a simple 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 your functions in terms of vectors and matrices you can pass a single structure containing both data and dimensions as an argument without needing additional function parameters. The structures are compatible with the vector and matrix formats used by blas routines.
All the functions are available for each of the standard data-types.
The versions for double have the prefix gsl_block,
gsl_vector and gsl_matrix. Similarly the versions for
single-precision float arrays have the prefix
gsl_block_float, gsl_vector_float and
gsl_matrix_float. The full list of available types is given
below,
gsl_block double
gsl_block_float float
gsl_block_long_double long double
gsl_block_int int
gsl_block_uint unsigned int
gsl_block_long long
gsl_block_ulong unsigned long
gsl_block_short short
gsl_block_ushort unsigned short
gsl_block_char char
gsl_block_uchar unsigned char
gsl_block_complex complex double
gsl_block_complex_float complex float
gsl_block_complex_long_double complex long double
Corresponding types exist for the gsl_vector and
gsl_matrix functions.
For consistency all memory is allocated through a gsl_block
structure. The structure contains two components, the size of an area of
memory and a pointer to the memory. The gsl_block structure looks
like this,
typedef struct
{
size_t size;
double * data;
} gsl_block;
Vectors and matrices are made by slicing an underlying block. A slice is a set of elements formed from an initial offset and a combination of indices and step-sizes. In the case of a matrix the step-size for the column index represents the row-length. The step-size for a vector is known as the stride.
The functions for allocating and deallocating blocks are defined in gsl_block.h
The functions for allocating memory to a block follow the style of
malloc and free. In addition they also perform their own
error checking. If there is insufficient memory available to allocate a
block then the functions call the GSL error handler (with an error
number of GSL_ENOMEM) in addition to returning a null
pointer. Thus if you use the library error handler to abort your program
then it isn't necessary to check every alloc.
This function allocates memory for a block of n double-precision elements, returning a pointer to the block struct. The block is not initialized and so the values of its elements are undefined. Use the function
gsl_block_callocif you want to ensure that all the elements are initialized to zero.A null pointer is returned if insufficient memory is available to create the block.
This function allocates memory for a block and initializes all the elements of the block to zero.
This function frees the memory used by a block b previously allocated with
gsl_block_allocorgsl_block_calloc.
The library provides functions for reading and writing blocks to a file as binary data or formatted text.
This function writes the elements of the block b to the stream stream in binary format. The return value is 0 for success and
GSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
This function reads into the block b from the open stream stream in binary format. The block b must be preallocated with the correct length since the function uses the size of b to determine how many bytes to read. The return value is 0 for success and
GSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
This function writes the elements of the block b line-by-line to the stream stream using the format specifier format, which should be one of the
%g,%eor%fformats for floating point numbers and%dfor integers. The function returns 0 for success andGSL_EFAILEDif there was a problem writing to the file.
This function reads formatted data from the stream stream into the block b. The block b must be preallocated with the correct length since the function uses the size of b to determine how many numbers to read. The function returns 0 for success and
GSL_EFAILEDif there was a problem reading from the file.
The following program shows how to allocate a block,
#include <stdio.h>
#include <gsl/gsl_block.h>
int
main (void)
{
gsl_block * b = gsl_block_alloc (100);
printf ("length of block = %u\n", b->size);
printf ("block data address = %#x\n", b->data);
gsl_block_free (b);
return 0;
}
Here is the output from the program,
length of block = 100
block data address = 0x804b0d8
Vectors are defined by a gsl_vector structure which describes a
slice of a block. Different vectors can be created which point to the
same block. A vector slice is a set of equally-spaced elements of an
area of memory.
The gsl_vector structure contains five components, the
size, the stride, a pointer to the memory where the elements
are stored, data, a pointer to the block owned by the vector,
block, if any, and an ownership flag, owner. The structure
is very simple and looks like this,
typedef struct
{
size_t size;
size_t stride;
double * data;
gsl_block * block;
int owner;
} gsl_vector;
The size is simply the number of vector elements. The range of
valid indices runs from 0 to size-1. The stride is the
step-size from one element to the next in physical memory, measured in
units of the appropriate datatype. The pointer data gives the
location of the first element of the vector in memory. The pointer
block stores the location of the memory block in which the vector
elements are located (if any). If the vector owns this block then the
owner field is set to one and the block will be deallocated when the
vector is freed. If the vector points to a block owned by another
object then the owner field is zero and any underlying block will not be
deallocated with the vector.
The functions for allocating and accessing vectors are defined in gsl_vector.h
The functions for allocating memory to a vector follow the style of
malloc and free. In addition they also perform their own
error checking. If there is insufficient memory available to allocate a
vector then the functions call the GSL error handler (with an error
number of GSL_ENOMEM) in addition to returning a null
pointer. Thus if you use the library error handler to abort your program
then it isn't necessary to check every alloc.
This function creates a vector of length n, returning a pointer to a newly initialized vector struct. A new block is allocated for the elements of the vector, and stored in the block component of the vector struct. The block is “owned” by the vector, and will be deallocated when the vector is deallocated.
This function allocates memory for a vector of length n and initializes all the elements of the vector to zero.
This function frees a previously allocated vector v. If the vector was created using
gsl_vector_allocthen the block underlying the vector will also be deallocated. If the vector has been created from another object then the memory is still owned by that object and will not be deallocated.
Unlike fortran compilers, C compilers do not usually provide
support for range checking of vectors and matrices. Range checking is
available in the GNU C Compiler bounds-checking extension, but it is not
part of the default installation of GCC. The functions
gsl_vector_get and gsl_vector_set can perform portable
range checking for you and report an error if you attempt to access
elements outside the allowed range.
The functions for accessing the elements of a vector or matrix are
defined in gsl_vector.h and declared extern inline to
eliminate function-call overhead. You must compile your program with
the macro HAVE_INLINE defined to use these functions.
If necessary you can turn off range checking completely without
modifying any source files by recompiling your program with the
preprocessor definition GSL_RANGE_CHECK_OFF. Provided your
compiler supports inline functions the effect of turning off range
checking is to replace calls to gsl_vector_get(v,i) by
v->data[i*v->stride] and calls to gsl_vector_set(v,i,x) by
v->data[i*v->stride]=x. Thus there should be no performance
penalty for using the range checking functions when range checking is
turned off.
This function returns the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and 0 is returned.
This function sets the value of the i-th element of a vector v to x. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked.
These functions return a pointer to the i-th element of a vector v. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and a null pointer is returned.
This function sets all the elements of the vector v to the value x.
This function sets all the elements of the vector v to zero.
This function makes a basis vector by setting all the elements of the vector v to zero except for the i-th element which is set to one.
The library provides functions for reading and writing vectors to a file as binary data or formatted text.
This function writes the elements of the vector v to the stream stream in binary format. The return value is 0 for success and
GSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
This function reads into the vector v from the open stream stream in binary format. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many bytes to read. The return value is 0 for success and
GSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
This function writes the elements of the vector v line-by-line to the stream stream using the format specifier format, which should be one of the
%g,%eor%fformats for floating point numbers and%dfor integers. The function returns 0 for success andGSL_EFAILEDif there was a problem writing to the file.
This function reads formatted data from the stream stream into the vector v. The vector v must be preallocated with the correct length since the function uses the size of v to determine how many numbers to read. The function returns 0 for success and
GSL_EFAILEDif there was a problem reading from the file.
In addition to creating vectors from slices of blocks it is also possible to slice vectors and create vector views. For example, a subvector of another vector can be described with a view, or two views can be made which provide access to the even and odd elements of a vector.
A vector view is a temporary object, stored on the stack, which can be
used to operate on a subset of vector elements. Vector views can be
defined for both constant and non-constant vectors, using separate types
that preserve constness. A vector view has the type
gsl_vector_view and a constant vector view has the type
gsl_vector_const_view. In both cases the elements of the view
can be accessed as a gsl_vector using the vector component
of the view object. A pointer to a vector of type gsl_vector *
or const gsl_vector * can be obtained by taking the address of
this component with the & operator.
When using this pointer it is important to ensure that the view itself
remains in scope—the simplest way to do so is by always writing the
pointer as &view.vector, and never storing this value
in another variable.
These functions return a vector view of a subvector of another vector v. The start of the new vector is offset by offset elements from the start of the original vector. The new vector has n elements. Mathematically, the i-th element of the new vector v' is given by,
v'(i) = v->data[(offset + i)*v->stride]where the index i runs from 0 to
n-1.The
datapointer of the returned vector struct is set to null if the combined parameters (offset,n) overrun the end of the original vector.The new vector is only a view of the block underlying the original vector, v. The block containing the elements of v is not owned by the new vector. When the view goes out of scope the original vector v and its block will continue to exist. The original memory can only be deallocated by freeing the original vector. Of course, the original vector should not be deallocated while the view is still in use.
The function
gsl_vector_const_subvectoris equivalent togsl_vector_subvectorbut can be used for vectors which are declaredconst.
These functions return a vector view of a subvector of another vector v with an additional stride argument. The subvector is formed in the same way as for
gsl_vector_subvectorbut the new vector has n elements with a step-size of stride from one element to the next in the original vector. Mathematically, the i-th element of the new vector v' is given by,v'(i) = v->data[(offset + i*stride)*v->stride]where the index i runs from 0 to
n-1.Note that subvector views give direct access to the underlying elements of the original vector. For example, the following code will zero the even elements of the vector
vof lengthn, while leaving the odd elements untouched,gsl_vector_view v_even = gsl_vector_subvector_with_stride (v, 0, 2, n/2); gsl_vector_set_zero (&v_even.vector);A vector view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using
&view.vector. For example, the following code computes the norm of the odd elements ofvusing the blas routine dnrm2,gsl_vector_view v_odd = gsl_vector_subvector_with_stride (v, 1, 2, n/2); double r = gsl_blas_dnrm2 (&v_odd.vector);The function
gsl_vector_const_subvector_with_strideis equivalent togsl_vector_subvector_with_stridebut can be used for vectors which are declaredconst.
These functions return a vector view of the real parts of the complex vector v.
The function
gsl_vector_complex_const_realis equivalent togsl_vector_complex_realbut can be used for vectors which are declaredconst.
These functions return a vector view of the imaginary parts of the complex vector v.
The function
gsl_vector_complex_const_imagis equivalent togsl_vector_complex_imagbut can be used for vectors which are declaredconst.
These functions return a vector view of an array. The start of the new vector is given by base and has n elements. Mathematically, the i-th element of the new vector v' is given by,
v'(i) = base[i]where the index i runs from 0 to
n-1.The array containing the elements of v is not owned by the new vector view. When the view goes out of scope the original array will continue to exist. The original memory can only be deallocated by freeing the original pointer base. Of course, the original array should not be deallocated while the view is still in use.
The function
gsl_vector_const_view_arrayis equivalent togsl_vector_view_arraybut can be used for arrays which are declaredconst.
These functions return a vector view of an array base with an additional stride argument. The subvector is formed in the same way as for
gsl_vector_view_arraybut the new vector has n elements with a step-size of stride from one element to the next in the original array. Mathematically, the i-th element of the new vector v' is given by,v'(i) = base[i*stride]where the index i runs from 0 to
n-1.Note that the view gives direct access to the underlying elements of the original array. A vector view can be passed to any subroutine which takes a vector argument just as a directly allocated vector would be, using
&view.vector.The function
gsl_vector_const_view_array_with_strideis equivalent togsl_vector_view_array_with_stridebut can be used for arrays which are declaredconst.
Common operations on vectors such as addition and multiplication are available in the blas part of the library (see BLAS Support). However, it is useful to have a small number of utility functions which do not require the full blas code. The following functions fall into this category.
This function copies the elements of the vector src into the vector dest. The two vectors must have the same length.
This function exchanges the elements of the vectors v and w by copying. The two vectors must have the same length.
The following function can be used to exchange, or permute, the elements of a vector.
This function exchanges the i-th and j-th elements of the vector v in-place.
This function reverses the order of the elements of the vector v.
The following operations are only defined for real vectors.
This function adds the elements of vector b to the elements of vector a, a'_i = a_i + b_i. The two vectors must have the same length.
This function subtracts the elements of vector b from the elements of vector a, a'_i = a_i - b_i. The two vectors must have the same length.
This function multiplies the elements of vector a by the elements of vector b, a'_i = a_i * b_i. The two vectors must have the same length.
This function divides the elements of vector a by the elements of vector b, a'_i = a_i / b_i. The two vectors must have the same length.
This function multiplies the elements of vector a by the constant factor x, a'_i = x a_i.
This function adds the constant value x to the elements of the vector a, a'_i = a_i + x.
This function returns the maximum value in the vector v.
This function returns the minimum value in the vector v.
This function returns the minimum and maximum values in the vector v, storing them in min_out and max_out.
This function returns the index of the maximum value in the vector v. When there are several equal maximum elements then the lowest index is returned.
This function returns the index of the minimum value in the vector v. When there are several equal minimum elements then the lowest index is returned.
This function returns the indices of the minimum and maximum values in the vector v, storing them in imin and imax. When there are several equal minimum or maximum elements then the lowest indices are returned.
This function returns 1 if all the elements of the vector v are zero, and 0 otherwise.
This program shows how to allocate, initialize and read from a vector
using the functions gsl_vector_alloc, gsl_vector_set and
gsl_vector_get.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (3);
for (i = 0; i < 3; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
for (i = 0; i < 100; i++)
{
printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
}
return 0;
}
Here is the output from the program. The final loop attempts to read
outside the range of the vector v, and the error is trapped by
the range-checking code in gsl_vector_get.
$ ./a.out
v_0 = 1.23
v_1 = 2.23
v_2 = 3.23
gsl: vector_source.c:12: ERROR: index out of range
Default GSL error handler invoked.
Aborted (core dumped)
The next program shows how to write a vector to a file.
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (100);
for (i = 0; i < 100; i++)
{
gsl_vector_set (v, i, 1.23 + i);
}
{
FILE * f = fopen ("test.dat", "w");
gsl_vector_fprintf (f, v, "%.5g");
fclose (f);
}
return 0;
}
After running this program the file test.dat should contain the
elements of v, written using the format specifier
%.5g. The vector could then be read back in using the function
gsl_vector_fscanf (f, v) as follows:
#include <stdio.h>
#include <gsl/gsl_vector.h>
int
main (void)
{
int i;
gsl_vector * v = gsl_vector_alloc (10);
{
FILE * f = fopen ("test.dat", "r");
gsl_vector_fscanf (f, v);
fclose (f);
}
for (i = 0; i < 10; i++)
{
printf ("%g\n", gsl_vector_get(v, i));
}
return 0;
}
Matrices are defined by a gsl_matrix structure which describes a
generalized slice of a block. Like a vector it represents a set of
elements in an area of memory, but uses two indices instead of one.
The gsl_matrix structure contains six components, the two
dimensions of the matrix, a physical dimension, a pointer to the memory
where the elements of the matrix are stored, data, a pointer to
the block owned by the matrix block, if any, and an ownership
flag, owner. The physical dimension determines the memory layout
and can differ from the matrix dimension to allow the use of
submatrices. The gsl_matrix structure is very simple and looks
like this,
typedef struct
{
size_t size1;
size_t size2;
size_t tda;
double * data;
gsl_block * block;
int owner;
} gsl_matrix;
Matrices are stored in row-major order, meaning that each row of
elements forms a contiguous block in memory. This is the standard
“C-language ordering” of two-dimensional arrays. Note that fortran
stores arrays in column-major order. The number of rows is size1.
The range of valid row indices runs from 0 to size1-1. Similarly
size2 is the number of columns. The range of valid column indices
runs from 0 to size2-1. The physical row dimension tda, or
trailing dimension, specifies the size of a row of the matrix as
laid out in memory.
For example, in the following matrix size1 is 3, size2 is 4, and tda is 8. The physical memory layout of the matrix begins in the top left hand-corner and proceeds from left to right along each row in turn.
00 01 02 03 XX XX XX XX
10 11 12 13 XX XX XX XX
20 21 22 23 XX XX XX XX
Each unused memory location is represented by “XX”. The
pointer data gives the location of the first element of the matrix
in memory. The pointer block stores the location of the memory
block in which the elements of the matrix are located (if any). If the
matrix owns this block then the owner field is set to one and the
block will be deallocated when the matrix is freed. If the matrix is
only a slice of a block owned by another object then the owner field is
zero and any underlying block will not be freed.
The functions for allocating and accessing matrices are defined in gsl_matrix.h
The functions for allocating memory to a matrix follow the style of
malloc and free. They also perform their own error
checking. If there is insufficient memory available to allocate a vector
then the functions call the GSL error handler (with an error number of
GSL_ENOMEM) in addition to returning a null pointer. Thus if you
use the library error handler to abort your program then it isn't
necessary to check every alloc.
This function creates a matrix of size n1 rows by n2 columns, returning a pointer to a newly initialized matrix struct. A new block is allocated for the elements of the matrix, and stored in the block component of the matrix struct. The block is “owned” by the matrix, and will be deallocated when the matrix is deallocated.
This function allocates memory for a matrix of size n1 rows by n2 columns and initializes all the elements of the matrix to zero.
This function frees a previously allocated matrix m. If the matrix was created using
gsl_matrix_allocthen the block underlying the matrix will also be deallocated. If the matrix has been created from another object then the memory is still owned by that object and will not be deallocated.
The functions for accessing the elements of a matrix use the same range
checking system as vectors. You can turn off range checking by recompiling
your program with the preprocessor definition
GSL_RANGE_CHECK_OFF.
The elements of the matrix are stored in “C-order”, where the second
index moves continuously through memory. More precisely, the element
accessed by the function gsl_matrix_get(m,i,j) and
gsl_matrix_set(m,i,j,x) is
m->data[i * m->tda + j]
where tda is the physical row-length of the matrix.
This function returns the (i,j)-th element of a matrix m. If i or j lie outside the allowed range of 0 to n1-1 and 0 to n2-1 then the error handler is invoked and 0 is returned.
This function sets the value of the (i,j)-th element of a matrix m to x. If i or j lies outside the allowed range of 0 to n1-1 and 0 to n2-1 then the error handler is invoked.
These functions return a pointer to the (i,j)-th element of a matrix m. If i or j lie outside the allowed range of 0 to n1-1 and 0 to n2-1 then the error handler is invoked and a null pointer is returned.
This function sets all the elements of the matrix m to the value x.
This function sets all the elements of the matrix m to zero.
This function sets the elements of the matrix m to the corresponding elements of the identity matrix, m(i,j) = \delta(i,j), i.e. a unit diagonal with all off-diagonal elements zero. This applies to both square and rectangular matrices.
The library provides functions for reading and writing matrices to a file as binary data or formatted text.
This function writes the elements of the matrix m to the stream stream in binary format. The return value is 0 for success and
GSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
This function reads into the matrix m from the open stream stream in binary format. The matrix m must be preallocated with the correct dimensions since the function uses the size of m to determine how many bytes to read. The return value is 0 for success and
GSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
This function writes the elements of the matrix m line-by-line to the stream stream using the format specifier format, which should be one of the
%g,%eor%fformats for floating point numbers and%dfor integers. The function returns 0 for success andGSL_EFAILEDif there was a problem writing to the file.
This function reads formatted data from the stream stream into the matrix m. The matrix m must be preallocated with the correct dimensions since the function uses the size of m to determine how many numbers to read. The function returns 0 for success and
GSL_EFAILEDif there was a problem reading from the file.
A matrix view is a temporary object, stored on the stack, which can be
used to operate on a subset of matrix elements. Matrix views can be
defined for both constant and non-constant matrices using separate types
that preserve constness. A matrix view has the type
gsl_matrix_view and a constant matrix view has the type
gsl_matrix_const_view. In both cases the elements of the view
can by accessed using the matrix component of the view object. A
pointer gsl_matrix * or const gsl_matrix * can be obtained
by taking the address of the matrix component with the &
operator. In addition to matrix views it is also possible to create
vector views of a matrix, such as row or column views.
These functions return a matrix view of a submatrix of the matrix m. The upper-left element of the submatrix is the element (k1,k2) of the original matrix. The submatrix has n1 rows and n2 columns. The physical number of columns in memory given by tda is unchanged. Mathematically, the (i,j)-th element of the new matrix is given by,
m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j]where the index i runs from 0 to
n1-1and the index j runs from 0 ton2-1.The
datapointer of the returned matrix struct is set to null if the combined parameters (i,j,n1,n2,tda) overrun the ends of the original matrix.The new matrix view is only a view of the block underlying the existing matrix, m. The block containing the elements of m is not owned by the new matrix view. When the view goes out of scope the original matrix m and its block will continue to exist. The original memory can only be deallocated by freeing the original matrix. Of course, the original matrix should not be deallocated while the view is still in use.
The function
gsl_matrix_const_submatrixis equivalent togsl_matrix_submatrixbut can be used for matrices which are declaredconst.
These functions return a matrix view of the array base. The matrix has n1 rows and n2 columns. The physical number of columns in memory is also given by n2. Mathematically, the (i,j)-th element of the new matrix is given by,
m'(i,j) = base[i*n2 + j]where the index i runs from 0 to
n1-1and the index j runs from 0 ton2-1.The new matrix is only a view of the array base. When the view goes out of scope the original array base will continue to exist. The original memory can only be deallocated by freeing the original array. Of course, the original array should not be deallocated while the view is still in use.
The function
gsl_matrix_const_view_arrayis equivalent togsl_matrix_view_arraybut can be used for matrices which are declaredconst.
These functions return a matrix view of the array base with a physical number of columns tda which may differ from the corresponding dimension of the matrix. The matrix has n1 rows and n2 columns, and the physical number of columns in memory is given by tda. Mathematically, the (i,j)-th element of the new matrix is given by,
m'(i,j) = base[i*tda + j]where the index i runs from 0 to
n1-1and the index j runs from 0 ton2-1.The new matrix is only a view of the array base. When the view goes out of scope the original array base will continue to exist. The original memory can only be deallocated by freeing the original array. Of course, the original array should not be deallocated while the view is still in use.
The function
gsl_matrix_const_view_array_with_tdais equivalent togsl_matrix_view_array_with_tdabut can be used for matrices which are declaredconst.
These functions return a matrix view of the vector v. The matrix has n1 rows and n2 columns. The vector must have unit stride. The physical number of columns in memory is also given by n2. Mathematically, the (i,j)-th element of the new matrix is given by,
m'(i,j) = v->data[i*n2 + j]where the index i runs from 0 to
n1-1and the index j runs from 0 ton2-1.The new matrix is only a view of the vector v. When the view goes out of scope the original vector v will continue to exist. The original memory can only be deallocated by freeing the original vector. Of course, the original vector should not be deallocated while the view is still in use.
The function
gsl_matrix_const_view_vectoris equivalent togsl_matrix_view_vectorbut can be used for matrices which are declaredconst.
These functions return a matrix view of the vector v with a physical number of columns tda which may differ from the corresponding matrix dimension. The vector must have unit stride. The matrix has n1 rows and n2 columns, and the physical number of columns in memory is given by tda. Mathematically, the (i,j)-th element of the new matrix is given by,
m'(i,j) = v->data[i*tda + j]where the index i runs from 0 to
n1-1and the index j runs from 0 ton2-1.The new matrix is only a view of the vector v. When the view goes out of scope the original vector v will continue to exist. The original memory can only be deallocated by freeing the original vector. Of course, the original vector should not be deallocated while the view is still in use.
The function
gsl_matrix_const_view_vector_with_tdais equivalent togsl_matrix_view_vector_with_tdabut can be used for matrices which are declaredconst.
In general there are two ways to access an object, by reference or by copying. The functions described in this section create vector views which allow access to a row or column of a matrix by reference. Modifying elements of the view is equivalent to modifying the matrix, since both the vector view and the matrix point to the same memory block.
These functions return a vector view of the i-th row of the matrix m. The
datapointer of the new vector is set to null if i is out of range.The function
gsl_vector_const_rowis equivalent togsl_matrix_rowbut can be used for matrices which are declaredconst.
These functions return a vector view of the j-th column of the matrix m. The
datapointer of the new vector is set to null if j is out of range.The function
gsl_vector_const_columnis equivalent togsl_matrix_columnbut can be used for matrices which are declaredconst.
These functions returns a vector view of the diagonal of the matrix m. The matrix m is not required to be square. For a rectangular matrix the length of the diagonal is the same as the smaller dimension of the matrix.
The function
gsl_matrix_const_diagonalis equivalent togsl_matrix_diagonalbut can be used for matrices which are declaredconst.
These functions return a vector view of the k-th subdiagonal of the matrix m. The matrix m is not required to be square. The diagonal of the matrix corresponds to k = 0.
The function
gsl_matrix_const_subdiagonalis equivalent togsl_matrix_subdiagonalbut can be used for matrices which are declaredconst.
These functions return a vector view of the k-th superdiagonal of the matrix m. The matrix m is not required to be square. The diagonal of the matrix corresponds to k = 0.
The function
gsl_matrix_const_superdiagonalis equivalent togsl_matrix_superdiagonalbut can be used for matrices which are declaredconst.
This function copies the elements of the matrix src into the matrix dest. The two matrices must have the same size.
This function exchanges the elements of the matrices m1 and m2 by copying. The two matrices must have the same size.
The functions described in this section copy a row or column of a matrix
into a vector. This allows the elements of the vector and the matrix to
be modified independently. Note that if the matrix and the vector point
to overlapping regions of memory then the result will be undefined. The
same effect can be achieved with more generality using
gsl_vector_memcpy with vector views of rows and columns.
This function copies the elements of the i-th row of the matrix m into the vector v. The length of the vector must be the same as the length of the row.
This function copies the elements of the j-th column of the matrix m into the vector v. The length of the vector must be the same as the length of the column.
This function copies the elements of the vector v into the i-th row of the matrix m. The length of the vector must be the same as the length of the row.
This function copies the elements of the vector v into the j-th column of the matrix m. The length of the vector must be the same as the length of the column.
The following functions can be used to exchange the rows and columns of a matrix.
This function exchanges the i-th and j-th rows of the matrix m in-place.
This function exchanges the i-th and j-th columns of the matrix m in-place.
This function exchanges the i-th row and j-th column of the matrix m in-place. The matrix must be square for this operation to be possible.
This function makes the matrix dest the transpose of the matrix src by copying the elements of src into dest. This function works for all matrices provided that the dimensions of the matrix dest match the transposed dimensions of the matrix src.
This function replaces the matrix m by its transpose by copying the elements of the matrix in-place. The matrix must be square for this operation to be possible.
The following operations are defined for real and complex matrices.
This function adds the elements of matrix b to the elements of matrix a, a'(i,j) = a(i,j) + b(i,j). The two matrices must have the same dimensions.
This function subtracts the elements of matrix b from the elements of matrix a, a'(i,j) = a(i,j) - b(i,j). The two matrices must have the same dimensions.
This function multiplies the elements of matrix a by the elements of matrix b, a'(i,j) = a(i,j) * b(i,j). The two matrices must have the same dimensions.
This function divides the elements of matrix a by the elements of matrix b, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the same dimensions.
This function multiplies the elements of matrix a by the constant factor x, a'(i,j) = x a(i,j).
This function adds the constant value x to the elements of the matrix a, a'(i,j) = a(i,j) + x.
The following operations are only defined for real matrices.
This function returns the maximum value in the matrix m.
This function returns the minimum value in the matrix m.
This function returns the minimum and maximum values in the matrix m, storing them in min_out and max_out.
This function returns the indices of the maximum value in the matrix m, storing them in imax and jmax. When there are several equal maximum elements then the first element found is returned, searching in row-major order.
This function returns the indices of the minimum value in the matrix m, storing them in imin and jmin. When there are several equal minimum elements then the first element found is returned, searching in row-major order.
This function returns the indices of the minimum and maximum values in the matrix m, storing them in (imin,jmin) and (imax,jmax). When there are several equal minimum or maximum elements then the first elements found are returned, searching in row-major order.
This function returns 1 if all the elements of the matrix m are zero, and 0 otherwise.
The program below shows how to allocate, initialize and read from a matrix
using the functions gsl_matrix_alloc, gsl_matrix_set and
gsl_matrix_get.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int
main (void)
{
int i, j;
gsl_matrix * m = gsl_matrix_alloc (10, 3);
for (i = 0; i < 10; i++)
for (j = 0; j < 3; j++)
gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
for (i = 0; i < 100; i++)
for (j = 0; j < 3; j++)
printf ("m(%d,%d) = %g\n", i, j,
gsl_matrix_get (m, i, j));
return 0;
}
Here is the output from the program. The final loop attempts to read
outside the range of the matrix m, and the error is trapped by
the range-checking code in gsl_matrix_get.
$ ./a.out
m(0,0) = 0.23
m(0,1) = 1.23
m(0,2) = 2.23
m(1,0) = 100.23
m(1,1) = 101.23
m(1,2) = 102.23
...
m(9,2) = 902.23
gsl: matrix_source.c:13: ERROR: first index out of range
Default GSL error handler invoked.
Aborted (core dumped)
The next program shows how to write a matrix to a file.
#include <stdio.h>
#include <gsl/gsl_matrix.h>
int
main (void)
{
int i, j, k = 0;
gsl_matrix * m = gsl_matrix_alloc (100, 100);
gsl_matrix * a = gsl_matrix_alloc (100, 100);
for (i = 0; i < 100; i++)
for (j = 0; j < 100; j++)
gsl_matrix_set (m, i, j, 0.23 + i + j);
{
FILE * f = fopen ("test.dat", "wb");
gsl_matrix_fwrite (f, m);
fclose (f);
}
{
FILE * f = fopen ("test.dat", "rb");
gsl_matrix_fread (f, a);
fclose (f);
}
for (i = 0; i < 100; i++)
for (j = 0; j < 100; j++)
{
double mij = gsl_matrix_get (m, i, j);
double aij = gsl_matrix_get (a, i, j);
if (mij != aij) k++;
}
printf ("differences = %d (should be zero)\n", k);
return (k > 0);
}
After running this program the file test.dat should contain the
elements of m, written in binary format. The matrix which is read
back in using the function gsl_matrix_fread should be exactly
equal to the original matrix.
The following program demonstrates the use of vector views. The program computes the column norms of a matrix.
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
int
main (void)
{
size_t i,j;
gsl_matrix *m = gsl_matrix_alloc (10, 10);
for (i = 0; i < 10; i++)
for (j = 0; j < 10; j++)
gsl_matrix_set (m, i, j, sin (i) + cos (j));
for (j = 0; j < 10; j++)
{
gsl_vector_view column = gsl_matrix_column (m, j);
double d;
d = gsl_blas_dnrm2 (&column.vector);
printf ("matrix column %d, norm = %g\n", j, d);
}
gsl_matrix_free (m);
return 0;
}
Here is the output of the program,
$ ./a.out
matrix column 0, norm = 4.31461
matrix column 1, norm = 3.1205
matrix column 2, norm = 2.19316
matrix column 3, norm = 3.26114
matrix column 4, norm = 2.53416
matrix column 5, norm = 2.57281
matrix column 6, norm = 4.20469
matrix column 7, norm = 3.65202
matrix column 8, norm = 2.08524
matrix column 9, norm = 3.07313
The results can be confirmed using gnu octave,
$ octave
GNU Octave, version 2.0.16.92
octave> m = sin(0:9)' * ones(1,10)
+ ones(10,1) * cos(0:9);
octave> sqrt(sum(m.^2))
ans =
4.3146 3.1205 2.1932 3.2611 2.5342 2.5728
4.2047 3.6520 2.0852 3.0731
The block, vector and matrix objects in GSL follow the valarray
model of C++. A description of this model can be found in the following
reference,
This chapter describes functions for creating and manipulating permutations. A permutation p is represented by an array of n integers in the range 0 to n-1, where each value p_i occurs once and only once. The application of a permutation p to a vector v yields a new vector v' where v'_i = v_{p_i}. For example, the array (0,1,3,2) represents a permutation which exchanges the last two elements of a four element vector. The corresponding identity permutation is (0,1,2,3).
Note that the permutations produced by the linear algebra routines correspond to the exchange of matrix columns, and so should be considered as applying to row-vectors in the form v' = v P rather than column-vectors, when permuting the elements of a vector.
The functions described in this chapter are defined in the header file gsl_permutation.h.
A permutation is defined by a structure containing two components, the size
of the permutation and a pointer to the permutation array. The elements
of the permutation array are all of type size_t. The
gsl_permutation structure looks like this,
typedef struct
{
size_t size;
size_t * data;
} gsl_permutation;
This function allocates memory for a new permutation of size n. The permutation is not initialized and its elements are undefined. Use the function
gsl_permutation_callocif you want to create a permutation which is initialized to the identity. A null pointer is returned if insufficient memory is available to create the permutation.
This function allocates memory for a new permutation of size n and initializes it to the identity. A null pointer is returned if insufficient memory is available to create the permutation.
This function initializes the permutation p to the identity, i.e. (0,1,2,...,n-1).
This function frees all the memory used by the permutation p.
This function copies the elements of the permutation src into the permutation dest. The two permutations must have the same size.
The following functions can be used to access and manipulate permutations.
This function returns the value of the i-th element of the permutation p. If i lies outside the allowed range of 0 to n-1 then the error handler is invoked and 0 is returned.
This function exchanges the i-th and j-th elements of the permutation p.
This function returns the size of the permutation p.
This function returns a pointer to the array of elements in the permutation p.
This function checks that the permutation p is valid. The n elements should contain each of the numbers 0 to n-1 once and only once.
This function computes the inverse of the permutation p, storing the result in inv.
This function advances the permutation p to the next permutation in lexicographic order and returns
GSL_SUCCESS. If no further permutations are available it returnsGSL_FAILUREand leaves p unmodified. Starting with the identity permutation and repeatedly applying this function will iterate through all possible permutations of a given order.
This function steps backwards from the permutation p to the previous permutation in lexicographic order, returning
GSL_SUCCESS. If no previous permutation is available it returnsGSL_FAILUREand leaves p unmodified.
This function applies the permutation p to the array data of size n with stride stride.
This function applies the inverse of the permutation p to the array data of size n with stride stride.
This function applies the permutation p to the elements of the vector v, considered as a row-vector acted on by a permutation matrix from the right, v' = v P. The j-th column of the permutation matrix P is given by the p_j-th column of the identity matrix. The permutation p and the vector v must have the same length.
This function applies the inverse of the permutation p to the elements of the vector v, considered as a row-vector acted on by an inverse permutation matrix from the right, v' = v P^T. Note that for permutation matrices the inverse is the same as the transpose. The j-th column of the permutation matrix P is given by the p_j-th column of the identity matrix. The permutation p and the vector v must have the same length.
This function combines the two permutations pa and pb into a single permutation p, where p = pa . pb. The permutation p is equivalent to applying pb first and then pa.
The library provides functions for reading and writing permutations to a file as binary data or formatted text.
This function writes the elements of the permutation p to the stream stream in binary format. The function returns
GSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
This function reads into the permutation p from the open stream stream in binary format. The permutation p must be preallocated with the correct length since the function uses the size of p to determine how many bytes to read. The function returns
GSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
This function writes the elements of the permutation p line-by-line to the stream stream using the format specifier format, which should be suitable for a type of size_t. On a GNU system the type modifier
Zrepresentssize_t, so"%Zu\n"is a suitable format. The function returnsGSL_EFAILEDif there was a problem writing to the file.
This function reads formatted data from the stream stream into the permutation p. The permutation p must be preallocated with the correct length since the function uses the size of p to determine how many numbers to read. The function returns
GSL_EFAILEDif there was a problem reading from the file.
A permutation can be represented in both linear and cyclic notations. The functions described in this section convert between the two forms. The linear notation is an index mapping, and has already been described above. The cyclic notation expresses a permutation as a series of circular rearrangements of groups of elements, or cycles.
For example, under the cycle (1 2 3), 1 is replaced by 2, 2 is replaced by 3 and 3 is replaced by 1 in a circular fashion. Cycles of different sets of elements can be combined independently, for example (1 2 3) (4 5) combines the cycle (1 2 3) with the cycle (4 5), which is an exchange of elements 4 and 5. A cycle of length one represents an element which is unchanged by the permutation and is referred to as a singleton.
It can be shown that every permutation can be decomposed into combinations of cycles. The decomposition is not unique, but can always be rearranged into a standard canonical form by a reordering of elements. The library uses the canonical form defined in Knuth's Art of Computer Programming (Vol 1, 3rd Ed, 1997) Section 1.3.3, p.178.
The procedure for obtaining the canonical form given by Knuth is,
For example, the linear representation (2 4 3 0 1) is represented as (1 4) (0 2 3) in canonical form. The permutation corresponds to an exchange of elements 1 and 4, and rotation of elements 0, 2 and 3.
The important property of the canonical form is that it can be reconstructed from the contents of each cycle without the brackets. In addition, by removing the brackets it can be considered as a linear representation of a different permutation. In the example given above the permutation (2 4 3 0 1) would become (1 4 0 2 3). This mapping has many applications in the theory of permutations.
This function computes the canonical form of the permutation p and stores it in the output argument q.
This function converts a permutation q in canonical form back into linear form storing it in the output argument p.
This function counts the number of inversions in the permutation p. An inversion is any pair of elements that are not in order. For example, the permutation 2031 has three inversions, corresponding to the pairs (2,0) (2,1) and (3,1). The identity permutation has no inversions.
This function counts the number of cycles in the permutation p, given in linear form.
This function counts the number of cycles in the permutation q, given in canonical form.
The example program below creates a random permutation (by shuffling the elements of the identity) and finds its inverse.
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
int
main (void)
{
const size_t N = 10;
const gsl_rng_type * T;
gsl_rng * r;
gsl_permutation * p = gsl_permutation_alloc (N);
gsl_permutation * q = gsl_permutation_alloc (N);
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
printf ("initial permutation:");
gsl_permutation_init (p);
gsl_permutation_fprintf (stdout, p, " %u");
printf ("\n");
printf (" random permutation:");
gsl_ran_shuffle (r, p->data, N, sizeof(size_t));
gsl_permutation_fprintf (stdout, p, " %u");
printf ("\n");
printf ("inverse permutation:");
gsl_permutation_inverse (q, p);
gsl_permutation_fprintf (stdout, q, " %u");
printf ("\n");
return 0;
}
Here is the output from the program,
$ ./a.out
initial permutation: 0 1 2 3 4 5 6 7 8 9
random permutation: 1 3 5 2 7 6 0 4 9 8
inverse permutation: 6 0 3 1 7 2 5 4 9 8
The random permutation p[i] and its inverse q[i] are
related through the identity p[q[i]] = i, which can be verified
from the output.
The next example program steps forwards through all possible third order permutations, starting from the identity,
#include <stdio.h>
#include <gsl/gsl_permutation.h>
int
main (void)
{
gsl_permutation * p = gsl_permutation_alloc (3);
gsl_permutation_init (p);
do
{
gsl_permutation_fprintf (stdout, p, " %u");
printf ("\n");
}
while (gsl_permutation_next(p) == GSL_SUCCESS);
return 0;
}
Here is the output from the program,
$ ./a.out
0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0
The permutations are generated in lexicographic order. To reverse the
sequence, begin with the final permutation (which is the reverse of the
identity) and replace gsl_permutation_next with
gsl_permutation_prev.
The subject of permutations is covered extensively in Knuth's Sorting and Searching,
For the definition of the canonical form see,
This chapter describes functions for creating and manipulating combinations. A combination c is represented by an array of k integers in the range 0 to n-1, where each value c_i occurs at most once. The combination c corresponds to indices of k elements chosen from an n element vector. Combinations are useful for iterating over all k-element subsets of a set.
The functions described in this chapter are defined in the header file gsl_combination.h.
A combination is defined by a structure containing three components, the
values of n and k, and a pointer to the combination array.
The elements of the combination array are all of type size_t, and
are stored in increasing order. The gsl_combination structure
looks like this,
typedef struct
{
size_t n;
size_t k;
size_t *data;
} gsl_combination;
This function allocates memory for a new combination with parameters n, k. The combination is not initialized and its elements are undefined. Use the function
gsl_combination_callocif you want to create a combination which is initialized to the lexicographically first combination. A null pointer is returned if insufficient memory is available to create the combination.
This function allocates memory for a new combination with parameters n, k and initializes it to the lexicographically first combination. A null pointer is returned if insufficient memory is available to create the combination.
This function initializes the combination c to the lexicographically first combination, i.e. (0,1,2,...,k-1).
This function initializes the combination c to the lexicographically last combination, i.e. (n-k,n-k+1,...,n-1).
This function frees all the memory used by the combination c.
This function copies the elements of the combination src into the combination dest. The two combinations must have the same size.
The following function can be used to access the elements of a combination.
This function returns the value of the i-th element of the combination c. If i lies outside the allowed range of 0 to k-1 then the error handler is invoked and 0 is returned.
This function returns the range (n) of the combination c.
This function returns the number of elements (k) in the combination c.
This function returns a pointer to the array of elements in the combination c.
This function checks that the combination c is valid. The k elements should lie in the range 0 to n-1, with each value occurring once at most and in increasing order.
This function advances the combination c to the next combination in lexicographic order and returns
GSL_SUCCESS. If no further combinations are available it returnsGSL_FAILUREand leaves c unmodified. Starting with the first combination and repeatedly applying this function will iterate through all possible combinations of a given order.
This function steps backwards from the combination c to the previous combination in lexicographic order, returning
GSL_SUCCESS. If no previous combination is available it returnsGSL_FAILUREand leaves c unmodified.
The library provides functions for reading and writing combinations to a file as binary data or formatted text.
This function writes the elements of the combination c to the stream stream in binary format. The function returns
GSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
This function reads elements from the open stream stream into the combination c in binary format. The combination c must be preallocated with correct values of n and k since the function uses the size of c to determine how many bytes to read. The function returns
GSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
This function writes the elements of the combination c line-by-line to the stream stream using the format specifier format, which should be suitable for a type of size_t. On a GNU system the type modifier
Zrepresentssize_t, so"%Zu\n"is a suitable format. The function returnsGSL_EFAILEDif there was a problem writing to the file.
This function reads formatted data from the stream stream into the combination c. The combination c must be preallocated with correct values of n and k since the function uses the size of c to determine how many numbers to read. The function returns
GSL_EFAILEDif there was a problem reading from the file.
The example program below prints all subsets of the set {0,1,2,3} ordered by size. Subsets of the same size are ordered lexicographically.
#include <stdio.h>
#include <gsl/gsl_combination.h>
int
main (void)
{
gsl_combination * c;
size_t i;
printf ("All subsets of {0,1,2,3} by size:\n") ;
for (i = 0; i <= 4; i++)
{
c = gsl_combination_calloc (4, i);
do
{
printf ("{");
gsl_combination_fprintf (stdout, c, " %u");
printf (" }\n");
}
while (gsl_combination_next (c) == GSL_SUCCESS);
gsl_combination_free (c);
}
return 0;
}
Here is the output from the program,
$ ./a.out
All subsets of {0,1,2,3} by size:
{ }
{ 0 }
{ 1 }
{ 2 }
{ 3 }
{ 0 1 }
{ 0 2 }
{ 0 3 }
{ 1 2 }
{ 1 3 }
{ 2 3 }
{ 0 1 2 }
{ 0 1 3 }
{ 0 2 3 }
{ 1 2 3 }
{ 0 1 2 3 }
All 16 subsets are generated, and the subsets of each size are sorted lexicographically.
Further information on combinations can be found in,
This chapter describes functions for sorting data, both directly and indirectly (using an index). All the functions use the heapsort algorithm. Heapsort is an O(N \log N) algorithm which operates in-place and does not require any additional storage. It also provides consistent performance, the running time for its worst-case (ordered data) being not significantly longer than the average and best cases. Note that the heapsort algorithm does not preserve the relative ordering of equal elements—it is an unstable sort. However the resulting order of equal elements will be consistent across different platforms when using these functions.
The following function provides a simple alternative to the standard
library function qsort. It is intended for systems lacking
qsort, not as a replacement for it. The function qsort
should be used whenever possible, as it will be faster and can provide
stable ordering of equal elements. Documentation for qsort is
available in the GNU C Library Reference Manual.
The functions described in this section are defined in the header file gsl_heapsort.h.
This function sorts the count elements of the array array, each of size size, into ascending order using the comparison function compare. The type of the comparison function is defined by,
int (*gsl_comparison_fn_t) (const void * a, const void * b)A comparison function should return a negative integer if the first argument is less than the second argument,
0if the two arguments are equal and a positive integer if the first argument is greater than the second argument.For example, the following function can be used to sort doubles into ascending numerical order.
int compare_doubles (const double * a, const double * b) { if (*a > *b) return 1; else if (*a < *b) return -1; else return 0; }The appropriate function call to perform the sort is,
gsl_heapsort (array, count, sizeof(double), compare_doubles);Note that unlike
qsortthe heapsort algorithm cannot be made into a stable sort by pointer arithmetic. The trick of comparing pointers for equal elements in the comparison function does not work for the heapsort algorithm. The heapsort algorithm performs an internal rearrangement of the data which destroys its initial ordering.
This function indirectly sorts the count elements of the array array, each of size size, into ascending order using the comparison function compare. The resulting permutation is stored in p, an array of length n. The elements of p give the index of the array element which would have been stored in that position if the array had been sorted in place. The first element of p gives the index of the least element in array, and the last element of p gives the index of the greatest element in array. The array itself is not changed.
The following functions will sort the elements of an array or vector,
either directly or indirectly. They are defined for all real and integer
types using the normal suffix rules. For example, the float
versions of the array functions are gsl_sort_float and
gsl_sort_float_index. The corresponding vector functions are
gsl_sort_vector_float and gsl_sort_vector_float_index. The
prototypes are available in the header files gsl_sort_float.h
gsl_sort_vector_float.h. The complete set of prototypes can be
included using the header files gsl_sort.h and
gsl_sort_vector.h.
There are no functions for sorting complex arrays or vectors, since the ordering of complex numbers is not uniquely defined. To sort a complex vector by magnitude compute a real vector containing the magnitudes of the complex elements, and sort this vector indirectly. The resulting index gives the appropriate ordering of the original complex vector.
This function sorts the n elements of the array data with stride stride into ascending numerical order.
This function sorts the elements of the vector v into ascending numerical order.
This function indirectly sorts the n elements of the array data with stride stride into ascending order, storing the resulting permutation in p. The array p must be allocated with a sufficient length to store the n elements of the permutation. The elements of p give the index of the array element which would have been stored in that position if the array had been sorted in place. The array data is not changed.
This function indirectly sorts the elements of the vector v into ascending order, storing the resulting permutation in p. The elements of p give the index of the vector element which would have been stored in that position if the vector had been sorted in place. The first element of p gives the index of the least element in v, and the last element of p gives the index of the greatest element in v. The vector v is not changed.
The functions described in this section select the k smallest or largest elements of a data set of size N. The routines use an O(kN) direct insertion algorithm which is suited to subsets that are small compared with the total size of the dataset. For example, the routines are useful for selecting the 10 largest values from one million data points, but not for selecting the largest 100,000 values. If the subset is a significant part of the total dataset it may be faster to sort all the elements of the dataset directly with an O(N \log N) algorithm and obtain the smallest or largest values that way.
This function copies the k smallest elements of the array src, of size n and stride stride, in ascending numerical order into the array dest. The size k of the subset must be less than or equal to n. The data src is not modified by this operation.
This function copies the k largest elements of the array src, of size n and stride stride, in descending numerical order into the array dest. k must be less than or equal to n. The data src is not modified by this operation.
These functions copy the k smallest or largest elements of the vector v into the array dest. k must be less than or equal to the length of the vector v.
The following functions find the indices of the k smallest or largest elements of a dataset,
This function stores the indices of the k smallest elements of the array src, of size n and stride stride, in the array p. The indices are chosen so that the corresponding data is in ascending numerical order. k must be less than or equal to n. The data src is not modified by this operation.
This function stores the indices of the k largest elements of the array src, of size n and stride stride, in the array p. The indices are chosen so that the corresponding data is in descending numerical order. k must be less than or equal to n. The data src is not modified by this operation.
These functions store the indices of the k smallest or largest elements of the vector v in the array p. k must be less than or equal to the length of the vector v.
The rank of an element is its order in the sorted data. The rank is the inverse of the index permutation, p. It can be computed using the following algorithm,
for (i = 0; i < p->size; i++)
{
size_t pi = p->data[i];
rank->data[pi] = i;
}
This can be computed directly from the function
gsl_permutation_inverse(rank,p).
The following function will print the rank of each element of the vector v,
void
print_rank (gsl_vector * v)
{
size_t i;
size_t n = v->size;
gsl_permutation * perm = gsl_permutation_alloc(n);
gsl_permutation * rank = gsl_permutation_alloc(n);
gsl_sort_vector_index (perm, v);
gsl_permutation_inverse (rank, perm);
for (i = 0; i < n; i++)
{
double vi = gsl_vector_get(v, i);
printf ("element = %d, value = %g, rank = %d\n",
i, vi, rank->data[i]);
}
gsl_permutation_free (perm);
gsl_permutation_free (rank);
}
The following example shows how to use the permutation p to print the elements of the vector v in ascending order,
gsl_sort_vector_index (p, v);
for (i = 0; i < v->size; i++)
{
double vpi = gsl_vector_get (v, p->data[i]);
printf ("order = %d, value = %g\n", i, vpi);
}
The next example uses the function gsl_sort_smallest to select
the 5 smallest numbers from 100000 uniform random variates stored in an
array,
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>
int
main (void)
{
const gsl_rng_type * T;
gsl_rng * r;
size_t i, k = 5, N = 100000;
double * x = malloc (N * sizeof(double));
double * small = malloc (k * sizeof(double));
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
for (i = 0; i < N; i++)
{
x[i] = gsl_rng_uniform(r);
}
gsl_sort_smallest (small, k, x, 1, N);
printf ("%d smallest values from %d\n", k, N);
for (i = 0; i < k; i++)
{
printf ("%d: %.18f\n", i, small[i]);
}
return 0;
}
The output lists the 5 smallest values, in ascending order,
$ ./a.out
5 smallest values from 100000
0: 0.000003489200025797
1: 0.000008199829608202
2: 0.000008953968062997
3: 0.000010712770745158
4: 0.000033531803637743
The subject of sorting is covered extensively in Knuth's Sorting and Searching,
The Heapsort algorithm is described in the following book,
The Basic Linear Algebra Subprograms (blas) define a set of fundamental operations on vectors and matrices which can be used to create optimized higher-level linear algebra functionality.
The library provides a low-level layer which corresponds directly to the
C-language blas standard, referred to here as “cblas”, and a
higher-level interface for operations on GSL vectors and matrices.
Users who are interested in simple operations on GSL vector and matrix
objects should use the high-level layer, which is declared in the file
gsl_blas.h. This should satisfy the needs of most users. Note
that GSL matrices are implemented using dense-storage so the interface
only includes the corresponding dense-storage blas functions. The full
blas functionality for band-format and packed-format matrices is
available through the low-level cblas interface.
The interface for the gsl_cblas layer is specified in the file
gsl_cblas.h. This interface corresponds to the blas Technical
Forum's draft standard for the C interface to legacy blas
implementations. Users who have access to other conforming cblas
implementations can use these in place of the version provided by the
library. Note that users who have only a Fortran blas library can
use a cblas conformant wrapper to convert it into a cblas
library. A reference cblas wrapper for legacy Fortran
implementations exists as part of the draft cblas standard and can
be obtained from Netlib. The complete set of cblas functions is
listed in an appendix (see GSL CBLAS Library).
There are three levels of blas operations,
Each routine has a name which specifies the operation, the type of matrices involved and their precisions. Some of the most common operations and their names are given below,
The types of matrices are,
Each operation is defined for four precisions,
Thus, for example, the name sgemm stands for “single-precision general matrix-matrix multiply” and zgemm stands for “double-precision complex matrix-matrix multiply”.
GSL provides dense vector and matrix objects, based on the relevant
built-in types. The library provides an interface to the blas
operations which apply to these objects. The interface to this
functionality is given in the file gsl_blas.h.
This function computes the sum \alpha + x^T y for the vectors x and y, returning the result in result.
These functions compute the scalar product x^T y for the vectors x and y, returning the result in result.
These functions compute the complex scalar product x^T y for the vectors x and y, returning the result in result