Added script front-end for primer-design code

[htsworkflow.git] / htswanalysis / MACS / lib / gsl / gsl-1.11 / doc / gsl-ref.info-1

This is gsl-ref.info, produced by makeinfo version 4.8 from
gsl-ref.texi.

INFO-DIR-SECTION Scientific software
START-INFO-DIR-ENTRY
* gsl-ref: (gsl-ref).                   GNU Scientific Library - Reference
END-INFO-DIR-ENTRY

   Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007 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 the freedom to copy and modify
this GNU Manual."

\1f
File: gsl-ref.info,  Node: Top,  Next: Introduction,  Prev: (dir),  Up: (dir)

GSL
***

This file documents the GNU Scientific Library (GSL), a collection of
numerical routines for scientific computing.  It corresponds to release
1.11 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.

   A Japanese translation of this manual is available from the GSL
project homepage thanks to Daisuke Tominaga.

   Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007 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 the freedom to copy and modify
this GNU Manual."


* Menu:

* Introduction::
* Using the library::
* Error Handling::
* Mathematical Functions::
* Complex Numbers::
* Polynomials::
* Special Functions::
* Vectors and Matrices::
* Permutations::
* Combinations::
* Sorting::
* BLAS Support::
* Linear Algebra::
* Eigensystems::
* Fast Fourier Transforms::
* Numerical Integration::
* Random Number Generation::
* Quasi-Random Sequences::
* Random Number Distributions::
* Statistics::
* Histograms::
* N-tuples::
* Monte Carlo Integration::
* Simulated Annealing::
* Ordinary Differential Equations::
* Interpolation::
* Numerical Differentiation::
* Chebyshev Approximations::
* Series Acceleration::
* Wavelet Transforms::
* Discrete Hankel Transforms::
* One dimensional Root-Finding::
* One dimensional Minimization::
* Multidimensional Root-Finding::
* Multidimensional Minimization::
* Least-Squares Fitting::
* Nonlinear Least-Squares Fitting::
* Basis Splines::
* Physical Constants::
* IEEE floating-point arithmetic::
* Debugging Numerical Programs::
* Contributors to GSL::
* Autoconf Macros::
* GSL CBLAS Library::
* Free Software Needs Free Documentation::
* GNU General Public License::
* GNU Free Documentation License::
* Function Index::
* Variable Index::
* Type Index::
* Concept Index::

\1f
File: gsl-ref.info,  Node: Introduction,  Next: Using the library,  Prev: Top,  Up: Top

1 Introduction
**************

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.

* Menu:

* Routines available in GSL::
* GSL is Free Software::
* Obtaining GSL::
* No Warranty::
* Reporting Bugs::
* Further Information::
* Conventions used in this manual::

\1f
File: gsl-ref.info,  Node: Routines available in GSL,  Next: GSL is Free Software,  Up: Introduction

1.1 Routines available in GSL
=============================

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               
     Basis Splines                    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.

\1f
File: gsl-ref.info,  Node: GSL is Free Software,  Next: Obtaining GSL,  Prev: Routines available in GSL,  Up: Introduction

1.2 GSL is Free Software
========================

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
(*note GNU General Public License::).  Further information about this
license is available from the GNU Project webpage `Frequently Asked
Questions about the GNU GPL',

     `http://www.gnu.org/copyleft/gpl-faq.html'

The Free Software Foundation also operates a license consulting service
for commercial users (contact details available from
`http://www.fsf.org/').

\1f
File: gsl-ref.info,  Node: Obtaining GSL,  Next: No Warranty,  Prev: GSL is Free Software,  Up: Introduction

1.3 Obtaining GSL
=================

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,

     `http://www.gnu.org/software/gsl/'

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.

\1f
File: gsl-ref.info,  Node: No Warranty,  Next: Reporting Bugs,  Prev: Obtaining GSL,  Up: Introduction

1.4 No Warranty
===============

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 (*note GNU General
Public License::).

\1f
File: gsl-ref.info,  Node: Reporting Bugs,  Next: Further Information,  Prev: No Warranty,  Up: Introduction

1.5 Reporting Bugs
==================

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:

   * The version number of GSL

   * The hardware and operating system

   * The compiler used, including version number and compilation options

   * A description of the bug behavior

   * A short program which exercises the bug

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.

\1f
File: gsl-ref.info,  Node: Further Information,  Next: Conventions used in this manual,  Prev: Reporting Bugs,  Up: Introduction

1.6 Further Information
=======================

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/'".

\1f
File: gsl-ref.info,  Node: Conventions used in this manual,  Prev: Further Information,  Up: Introduction

1.7 Conventions used in this manual
===================================

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').

\1f
File: gsl-ref.info,  Node: Using the library,  Next: Error Handling,  Prev: Introduction,  Up: Top

2 Using the library
*******************

This chapter describes how to compile programs that use GSL, and
introduces its conventions.

* Menu:

* An Example Program::
* Compiling and Linking::
* Shared Libraries::
* ANSI C Compliance::
* Inline functions::
* Long double::
* Portability functions::
* Alternative optimized functions::
* Support for different numeric types::
* Compatibility with C++::
* Aliasing of arrays::
* Thread-safety::
* Deprecated Functions::
* Code Reuse::

\1f
File: gsl-ref.info,  Node: An Example Program,  Next: Compiling and Linking,  Up: Using the library

2.1 An Example Program
======================

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,(1)

     J0(5) = -1.775967713143382920e-01

The steps needed to compile this program are described in the following
sections.

   ---------- Footnotes ----------

   (1) The last few digits may vary slightly depending on the compiler
and platform used--this is normal.

\1f
File: gsl-ref.info,  Node: Compiling and Linking,  Next: Shared Libraries,  Prev: An Example Program,  Up: Using the library

2.2 Compiling and Linking
=========================

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.

* Menu:

* Linking programs with the library::
* Linking with an alternative BLAS library::

\1f
File: gsl-ref.info,  Node: Linking programs with the library,  Next: Linking with an alternative BLAS library,  Up: Compiling and Linking

2.2.1 Linking programs with the library
---------------------------------------

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.

\1f
File: gsl-ref.info,  Node: Linking with an alternative BLAS library,  Prev: Linking programs with the library,  Up: Compiling and Linking

2.2.2 Linking with an alternative BLAS library
----------------------------------------------

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 *Note BLAS Support::.

\1f
File: gsl-ref.info,  Node: Shared Libraries,  Next: ANSI C Compliance,  Prev: Compiling and Linking,  Up: Using the library

2.3 Shared Libraries
====================

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

\1f
File: gsl-ref.info,  Node: ANSI C Compliance,  Next: Inline functions,  Prev: Shared Libraries,  Up: Using the library

2.4 ANSI C Compliance
=====================

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_'.

\1f
File: gsl-ref.info,  Node: Inline functions,  Next: Long double,  Prev: ANSI C Compliance,  Up: Using the library

2.5 Inline functions
====================

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 *Note Autoconf Macros::.

\1f
File: gsl-ref.info,  Node: Long double,  Next: Portability functions,  Prev: Inline functions,  Up: Using the library

2.6 Long double
===============

In general, the algorithms in the library are written for double
precision only.  The `long double' type is not supported for actual
computation.

   One reason for this choice is that the precision of `long double' is
platform dependent.  The IEEE standard only specifies the minimum
precision of extended precision numbers, while the precision of
`double' is the same on all platforms.

   However, it is sometimes necessary to interact with external data in
long-double format, so the vector and matrix datatypes include
long-double versions.

   It should be noted that 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.

\1f
File: gsl-ref.info,  Node: Portability functions,  Next: Alternative optimized functions,  Prev: Long double,  Up: Using the library

2.7 Portability functions
=========================

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 *Note 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.

\1f
File: gsl-ref.info,  Node: Alternative optimized functions,  Next: Support for different numeric types,  Prev: Portability functions,  Up: Using the library

2.8 Alternative optimized functions
===================================

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.

\1f
File: gsl-ref.info,  Node: Support for different numeric types,  Next: Compatibility with C++,  Prev: Alternative optimized functions,  Up: Using the library

2.9 Support for different numeric types
=======================================

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

\1f
File: gsl-ref.info,  Node: Compatibility with C++,  Next: Aliasing of arrays,  Prev: Support for different numeric types,  Up: Using the library

2.10 Compatibility with C++
===========================

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'.

\1f
File: gsl-ref.info,  Node: Aliasing of arrays,  Next: Thread-safety,  Prev: Compatibility with C++,  Up: Using the library

2.11 Aliasing of arrays
=======================

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.

\1f
File: gsl-ref.info,  Node: Thread-safety,  Next: Deprecated Functions,  Prev: Aliasing of arrays,  Up: Using the library

2.12 Thread-safety
==================

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.

\1f
File: gsl-ref.info,  Node: Deprecated Functions,  Next: Code Reuse,  Prev: Thread-safety,  Up: Using the library

2.13 Deprecated Functions
=========================

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.

\1f
File: gsl-ref.info,  Node: Code Reuse,  Prev: Deprecated Functions,  Up: Using the library

2.14 Code Reuse
===============

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.

\1f
File: gsl-ref.info,  Node: Error Handling,  Next: Mathematical Functions,  Prev: Using the library,  Up: Top

3 Error Handling
****************

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'.

* Menu:

* Error Reporting::
* Error Codes::
* Error Handlers::
* Using GSL error reporting in your own functions::
* Error Reporting Examples::

\1f
File: gsl-ref.info,  Node: Error Reporting,  Next: Error Codes,  Up: Error Handling

3.1 Error Reporting
===================

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.

\1f
File: gsl-ref.info,  Node: Error Codes,  Next: Error Handlers,  Prev: Error Reporting,  Up: Error Handling

3.2 Error 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. Error codes above 1024 are reserved
for applications, and are not used by the library.  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,

 -- Macro: int GSL_EDOM
     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)

 -- Macro: int GSL_ERANGE
     Range error; used by mathematical functions when the result value
     is not representable because of overflow or underflow (like ERANGE
     in the C library)

 -- Macro: int GSL_ENOMEM
     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'.

 -- Macro: int GSL_EINVAL
     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'.

 -- Function: const char * gsl_strerror (const int GSL_ERRNO)
     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 error' for a
     status value of `GSL_ERANGE'.

\1f
File: gsl-ref.info,  Node: Error Handlers,  Next: Using GSL error reporting in your own functions,  Prev: Error Codes,  Up: Error Handling

3.3 Error Handlers
==================

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',

 -- Data Type: gsl_error_handler_t
     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 type
     `void'.  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',

 -- Function: gsl_error_handler_t * gsl_set_error_handler
          (gsl_error_handler_t * NEW_HANDLER)
     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 (`abort' on error) set the error
     handler to `NULL',

          old_handler = gsl_set_error_handler (NULL);

 -- Function: gsl_error_handler_t * gsl_set_error_handler_off ()
     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'.

\1f
File: gsl-ref.info,  Node: Using GSL error reporting in your own functions,  Next: Error Reporting Examples,  Prev: Error Handlers,  Up: Error Handling

3.4 Using GSL error reporting in your own functions
===================================================

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:

 -- Macro: GSL_ERROR (REASON, GSL_ERRNO)
     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);
       }

 -- Macro: GSL_ERROR_VAL (REASON, GSL_ERRNO, VALUE)
     This macro is the same as `GSL_ERROR' but 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);
       }

\1f
File: gsl-ref.info,  Node: Error Reporting Examples,  Prev: Using GSL error reporting in your own functions,  Up: Error Handling

3.5 Examples
============

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, stride, 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.

\1f
File: gsl-ref.info,  Node: Mathematical Functions,  Next: Complex Numbers,  Prev: Error Handling,  Up: Top

4 Mathematical Functions
************************

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'.

* Menu:

* Mathematical Constants::
* Infinities and Not-a-number::
* Elementary Functions::
* Small integer powers::
* Testing the Sign of Numbers::
* Testing for Odd and Even Numbers::
* Maximum and Minimum functions::
* Approximate Comparison of Floating Point Numbers::

\1f
File: gsl-ref.info,  Node: Mathematical Constants,  Next: Infinities and Not-a-number,  Up: Mathematical Functions

4.1 Mathematical Constants
==========================

The library ensures that the standard BSD mathematical constants are
defined. For reference, here is a list of the constants:

`M_E'
     The base of exponentials, e

`M_LOG2E'
     The base-2 logarithm of e, \log_2 (e)

`M_LOG10E'
     The base-10 logarithm of e, \log_10 (e)

`M_SQRT2'
     The square root of two, \sqrt 2

`M_SQRT1_2'
     The square root of one-half, \sqrt{1/2}

`M_SQRT3'
     The square root of three, \sqrt 3

`M_PI'
     The constant pi, \pi

`M_PI_2'
     Pi divided by two, \pi/2

`M_PI_4'
     Pi divided by four, \pi/4

`M_SQRTPI'
     The square root of pi, \sqrt\pi

`M_2_SQRTPI'
     Two divided by the square root of pi, 2/\sqrt\pi

`M_1_PI'
     The reciprocal of pi, 1/\pi

`M_2_PI'
     Twice the reciprocal of pi, 2/\pi

`M_LN10'
     The natural logarithm of ten, \ln(10)

`M_LN2'
     The natural logarithm of two, \ln(2)

`M_LNPI'
     The natural logarithm of pi, \ln(\pi)

`M_EULER'
     Euler's constant, \gamma


\1f
File: gsl-ref.info,  Node: Infinities and Not-a-number,  Next: Elementary Functions,  Prev: Mathematical Constants,  Up: Mathematical Functions

4.2 Infinities and Not-a-number
===============================

 -- Macro: GSL_POSINF
     This macro contains the IEEE representation of positive infinity,
     +\infty. It is computed from the expression `+1.0/0.0'.

 -- Macro: GSL_NEGINF
     This macro contains the IEEE representation of negative infinity,
     -\infty. It is computed from the expression `-1.0/0.0'.

 -- Macro: GSL_NAN
     This macro contains the IEEE representation of the Not-a-Number
     symbol, `NaN'. It is computed from the ratio `0.0/0.0'.

 -- Function: int gsl_isnan (const double X)
     This function returns 1 if X is not-a-number.

 -- Function: int gsl_isinf (const double X)
     This function returns +1 if X is positive infinity, -1 if X is
     negative infinity and 0 otherwise.

 -- Function: int gsl_finite (const double X)
     This function returns 1 if X is a real number, and 0 if it is
     infinite or not-a-number.

\1f
File: gsl-ref.info,  Node: Elementary Functions,  Next: Small integer powers,  Prev: Infinities and Not-a-number,  Up: Mathematical Functions

4.3 Elementary Functions
========================

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
(*note Portability functions::).

 -- Function: double gsl_log1p (const double X)
     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)'.

 -- Function: double gsl_expm1 (const double 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)'.

 -- Function: double gsl_hypot (const double X, const double Y)
     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)'.

 -- Function: double gsl_hypot3 (const double X, const double Y, const
          double Z)
     This function computes the value of \sqrt{x^2 + y^2 + z^2} in a
     way that avoids overflow.

 -- Function: double gsl_acosh (const double X)
     This function computes the value of \arccosh(x). It provides an
     alternative to the standard math function `acosh(x)'.

 -- Function: double gsl_asinh (const double X)
     This function computes the value of \arcsinh(x). It provides an
     alternative to the standard math function `asinh(x)'.

 -- Function: double gsl_atanh (const double X)
     This function computes the value of \arctanh(x). It provides an
     alternative to the standard math function `atanh(x)'.

 -- Function: double gsl_ldexp (double X, int E)
     This function computes the value of x * 2^e. It provides an
     alternative to the standard math function `ldexp(x,e)'.

 -- Function: double gsl_frexp (double X, int * 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)'.

\1f
File: gsl-ref.info,  Node: Small integer powers,  Next: Testing the Sign of Numbers,  Prev: Elementary Functions,  Up: Mathematical Functions

4.4 Small integer powers
========================

A common complaint about the standard C library is its lack of a
function for calculating (small) integer powers.  GSL provides some
simple functions to fill this gap.  For reasons of efficiency, these
functions do not check for overflow or underflow conditions.

 -- Function: double gsl_pow_int (double X, int N)
     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'.

 -- Function: double gsl_pow_2 (const double X)
 -- Function: double gsl_pow_3 (const double X)
 -- Function: double gsl_pow_4 (const double X)
 -- Function: double gsl_pow_5 (const double X)
 -- Function: double gsl_pow_6 (const double X)
 -- Function: double gsl_pow_7 (const double X)
 -- Function: double gsl_pow_8 (const double X)
 -- Function: double gsl_pow_9 (const double X)
     These functions can be used to compute small integer powers x^2,
     x^3, etc. efficiently. The functions will be inlined when
     `HAVE_INLINE' is defined, 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 */

\1f
File: gsl-ref.info,  Node: Testing the Sign of Numbers,  Next: Testing for Odd and Even Numbers,  Prev: Small integer powers,  Up: Mathematical Functions

4.5 Testing the Sign of Numbers
===============================

 -- Macro: GSL_SIGN (x)
     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).

\1f
File: gsl-ref.info,  Node: Testing for Odd and Even Numbers,  Next: Maximum and Minimum functions,  Prev: Testing the Sign of Numbers,  Up: Mathematical Functions

4.6 Testing for Odd and Even Numbers
====================================

 -- Macro: GSL_IS_ODD (n)
     This macro evaluates to 1 if N is odd and 0 if N is even. The
     argument N must be of integer type.

 -- Macro: GSL_IS_EVEN (n)
     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.

\1f
File: gsl-ref.info,  Node: Maximum and Minimum functions,  Next: Approximate Comparison of Floating Point Numbers,  Prev: Testing for Odd and Even Numbers,  Up: Mathematical Functions

4.7 Maximum and Minimum functions
=================================

 -- Macro: GSL_MAX (a, b)
     This macro returns the maximum of A and B. It is defined as `((a)
     > (b) ? (a):(b))'.

 -- Macro: GSL_MIN (a, b)
     This macro returns the minimum of A and B. It is defined as `((a)
     < (b) ? (a):(b))'.

 -- Function: extern inline double GSL_MAX_DBL (double A, double 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_MAX' will be automatically substituted.

 -- Function: extern inline double GSL_MIN_DBL (double A, double B)
     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_MIN' will be automatically substituted.

 -- Function: extern inline int GSL_MAX_INT (int A, int B)
 -- Function: extern inline int GSL_MIN_INT (int A, int B)
     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_MAX' or `GSL_MIN' will
     be automatically substituted.

 -- Function: extern inline long double GSL_MAX_LDBL (long double A,
          long double B)
 -- Function: extern inline long double GSL_MIN_LDBL (long double A,
          long double B)
     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_MAX' or `GSL_MIN' will
     be automatically substituted.

\1f
File: gsl-ref.info,  Node: Approximate Comparison of Floating Point Numbers,  Prev: Maximum and Minimum functions,  Up: Mathematical Functions

4.8 Approximate Comparison of Floating Point Numbers
====================================================

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).

 -- Function: int gsl_fcmp (double X, double Y, double EPSILON)
     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.

     Note that x and y are compared to relative accuracy, so this
     function is not suitable for testing whether a value is
     approximately zero.

     The implementation is based on the package `fcmp' by T.C. Belding.

\1f
File: gsl-ref.info,  Node: Complex Numbers,  Next: Polynomials,  Prev: Mathematical Functions,  Up: Top

5 Complex Numbers
*****************

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'.

* Menu:

* Complex numbers::
* Properties of complex numbers::
* Complex arithmetic operators::
* Elementary Complex Functions::
* Complex Trigonometric Functions::
* Inverse Complex Trigonometric Functions::
* Complex Hyperbolic Functions::
* Inverse Complex Hyperbolic Functions::
* Complex Number References and Further Reading::

   ---------- Footnotes ----------

   (1) Note that the first edition uses different definitions.

\1f
File: gsl-ref.info,  Node: Complex numbers,  Next: Properties of complex numbers,  Up: Complex Numbers

5.1 Complex numbers
===================

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.

 -- Function: gsl_complex gsl_complex_rect (double X, double Y)
     This function uses the rectangular cartesian components (X,Y) to
     return the complex number z = x + i y.  An inline version of this
     function is used when `HAVE_INLINE' is defined.

 -- Function: gsl_complex gsl_complex_polar (double R, double THETA)
     This function returns the complex number z = r \exp(i \theta) = r
     (\cos(\theta) + i \sin(\theta)) from the polar representation
     (R,THETA).

 -- Macro: GSL_REAL (Z)
 -- Macro: GSL_IMAG (Z)
     These macros return the real and imaginary parts of the complex
     number Z.

 -- Macro: GSL_SET_COMPLEX (ZP, X, Y)
     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.

 -- Macro: GSL_SET_REAL (ZP,X)
 -- Macro: GSL_SET_IMAG (ZP,Y)
     These macros allow the real and imaginary parts of the complex
     number pointed to by ZP to be set independently.

\1f
File: gsl-ref.info,  Node: Properties of complex numbers,  Next: Complex arithmetic operators,  Prev: Complex numbers,  Up: Complex Numbers

5.2 Properties of complex numbers
=================================

 -- Function: double gsl_complex_arg (gsl_complex Z)
     This function returns the argument of the complex number Z,
     \arg(z), where -\pi < \arg(z) <= \pi.

 -- Function: double gsl_complex_abs (gsl_complex Z)
     This function returns the magnitude of the complex number Z, |z|.

 -- Function: double gsl_complex_abs2 (gsl_complex Z)
     This function returns the squared magnitude of the complex number
     Z, |z|^2.

 -- Function: double gsl_complex_logabs (gsl_complex Z)
     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.

\1f
File: gsl-ref.info,  Node: Complex arithmetic operators,  Next: Elementary Complex Functions,  Prev: Properties of complex numbers,  Up: Complex Numbers

5.3 Complex arithmetic operators
================================

 -- Function: gsl_complex gsl_complex_add (gsl_complex A, gsl_complex B)
     This function returns the sum of the complex numbers A and B,
     z=a+b.

 -- Function: gsl_complex gsl_complex_sub (gsl_complex A, gsl_complex B)
     This function returns the difference of the complex numbers A and
     B, z=a-b.

 -- Function: gsl_complex gsl_complex_mul (gsl_complex A, gsl_complex B)
     This function returns the product of the complex numbers A and B,
     z=ab.

 -- Function: gsl_complex gsl_complex_div (gsl_complex A, gsl_complex B)
     This function returns the quotient of the complex numbers A and B,
     z=a/b.

 -- Function: gsl_complex gsl_complex_add_real (gsl_complex A, double X)
     This function returns the sum of the complex number A and the real
     number X, z=a+x.

 -- Function: gsl_complex gsl_complex_sub_real (gsl_complex A, double X)
     This function returns the difference of the complex number A and
     the real number X, z=a-x.

 -- Function: gsl_complex gsl_complex_mul_real (gsl_complex A, double X)
     This function returns the product of the complex number A and the
     real number X, z=ax.

 -- Function: gsl_complex gsl_complex_div_real (gsl_complex A, double X)
     This function returns the quotient of the complex number A and the
     real number X, z=a/x.

 -- Function: gsl_complex gsl_complex_add_imag (gsl_complex A, double Y)
     This function returns the sum of the complex number A and the
     imaginary number iY, z=a+iy.

 -- Function: gsl_complex gsl_complex_sub_imag (gsl_complex A, double Y)
     This function returns the difference of the complex number A and
     the imaginary number iY, z=a-iy.

 -- Function: gsl_complex gsl_complex_mul_imag (gsl_complex A, double Y)
     This function returns the product of the complex number A and the
     imaginary number iY, z=a*(iy).

 -- Function: gsl_complex gsl_complex_div_imag (gsl_complex A, double Y)
     This function returns the quotient of the complex number A and the
     imaginary number iY, z=a/(iy).

 -- Function: gsl_complex gsl_complex_conjugate (gsl_complex Z)
     This function returns the complex conjugate of the complex number
     Z, z^* = x - i y.

 -- Function: gsl_complex gsl_complex_inverse (gsl_complex Z)
     This function returns the inverse, or reciprocal, of the complex
     number Z, 1/z = (x - i y)/(x^2 + y^2).

 -- Function: gsl_complex gsl_complex_negative (gsl_complex Z)
     This function returns the negative of the complex number Z, -z =
     (-x) + i(-y).

\1f
File: gsl-ref.info,  Node: Elementary Complex Functions,  Next: Complex Trigonometric Functions,  Prev: Complex arithmetic operators,  Up: Complex Numbers

5.4 Elementary Complex Functions
================================

 -- Function: gsl_complex gsl_complex_sqrt (gsl_complex Z)
     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.

 -- Function: gsl_complex gsl_complex_sqrt_real (double X)
     This function returns the complex square root of the real number
     X, where X may be negative.

 -- Function: gsl_complex gsl_complex_pow (gsl_complex Z, gsl_complex A)
     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.

 -- Function: gsl_complex gsl_complex_pow_real (gsl_complex Z, double X)
     This function returns the complex number Z raised to the real
     power X, z^x.

 -- Function: gsl_complex gsl_complex_exp (gsl_complex Z)
     This function returns the complex exponential of the complex number
     Z, \exp(z).

 -- Function: gsl_complex gsl_complex_log (gsl_complex 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.

 -- Function: gsl_complex gsl_complex_log10 (gsl_complex Z)
     This function returns the complex base-10 logarithm of the complex
     number Z, \log_10 (z).

 -- Function: gsl_complex gsl_complex_log_b (gsl_complex Z, gsl_complex
          B)
     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).

\1f
File: gsl-ref.info,  Node: Complex Trigonometric Functions,  Next: Inverse Complex Trigonometric Functions,  Prev: Elementary Complex Functions,  Up: Complex Numbers

5.5 Complex Trigonometric Functions
===================================

 -- Function: gsl_complex gsl_complex_sin (gsl_complex Z)
     This function returns the complex sine of the complex number Z,
     \sin(z) = (\exp(iz) - \exp(-iz))/(2i).

 -- Function: gsl_complex gsl_complex_cos (gsl_complex Z)
     This function returns the complex cosine of the complex number Z,
     \cos(z) = (\exp(iz) + \exp(-iz))/2.

 -- Function: gsl_complex gsl_complex_tan (gsl_complex Z)
     This function returns the complex tangent of the complex number Z,
     \tan(z) = \sin(z)/\cos(z).

 -- Function: gsl_complex gsl_complex_sec (gsl_complex Z)
     This function returns the complex secant of the complex number Z,
     \sec(z) = 1/\cos(z).

 -- Function: gsl_complex gsl_complex_csc (gsl_complex Z)
     This function returns the complex cosecant of the complex number Z,
     \csc(z) = 1/\sin(z).

 -- Function: gsl_complex gsl_complex_cot (gsl_complex Z)
     This function returns the complex cotangent of the complex number
     Z, \cot(z) = 1/\tan(z).

\1f
File: gsl-ref.info,  Node: Inverse Complex Trigonometric Functions,  Next: Complex Hyperbolic Functions,  Prev: Complex Trigonometric Functions,  Up: Complex Numbers

5.6 Inverse Complex Trigonometric Functions
===========================================

 -- Function: gsl_complex gsl_complex_arcsin (gsl_complex 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.

 -- Function: gsl_complex gsl_complex_arcsin_real (double Z)
     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.

 -- Function: gsl_complex gsl_complex_arccos (gsl_complex Z)
     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.

 -- Function: gsl_complex gsl_complex_arccos_real (double Z)
     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.

 -- Function: gsl_complex gsl_complex_arctan (gsl_complex Z)
     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.

 -- Function: gsl_complex gsl_complex_arcsec (gsl_complex Z)
     This function returns the complex arcsecant of the complex number
     Z, \arcsec(z) = \arccos(1/z).

 -- Function: gsl_complex gsl_complex_arcsec_real (double Z)
     This function returns the complex arcsecant of the real number Z,
     \arcsec(z) = \arccos(1/z).

 -- Function: gsl_complex gsl_complex_arccsc (gsl_complex Z)
     This function returns the complex arccosecant of the complex
     number Z, \arccsc(z) = \arcsin(1/z).

 -- Function: gsl_complex gsl_complex_arccsc_real (double Z)
     This function returns the complex arccosecant of the real number Z,
     \arccsc(z) = \arcsin(1/z).

 -- Function: gsl_complex gsl_complex_arccot (gsl_complex Z)
     This function returns the complex arccotangent of the complex
     number Z, \arccot(z) = \arctan(1/z).

\1f
File: gsl-ref.info,  Node: Complex Hyperbolic Functions,  Next: Inverse Complex Hyperbolic Functions,  Prev: Inverse Complex Trigonometric Functions,  Up: Complex Numbers

5.7 Complex Hyperbolic Functions
================================

 -- Function: gsl_complex gsl_complex_sinh (gsl_complex Z)
     This function returns the complex hyperbolic sine of the complex
     number Z, \sinh(z) = (\exp(z) - \exp(-z))/2.

 -- Function: gsl_complex gsl_complex_cosh (gsl_complex Z)
     This function returns the complex hyperbolic cosine of the complex
     number Z, \cosh(z) = (\exp(z) + \exp(-z))/2.

 -- Function: gsl_complex gsl_complex_tanh (gsl_complex Z)
     This function returns the complex hyperbolic tangent of the
     complex number Z, \tanh(z) = \sinh(z)/\cosh(z).

 -- Function: gsl_complex gsl_complex_sech (gsl_complex Z)
     This function returns the complex hyperbolic secant of the complex
     number Z, \sech(z) = 1/\cosh(z).

 -- Function: gsl_complex gsl_complex_csch (gsl_complex Z)
     This function returns the complex hyperbolic cosecant of the
     complex number Z, \csch(z) = 1/\sinh(z).

 -- Function: gsl_complex gsl_complex_coth (gsl_complex Z)
     This function returns the complex hyperbolic cotangent of the
     complex number Z, \coth(z) = 1/\tanh(z).

\1f
File: gsl-ref.info,  Node: Inverse Complex Hyperbolic Functions,  Next: Complex Number References and Further Reading,  Prev: Complex Hyperbolic Functions,  Up: Complex Numbers

5.8 Inverse Complex Hyperbolic Functions
========================================

 -- Function: gsl_complex gsl_complex_arcsinh (gsl_complex 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.

 -- Function: gsl_complex gsl_complex_arccosh (gsl_complex Z)
     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.  Note that in this case we use the negative
     square root in formula 4.6.21 of Abramowitz & Stegun giving
     \arccosh(z)=\log(z-\sqrt{z^2-1}).

 -- Function: gsl_complex gsl_complex_arccosh_real (double Z)
     This function returns the complex hyperbolic arccosine of the real
     number Z, \arccosh(z).

 -- Function: gsl_complex gsl_complex_arctanh (gsl_complex 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.

 -- Function: gsl_complex gsl_complex_arctanh_real (double Z)
     This function returns the complex hyperbolic arctangent of the real
     number Z, \arctanh(z).

 -- Function: gsl_complex gsl_complex_arcsech (gsl_complex Z)
     This function returns the complex hyperbolic arcsecant of the
     complex number Z, \arcsech(z) = \arccosh(1/z).

 -- Function: gsl_complex gsl_complex_arccsch (gsl_complex Z)
     This function returns the complex hyperbolic arccosecant of the
     complex number Z, \arccsch(z) = \arcsin(1/z).

 -- Function: gsl_complex gsl_complex_arccoth (gsl_complex Z)
     This function returns the complex hyperbolic arccotangent of the
     complex number Z, \arccoth(z) = \arctanh(1/z).

\1f
File: gsl-ref.info,  Node: Complex Number References and Further Reading,  Prev: Inverse Complex Hyperbolic Functions,  Up: Complex Numbers

5.9 References and Further Reading
==================================

The implementations of the elementary and trigonometric functions are
based on the following papers,

     T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang,
     "Implementing Complex Elementary Functions Using Exception
     Handling", `ACM Transactions on Mathematical Software', Volume 20
     (1994), pp 215-244, Corrigenda, p553

     T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang,
     "Implementing the complex arcsin and arccosine functions using
     exception handling", `ACM Transactions on Mathematical Software',
     Volume 23 (1997) pp 299-335

The general formulas and details of branch cuts can be found in the
following books,

     Abramowitz and Stegun, `Handbook of Mathematical Functions',
     "Circular Functions in Terms of Real and Imaginary Parts", Formulas
     4.3.55-58, "Inverse Circular Functions in Terms of Real and
     Imaginary Parts", Formulas 4.4.37-39, "Hyperbolic Functions in
     Terms of Real and Imaginary Parts", Formulas 4.5.49-52, "Inverse
     Hyperbolic Functions--relation to Inverse Circular Functions",
     Formulas 4.6.14-19.

     Dave Gillespie, `Calc Manual', Free Software Foundation, ISBN
     1-882114-18-3

\1f
File: gsl-ref.info,  Node: Polynomials,  Next: Special Functions,  Prev: Complex Numbers,  Up: Top

6 Polynomials
*************

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'.

* Menu:

* Polynomial Evaluation::
* Divided Difference Representation of Polynomials::
* Quadratic Equations::
* Cubic Equations::
* General Polynomial Equations::
* Roots of Polynomials Examples::
* Roots of Polynomials References and Further Reading::

\1f
File: gsl-ref.info,  Node: Polynomial Evaluation,  Next: Divided Difference Representation of Polynomials,  Up: Polynomials

6.1 Polynomial Evaluation
=========================

The functions described here evaluate the polynomial c[0] + c[1] x +
c[2] x^2 + \dots + c[len-1] x^{len-1} using Horner's method for
stability. Inline versions of these functions are used when
`HAVE_INLINE' is defined.

 -- Function: double gsl_poly_eval (const double C[], const int LEN,
          const double X)
     This function evaluates a polynomial with real coefficients for
     the real variable X.

 -- Function: gsl_complex gsl_poly_complex_eval (const double C[],
          const int LEN, const gsl_complex Z)
     This function evaluates a polynomial with real coefficients for
     the complex variable Z.

 -- Function: gsl_complex gsl_complex_poly_complex_eval (const
          gsl_complex C[], const int LEN, const gsl_complex Z)
     This function evaluates a polynomial with complex coefficients for
     the complex variable Z.

\1f
File: gsl-ref.info,  Node: Divided Difference Representation of Polynomials,  Next: Quadratic Equations,  Prev: Polynomial Evaluation,  Up: Polynomials

6.2 Divided Difference Representation of Polynomials
====================================================

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.

 -- Function: int gsl_poly_dd_init (double DD[], const double XA[],
          const double YA[], size_t SIZE)
     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.

 -- Function: double gsl_poly_dd_eval (const double DD[], const double
          XA[], const size_t SIZE, const double X)
     This function evaluates the polynomial stored in
     divided-difference form in the arrays DD and XA of length SIZE at
     the point X.  An inline version of this function is used when
     `HAVE_INLINE' is defined.

 -- Function: int gsl_poly_dd_taylor (double C[], double XP, const
          double DD[], const double XA[], size_t SIZE, double W[])
     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.

\1f
File: gsl-ref.info,  Node: Quadratic Equations,  Next: Cubic Equations,  Prev: Divided Difference Representation of Polynomials,  Up: Polynomials

6.3 Quadratic Equations
=======================

 -- Function: int gsl_poly_solve_quadratic (double A, double B, double
          C, double * X0, double * X1)
     This function finds the real roots of the quadratic equation,

          a x^2 + b x + c = 0

     The 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.


 -- Function: int gsl_poly_complex_solve_quadratic (double A, double B,
          double C, gsl_complex * Z0, gsl_complex * Z1)
     This function finds the complex roots of the quadratic equation,

          a z^2 + b z + c = 0

     The 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.


\1f
File: gsl-ref.info,  Node: Cubic Equations,  Next: General Polynomial Equations,  Prev: Quadratic Equations,  Up: Polynomials

6.4 Cubic Equations
===================

 -- Function: int gsl_poly_solve_cubic (double A, double B, double C,
          double * X0, double * X1, double * X2)
     This function finds the real roots of the cubic equation,

          x^3 + a x^2 + b x + c = 0

     with 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.  As in the quadratic case,
     finite precision may cause equal or closely-spaced real roots to
     move off the real axis into the complex plane, leading to a
     discrete change in the number of real roots.

 -- Function: int gsl_poly_complex_solve_cubic (double A, double B,
          double C, gsl_complex * Z0, gsl_complex * Z1, gsl_complex *
          Z2)
     This function finds the complex roots of the cubic equation,

          z^3 + a z^2 + b z + c = 0

     The 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.


\1f
File: gsl-ref.info,  Node: General Polynomial Equations,  Next: Roots of Polynomials Examples,  Prev: Cubic Equations,  Up: Polynomials

6.5 General Polynomial Equations
================================

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.

 -- Function: gsl_poly_complex_workspace *
gsl_poly_complex_workspace_alloc (size_t N)
     This function allocates space for a `gsl_poly_complex_workspace'
     struct and a workspace suitable for solving a polynomial with N
     coefficients using the routine `gsl_poly_complex_solve'.

     The function returns a pointer to the newly allocated
     `gsl_poly_complex_workspace' if no errors were detected, and a null
     pointer in the case of error.

 -- Function: void gsl_poly_complex_workspace_free
          (gsl_poly_complex_workspace * W)
     This function frees all the memory associated with the workspace W.

 -- Function: int gsl_poly_complex_solve (const double * A, size_t N,
          gsl_poly_complex_workspace * W, gsl_complex_packed_ptr Z)
     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_SUCCESS' if all the roots are found. If
     the QR reduction does not converge, the error handler is invoked
     with an error code of `GSL_EFAILED'.  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).

\1f
File: gsl-ref.info,  Node: Roots of Polynomials Examples,  Next: Roots of Polynomials References and Further Reading,  Prev: General Polynomial Equations,  Up: Polynomials

6.6 Examples
============

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).

\1f
File: gsl-ref.info,  Node: Roots of Polynomials References and Further Reading,  Prev: Roots of Polynomials Examples,  Up: Polynomials

6.7 References and Further Reading
==================================

The balanced-QR method and its error analysis are described in the
following papers,

     R.S. Martin, G. Peters and J.H. Wilkinson, "The QR Algorithm for
     Real Hessenberg Matrices", `Numerische Mathematik', 14 (1970),
     219-231.

     B.N. Parlett and C. Reinsch, "Balancing a Matrix for Calculation of
     Eigenvalues and Eigenvectors", `Numerische Mathematik', 13 (1969),
     293-304.

     A. Edelman and H. Murakami, "Polynomial roots from companion matrix
     eigenvalues", `Mathematics of Computation', Vol. 64, No. 210
     (1995), 763-776.

The formulas for divided differences are given in Abramowitz and Stegun,

     Abramowitz and Stegun, `Handbook of Mathematical Functions',
     Sections 25.1.4 and 25.2.26.

\1f
File: gsl-ref.info,  Node: Special Functions,  Next: Vectors and Matrices,  Prev: Polynomials,  Up: Top

7 Special Functions
*******************

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'.

* Menu:

* Special Function Usage::
* The gsl_sf_result struct::
* Special Function Modes::
* Airy Functions and Derivatives::
* Bessel Functions::
* Clausen Functions::
* Coulomb Functions::
* Coupling Coefficients::
* Dawson Function::
* Debye Functions::
* Dilogarithm::
* Elementary Operations::
* Elliptic Integrals::
* Elliptic Functions (Jacobi)::
* Error Functions::
* Exponential Functions::
* Exponential Integrals::
* Fermi-Dirac Function::
* Gamma and Beta Functions::
* Gegenbauer Functions::
* Hypergeometric Functions::
* Laguerre Functions::
* Lambert W Functions::
* Legendre Functions and Spherical Harmonics::
* Logarithm and Related Functions::
* Mathieu Functions::
* Power Function::
* Psi (Digamma) Function::
* Synchrotron Functions::
* Transport Functions::
* Trigonometric Functions::
* Zeta Functions::
* Special Functions Examples::
* Special Functions References and Further Reading::

\1f
File: gsl-ref.info,  Node: Special Function Usage,  Next: The gsl_sf_result struct,  Up: Special Functions

7.1 Usage
=========

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'.

\1f
File: gsl-ref.info,  Node: The gsl_sf_result struct,  Next: Special Function Modes,  Prev: Special Function Usage,  Up: Special Functions

7.2 The gsl_sf_result struct
============================

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;

\1f
File: gsl-ref.info,  Node: Special Function Modes,  Next: Airy Functions and Derivatives,  Prev: The gsl_sf_result struct,  Up: Special Functions

7.3 Modes
=========

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_DOUBLE'
     Double-precision, a relative accuracy of approximately 2 * 10^-16.

`GSL_PREC_SINGLE'
     Single-precision, a relative accuracy of approximately 10^-7.

`GSL_PREC_APPROX'
     Approximate values, a relative accuracy of approximately 5 * 10^-4.

The approximate mode provides the fastest evaluation at the lowest
accuracy.

\1f
File: gsl-ref.info,  Node: Airy Functions and Derivatives,  Next: Bessel Functions,  Prev: Special Function Modes,  Up: Special Functions

7.4 Airy Functions and Derivatives
==================================

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'.

* Menu:

* Airy Functions::
* Derivatives of Airy Functions::
* Zeros of Airy Functions::
* Zeros of Derivatives of Airy Functions::

\1f
File: gsl-ref.info,  Node: Airy Functions,  Next: Derivatives of Airy Functions,  Up: Airy Functions and Derivatives

7.4.1 Airy Functions
--------------------

 -- Function: double gsl_sf_airy_Ai (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Ai_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the Airy function Ai(x) with an accuracy
     specified by MODE.

 -- Function: double gsl_sf_airy_Bi (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Bi_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the Airy function Bi(x) with an accuracy
     specified by MODE.

 -- Function: double gsl_sf_airy_Ai_scaled (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Ai_scaled_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_airy_Bi_scaled (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Bi_scaled_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Derivatives of Airy Functions,  Next: Zeros of Airy Functions,  Prev: Airy Functions,  Up: Airy Functions and Derivatives

7.4.2 Derivatives of Airy Functions
-----------------------------------

 -- Function: double gsl_sf_airy_Ai_deriv (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Ai_deriv_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the Airy function derivative Ai'(x) with an
     accuracy specified by MODE.

 -- Function: double gsl_sf_airy_Bi_deriv (double X, gsl_mode_t MODE)
 -- Function: int gsl_sf_airy_Bi_deriv_e (double X, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the Airy function derivative Bi'(x) with an
     accuracy specified by MODE.

 -- Function: double gsl_sf_airy_Ai_deriv_scaled (double X, gsl_mode_t
          MODE)
 -- Function: int gsl_sf_airy_Ai_deriv_scaled_e (double X, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_airy_Bi_deriv_scaled (double X, gsl_mode_t
          MODE)
 -- Function: int gsl_sf_airy_Bi_deriv_scaled_e (double X, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Zeros of Airy Functions,  Next: Zeros of Derivatives of Airy Functions,  Prev: Derivatives of Airy Functions,  Up: Airy Functions and Derivatives

7.4.3 Zeros of Airy Functions
-----------------------------

 -- Function: double gsl_sf_airy_zero_Ai (unsigned int S)
 -- Function: int gsl_sf_airy_zero_Ai_e (unsigned int S, gsl_sf_result
          * RESULT)
     These routines compute the location of the S-th zero of the Airy
     function Ai(x).

 -- Function: double gsl_sf_airy_zero_Bi (unsigned int S)
 -- Function: int gsl_sf_airy_zero_Bi_e (unsigned int S, gsl_sf_result
          * RESULT)
     These routines compute the location of the S-th zero of the Airy
     function Bi(x).

\1f
File: gsl-ref.info,  Node: Zeros of Derivatives of Airy Functions,  Prev: Zeros of Airy Functions,  Up: Airy Functions and Derivatives

7.4.4 Zeros of Derivatives of Airy Functions
--------------------------------------------

 -- Function: double gsl_sf_airy_zero_Ai_deriv (unsigned int S)
 -- Function: int gsl_sf_airy_zero_Ai_deriv_e (unsigned int S,
          gsl_sf_result * RESULT)
     These routines compute the location of the S-th zero of the Airy
     function derivative Ai'(x).

 -- Function: double gsl_sf_airy_zero_Bi_deriv (unsigned int S)
 -- Function: int gsl_sf_airy_zero_Bi_deriv_e (unsigned int S,
          gsl_sf_result * RESULT)
     These routines compute the location of the S-th zero of the Airy
     function derivative Bi'(x).

\1f
File: gsl-ref.info,  Node: Bessel Functions,  Next: Clausen Functions,  Prev: Airy Functions and Derivatives,  Up: Special Functions

7.5 Bessel Functions
====================

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'.

* Menu:

* Regular Cylindrical Bessel Functions::
* Irregular Cylindrical Bessel Functions::
* Regular Modified Cylindrical Bessel Functions::
* Irregular Modified Cylindrical Bessel Functions::
* Regular Spherical Bessel Functions::
* Irregular Spherical Bessel Functions::
* Regular Modified Spherical Bessel Functions::
* Irregular Modified Spherical Bessel Functions::
* Regular Bessel Function - Fractional Order::
* Irregular Bessel Functions - Fractional Order::
* Regular Modified Bessel Functions - Fractional Order::
* Irregular Modified Bessel Functions - Fractional Order::
* Zeros of Regular Bessel Functions::

\1f
File: gsl-ref.info,  Node: Regular Cylindrical Bessel Functions,  Next: Irregular Cylindrical Bessel Functions,  Up: Bessel Functions

7.5.1 Regular Cylindrical Bessel Functions
------------------------------------------

 -- Function: double gsl_sf_bessel_J0 (double X)
 -- Function: int gsl_sf_bessel_J0_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular cylindrical Bessel function of
     zeroth order, J_0(x).

 -- Function: double gsl_sf_bessel_J1 (double X)
 -- Function: int gsl_sf_bessel_J1_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular cylindrical Bessel function of
     first order, J_1(x).

 -- Function: double gsl_sf_bessel_Jn (int N, double X)
 -- Function: int gsl_sf_bessel_Jn_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the regular cylindrical Bessel function of
     order N, J_n(x).

 -- Function: int gsl_sf_bessel_Jn_array (int NMIN, int NMAX, double X,
          double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Irregular Cylindrical Bessel Functions,  Next: Regular Modified Cylindrical Bessel Functions,  Prev: Regular Cylindrical Bessel Functions,  Up: Bessel Functions

7.5.2 Irregular Cylindrical Bessel Functions
--------------------------------------------

 -- Function: double gsl_sf_bessel_Y0 (double X)
 -- Function: int gsl_sf_bessel_Y0_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular cylindrical Bessel function
     of zeroth order, Y_0(x), for x>0.

 -- Function: double gsl_sf_bessel_Y1 (double X)
 -- Function: int gsl_sf_bessel_Y1_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular cylindrical Bessel function
     of first order, Y_1(x), for x>0.

 -- Function: double gsl_sf_bessel_Yn (int N,double X)
 -- Function: int gsl_sf_bessel_Yn_e (int N,double X, gsl_sf_result *
          RESULT)
     These routines compute the irregular cylindrical Bessel function of
     order N, Y_n(x), for x>0.

 -- Function: int gsl_sf_bessel_Yn_array (int NMIN, int NMAX, double X,
          double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Regular Modified Cylindrical Bessel Functions,  Next: Irregular Modified Cylindrical Bessel Functions,  Prev: Irregular Cylindrical Bessel Functions,  Up: Bessel Functions

7.5.3 Regular Modified Cylindrical Bessel Functions
---------------------------------------------------

 -- Function: double gsl_sf_bessel_I0 (double X)
 -- Function: int gsl_sf_bessel_I0_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular modified cylindrical Bessel
     function of zeroth order, I_0(x).

 -- Function: double gsl_sf_bessel_I1 (double X)
 -- Function: int gsl_sf_bessel_I1_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular modified cylindrical Bessel
     function of first order, I_1(x).

 -- Function: double gsl_sf_bessel_In (int N, double X)
 -- Function: int gsl_sf_bessel_In_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the regular modified cylindrical Bessel
     function of order N, I_n(x).

 -- Function: int gsl_sf_bessel_In_array (int NMIN, int NMAX, double X,
          double RESULT_ARRAY[])
     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.

 -- Function: double gsl_sf_bessel_I0_scaled (double X)
 -- Function: int gsl_sf_bessel_I0_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled regular modified cylindrical
     Bessel function of zeroth order \exp(-|x|) I_0(x).

 -- Function: double gsl_sf_bessel_I1_scaled (double X)
 -- Function: int gsl_sf_bessel_I1_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled regular modified cylindrical
     Bessel function of first order \exp(-|x|) I_1(x).

 -- Function: double gsl_sf_bessel_In_scaled (int N, double X)
 -- Function: int gsl_sf_bessel_In_scaled_e (int N, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled regular modified cylindrical
     Bessel function of order N, \exp(-|x|) I_n(x)

 -- Function: int gsl_sf_bessel_In_scaled_array (int NMIN, int NMAX,
          double X, double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Irregular Modified Cylindrical Bessel Functions,  Next: Regular Spherical Bessel Functions,  Prev: Regular Modified Cylindrical Bessel Functions,  Up: Bessel Functions

7.5.4 Irregular Modified Cylindrical Bessel Functions
-----------------------------------------------------

 -- Function: double gsl_sf_bessel_K0 (double X)
 -- Function: int gsl_sf_bessel_K0_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular modified cylindrical Bessel
     function of zeroth order, K_0(x), for x > 0.

 -- Function: double gsl_sf_bessel_K1 (double X)
 -- Function: int gsl_sf_bessel_K1_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular modified cylindrical Bessel
     function of first order, K_1(x), for x > 0.

 -- Function: double gsl_sf_bessel_Kn (int N, double X)
 -- Function: int gsl_sf_bessel_Kn_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the irregular modified cylindrical Bessel
     function of order N, K_n(x), for x > 0.

 -- Function: int gsl_sf_bessel_Kn_array (int NMIN, int NMAX, double X,
          double RESULT_ARRAY[])
     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.

 -- Function: double gsl_sf_bessel_K0_scaled (double X)
 -- Function: int gsl_sf_bessel_K0_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled irregular modified cylindrical
     Bessel function of zeroth order \exp(x) K_0(x) for x>0.

 -- Function: double gsl_sf_bessel_K1_scaled (double X)
 -- Function: int gsl_sf_bessel_K1_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled irregular modified cylindrical
     Bessel function of first order \exp(x) K_1(x) for x>0.

 -- Function: double gsl_sf_bessel_Kn_scaled (int N, double X)
 -- Function: int gsl_sf_bessel_Kn_scaled_e (int N, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled irregular modified cylindrical
     Bessel function of order N, \exp(x) K_n(x), for x>0.

 -- Function: int gsl_sf_bessel_Kn_scaled_array (int NMIN, int NMAX,
          double X, double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Regular Spherical Bessel Functions,  Next: Irregular Spherical Bessel Functions,  Prev: Irregular Modified Cylindrical Bessel Functions,  Up: Bessel Functions

7.5.5 Regular Spherical Bessel Functions
----------------------------------------

 -- Function: double gsl_sf_bessel_j0 (double X)
 -- Function: int gsl_sf_bessel_j0_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular spherical Bessel function of
     zeroth order, j_0(x) = \sin(x)/x.

 -- Function: double gsl_sf_bessel_j1 (double X)
 -- Function: int gsl_sf_bessel_j1_e (double X, gsl_sf_result * RESULT)
     These routines compute the regular spherical Bessel function of
     first order, j_1(x) = (\sin(x)/x - \cos(x))/x.

 -- Function: double gsl_sf_bessel_j2 (double X)
 -- Function: int gsl_sf_bessel_j2_e (double X, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_bessel_jl (int L, double X)
 -- Function: int gsl_sf_bessel_jl_e (int L, double X, gsl_sf_result *
          RESULT)
     These routines compute the regular spherical Bessel function of
     order L, j_l(x), for l >= 0 and x >= 0.

 -- Function: int gsl_sf_bessel_jl_array (int LMAX, double X, double
          RESULT_ARRAY[])
     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.

 -- Function: int gsl_sf_bessel_jl_steed_array (int LMAX, double X,
          double * JL_X_ARRAY)
     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.

\1f
File: gsl-ref.info,  Node: Irregular Spherical Bessel Functions,  Next: Regular Modified Spherical Bessel Functions,  Prev: Regular Spherical Bessel Functions,  Up: Bessel Functions

7.5.6 Irregular Spherical Bessel Functions
------------------------------------------

 -- Function: double gsl_sf_bessel_y0 (double X)
 -- Function: int gsl_sf_bessel_y0_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular spherical Bessel function of
     zeroth order, y_0(x) = -\cos(x)/x.

 -- Function: double gsl_sf_bessel_y1 (double X)
 -- Function: int gsl_sf_bessel_y1_e (double X, gsl_sf_result * RESULT)
     These routines compute the irregular spherical Bessel function of
     first order, y_1(x) = -(\cos(x)/x + \sin(x))/x.

 -- Function: double gsl_sf_bessel_y2 (double X)
 -- Function: int gsl_sf_bessel_y2_e (double X, gsl_sf_result * RESULT)
     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).

 -- Function: double gsl_sf_bessel_yl (int L, double X)
 -- Function: int gsl_sf_bessel_yl_e (int L, double X, gsl_sf_result *
          RESULT)
     These routines compute the irregular spherical Bessel function of
     order L, y_l(x), for l >= 0.

 -- Function: int gsl_sf_bessel_yl_array (int LMAX, double X, double
          RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Regular Modified Spherical Bessel Functions,  Next: Irregular Modified Spherical Bessel Functions,  Prev: Irregular Spherical Bessel Functions,  Up: Bessel Functions

7.5.7 Regular Modified Spherical Bessel Functions
-------------------------------------------------

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)

 -- Function: double gsl_sf_bessel_i0_scaled (double X)
 -- Function: int gsl_sf_bessel_i0_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled regular modified spherical Bessel
     function of zeroth order, \exp(-|x|) i_0(x).

 -- Function: double gsl_sf_bessel_i1_scaled (double X)
 -- Function: int gsl_sf_bessel_i1_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled regular modified spherical Bessel
     function of first order, \exp(-|x|) i_1(x).

 -- Function: double gsl_sf_bessel_i2_scaled (double X)
 -- Function: int gsl_sf_bessel_i2_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled regular modified spherical Bessel
     function of second order,  \exp(-|x|) i_2(x)

 -- Function: double gsl_sf_bessel_il_scaled (int L, double X)
 -- Function: int gsl_sf_bessel_il_scaled_e (int L, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled regular modified spherical Bessel
     function of order L,  \exp(-|x|) i_l(x)

 -- Function: int gsl_sf_bessel_il_scaled_array (int LMAX, double X,
          double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Irregular Modified Spherical Bessel Functions,  Next: Regular Bessel Function - Fractional Order,  Prev: Regular Modified Spherical Bessel Functions,  Up: Bessel Functions

7.5.8 Irregular Modified Spherical Bessel Functions
---------------------------------------------------

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).

 -- Function: double gsl_sf_bessel_k0_scaled (double X)
 -- Function: int gsl_sf_bessel_k0_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled irregular modified spherical
     Bessel function of zeroth order, \exp(x) k_0(x), for x>0.

 -- Function: double gsl_sf_bessel_k1_scaled (double X)
 -- Function: int gsl_sf_bessel_k1_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled irregular modified spherical
     Bessel function of first order, \exp(x) k_1(x), for x>0.

 -- Function: double gsl_sf_bessel_k2_scaled (double X)
 -- Function: int gsl_sf_bessel_k2_scaled_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the scaled irregular modified spherical
     Bessel function of second order, \exp(x) k_2(x), for x>0.

 -- Function: double gsl_sf_bessel_kl_scaled (int L, double X)
 -- Function: int gsl_sf_bessel_kl_scaled_e (int L, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled irregular modified spherical
     Bessel function of order L, \exp(x) k_l(x), for x>0.

 -- Function: int gsl_sf_bessel_kl_scaled_array (int LMAX, double X,
          double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Regular Bessel Function - Fractional Order,  Next: Irregular Bessel Functions - Fractional Order,  Prev: Irregular Modified Spherical Bessel Functions,  Up: Bessel Functions

7.5.9 Regular Bessel Function--Fractional Order
-----------------------------------------------

 -- Function: double gsl_sf_bessel_Jnu (double NU, double X)
 -- Function: int gsl_sf_bessel_Jnu_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the regular cylindrical Bessel function of
     fractional order \nu, J_\nu(x).

 -- Function: int gsl_sf_bessel_sequence_Jnu_e (double NU, gsl_mode_t
          MODE, size_t SIZE, double V[])
     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).

\1f
File: gsl-ref.info,  Node: Irregular Bessel Functions - Fractional Order,  Next: Regular Modified Bessel Functions - Fractional Order,  Prev: Regular Bessel Function - Fractional Order,  Up: Bessel Functions

7.5.10 Irregular Bessel Functions--Fractional Order
---------------------------------------------------

 -- Function: double gsl_sf_bessel_Ynu (double NU, double X)
 -- Function: int gsl_sf_bessel_Ynu_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the irregular cylindrical Bessel function of
     fractional order \nu, Y_\nu(x).

\1f
File: gsl-ref.info,  Node: Regular Modified Bessel Functions - Fractional Order,  Next: Irregular Modified Bessel Functions - Fractional Order,  Prev: Irregular Bessel Functions - Fractional Order,  Up: Bessel Functions

7.5.11 Regular Modified Bessel Functions--Fractional Order
----------------------------------------------------------

 -- Function: double gsl_sf_bessel_Inu (double NU, double X)
 -- Function: int gsl_sf_bessel_Inu_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the regular modified Bessel function of
     fractional order \nu, I_\nu(x) for x>0, \nu>0.

 -- Function: double gsl_sf_bessel_Inu_scaled (double NU, double X)
 -- Function: int gsl_sf_bessel_Inu_scaled_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled regular modified Bessel function
     of fractional order \nu, \exp(-|x|)I_\nu(x) for x>0, \nu>0.

\1f
File: gsl-ref.info,  Node: Irregular Modified Bessel Functions - Fractional Order,  Next: Zeros of Regular Bessel Functions,  Prev: Regular Modified Bessel Functions - Fractional Order,  Up: Bessel Functions

7.5.12 Irregular Modified Bessel Functions--Fractional Order
------------------------------------------------------------

 -- Function: double gsl_sf_bessel_Knu (double NU, double X)
 -- Function: int gsl_sf_bessel_Knu_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the irregular modified Bessel function of
     fractional order \nu, K_\nu(x) for x>0, \nu>0.

 -- Function: double gsl_sf_bessel_lnKnu (double NU, double X)
 -- Function: int gsl_sf_bessel_lnKnu_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the logarithm of the irregular modified
     Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0,
     \nu>0.

 -- Function: double gsl_sf_bessel_Knu_scaled (double NU, double X)
 -- Function: int gsl_sf_bessel_Knu_scaled_e (double NU, double X,
          gsl_sf_result * RESULT)
     These routines compute the scaled irregular modified Bessel
     function of fractional order \nu, \exp(+|x|) K_\nu(x) for x>0,
     \nu>0.

\1f
File: gsl-ref.info,  Node: Zeros of Regular Bessel Functions,  Prev: Irregular Modified Bessel Functions - Fractional Order,  Up: Bessel Functions

7.5.13 Zeros of Regular Bessel Functions
----------------------------------------

 -- Function: double gsl_sf_bessel_zero_J0 (unsigned int S)
 -- Function: int gsl_sf_bessel_zero_J0_e (unsigned int S,
          gsl_sf_result * RESULT)
     These routines compute the location of the S-th positive zero of
     the Bessel function J_0(x).

 -- Function: double gsl_sf_bessel_zero_J1 (unsigned int S)
 -- Function: int gsl_sf_bessel_zero_J1_e (unsigned int S,
          gsl_sf_result * RESULT)
     These routines compute the location of the S-th positive zero of
     the Bessel function J_1(x).

 -- Function: double gsl_sf_bessel_zero_Jnu (double NU, unsigned int S)
 -- Function: int gsl_sf_bessel_zero_Jnu_e (double NU, unsigned int S,
          gsl_sf_result * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Clausen Functions,  Next: Coulomb Functions,  Prev: Bessel Functions,  Up: Special Functions

7.6 Clausen Functions
=====================

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'.

 -- Function: double gsl_sf_clausen (double X)
 -- Function: int gsl_sf_clausen_e (double X, gsl_sf_result * RESULT)
     These routines compute the Clausen integral Cl_2(x).

\1f
File: gsl-ref.info,  Node: Coulomb Functions,  Next: Coupling Coefficients,  Prev: Clausen Functions,  Up: Special Functions

7.7 Coulomb Functions
=====================

The prototypes of the Coulomb functions are declared in the header file
`gsl_sf_coulomb.h'.  Both bound state and scattering solutions are
available.

* Menu:

* Normalized Hydrogenic Bound States::
* Coulomb Wave Functions::
* Coulomb Wave Function Normalization Constant::

\1f
File: gsl-ref.info,  Node: Normalized Hydrogenic Bound States,  Next: Coulomb Wave Functions,  Up: Coulomb Functions

7.7.1 Normalized Hydrogenic Bound States
----------------------------------------

 -- Function: double gsl_sf_hydrogenicR_1 (double Z, double R)
 -- Function: int gsl_sf_hydrogenicR_1_e (double Z, double R,
          gsl_sf_result * RESULT)
     These routines compute the lowest-order normalized hydrogenic bound
     state radial wavefunction R_1 := 2Z \sqrt{Z} \exp(-Z r).

 -- Function: double gsl_sf_hydrogenicR (int N, int L, double Z, double
          R)
 -- Function: int gsl_sf_hydrogenicR_e (int N, int L, double Z, double
          R, gsl_sf_result * RESULT)
     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 (*note
     Laguerre Functions::).  The normalization is chosen such that the
     wavefunction \psi is given by \psi(n,l,r) = R_n Y_{lm}.

\1f
File: gsl-ref.info,  Node: Coulomb Wave Functions,  Next: Coulomb Wave Function Normalization Constant,  Prev: Normalized Hydrogenic Bound States,  Up: Coulomb Functions

7.7.2 Coulomb Wave Functions
----------------------------

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)


 -- Function: int gsl_sf_coulomb_wave_FG_e (double ETA, double X,
          double L_F, int K, gsl_sf_result * F, gsl_sf_result * FP,
          gsl_sf_result * G, gsl_sf_result * GP, double * EXP_F, double
          * 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_EOVRFLW' is
     returned and scaling exponents are stored in the modifiable
     parameters EXP_F, EXP_G.

 -- Function: int gsl_sf_coulomb_wave_F_array (double L_MIN, int KMAX,
          double ETA, double X, double FC_ARRAY[], double * F_EXPONENT)
     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.

 -- Function: int gsl_sf_coulomb_wave_FG_array (double L_MIN, int KMAX,
          double ETA, double X, double FC_ARRAY[], double GC_ARRAY[],
          double * F_EXPONENT, double * G_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.

 -- Function: int gsl_sf_coulomb_wave_FGp_array (double L_MIN, int
          KMAX, double ETA, double X, double FC_ARRAY[], double
          FCP_ARRAY[], double GC_ARRAY[], double GCP_ARRAY[], double *
          F_EXPONENT, double * 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.

 -- Function: int gsl_sf_coulomb_wave_sphF_array (double L_MIN, int
          KMAX, double ETA, double X, double FC_ARRAY[], double
          F_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.

\1f
File: gsl-ref.info,  Node: Coulomb Wave Function Normalization Constant,  Prev: Coulomb Wave Functions,  Up: Coulomb Functions

7.7.3 Coulomb Wave Function Normalization Constant
--------------------------------------------------

The Coulomb wave function normalization constant is defined in
Abramowitz 14.1.7.

 -- Function: int gsl_sf_coulomb_CL_e (double L, double ETA,
          gsl_sf_result * RESULT)
     This function computes the Coulomb wave function normalization
     constant C_L(\eta) for L > -1.

 -- Function: int gsl_sf_coulomb_CL_array (double LMIN, int KMAX,
          double ETA, double CL[])
     This function computes the Coulomb wave function normalization
     constant C_L(\eta) for L = Lmin \dots Lmin + kmax, Lmin > -1.

\1f
File: gsl-ref.info,  Node: Coupling Coefficients,  Next: Dawson Function,  Prev: Coulomb Functions,  Up: Special Functions

7.8 Coupling Coefficients
=========================

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'.

* Menu:

* 3-j Symbols::
* 6-j Symbols::
* 9-j Symbols::

\1f
File: gsl-ref.info,  Node: 3-j Symbols,  Next: 6-j Symbols,  Up: Coupling Coefficients

7.8.1 3-j Symbols
-----------------

 -- Function: double gsl_sf_coupling_3j (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC)
 -- Function: int gsl_sf_coupling_3j_e (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC, gsl_sf_result *
          RESULT)
     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.

\1f
File: gsl-ref.info,  Node: 6-j Symbols,  Next: 9-j Symbols,  Prev: 3-j Symbols,  Up: Coupling Coefficients

7.8.2 6-j Symbols
-----------------

 -- Function: double gsl_sf_coupling_6j (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF)
 -- Function: int gsl_sf_coupling_6j_e (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, gsl_sf_result *
          RESULT)
     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.

\1f
File: gsl-ref.info,  Node: 9-j Symbols,  Prev: 6-j Symbols,  Up: Coupling Coefficients

7.8.3 9-j Symbols
-----------------

 -- Function: double gsl_sf_coupling_9j (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int
          TWO_JH, int TWO_JI)
 -- Function: int gsl_sf_coupling_9j_e (int TWO_JA, int TWO_JB, int
          TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int
          TWO_JH, int TWO_JI, gsl_sf_result * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Dawson Function,  Next: Debye Functions,  Prev: Coupling Coefficients,  Up: Special Functions

7.9 Dawson Function
===================

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'.

 -- Function: double gsl_sf_dawson (double X)
 -- Function: int gsl_sf_dawson_e (double X, gsl_sf_result * RESULT)
     These routines compute the value of Dawson's integral for X.

\1f
File: gsl-ref.info,  Node: Debye Functions,  Next: Dilogarithm,  Prev: Dawson Function,  Up: Special Functions

7.10 Debye Functions
====================

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'.

 -- Function: double gsl_sf_debye_1 (double X)
 -- Function: int gsl_sf_debye_1_e (double X, gsl_sf_result * RESULT)
     These routines compute the first-order Debye function D_1(x) =
     (1/x) \int_0^x dt (t/(e^t - 1)).

 -- Function: double gsl_sf_debye_2 (double X)
 -- Function: int gsl_sf_debye_2_e (double X, gsl_sf_result * RESULT)
     These routines compute the second-order Debye function D_2(x) =
     (2/x^2) \int_0^x dt (t^2/(e^t - 1)).

 -- Function: double gsl_sf_debye_3 (double X)
 -- Function: int gsl_sf_debye_3_e (double X, gsl_sf_result * RESULT)
     These routines compute the third-order Debye function D_3(x) =
     (3/x^3) \int_0^x dt (t^3/(e^t - 1)).

 -- Function: double gsl_sf_debye_4 (double X)
 -- Function: int gsl_sf_debye_4_e (double X, gsl_sf_result * RESULT)
     These routines compute the fourth-order Debye function D_4(x) =
     (4/x^4) \int_0^x dt (t^4/(e^t - 1)).

 -- Function: double gsl_sf_debye_5 (double X)
 -- Function: int gsl_sf_debye_5_e (double X, gsl_sf_result * RESULT)
     These routines compute the fifth-order Debye function D_5(x) =
     (5/x^5) \int_0^x dt (t^5/(e^t - 1)).

 -- Function: double gsl_sf_debye_6 (double X)
 -- Function: int gsl_sf_debye_6_e (double X, gsl_sf_result * RESULT)
     These routines compute the sixth-order Debye function D_6(x) =
     (6/x^6) \int_0^x dt (t^6/(e^t - 1)).

\1f
File: gsl-ref.info,  Node: Dilogarithm,  Next: Elementary Operations,  Prev: Debye Functions,  Up: Special Functions

7.11 Dilogarithm
================

The functions described in this section are declared in the header file
`gsl_sf_dilog.h'.

* Menu:

* Real Argument::
* Complex Argument::

\1f
File: gsl-ref.info,  Node: Real Argument,  Next: Complex Argument,  Up: Dilogarithm

7.11.1 Real Argument
--------------------

 -- Function: double gsl_sf_dilog (double X)
 -- Function: int gsl_sf_dilog_e (double X, gsl_sf_result * RESULT)
     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.

     Note that Abramowitz & Stegun refer to the Spence integral
     S(x)=Li_2(1-x) as the dilogarithm rather than Li_2(x).

\1f
File: gsl-ref.info,  Node: Complex Argument,  Prev: Real Argument,  Up: Dilogarithm

7.11.2 Complex Argument
-----------------------

 -- Function: int gsl_sf_complex_dilog_e (double R, double THETA,
          gsl_sf_result * RESULT_RE, gsl_sf_result * RESULT_IM)
     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.

\1f
File: gsl-ref.info,  Node: Elementary Operations,  Next: Elliptic Integrals,  Prev: Dilogarithm,  Up: Special Functions

7.12 Elementary Operations
==========================

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'.

 -- Function: int gsl_sf_multiply_e (double X, double Y, gsl_sf_result
          * RESULT)
     This function multiplies X and Y storing the product and its
     associated error in RESULT.

 -- Function: int gsl_sf_multiply_err_e (double X, double DX, double Y,
          double DY, gsl_sf_result * 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.

\1f
File: gsl-ref.info,  Node: Elliptic Integrals,  Next: Elliptic Functions (Jacobi),  Prev: Elementary Operations,  Up: Special Functions

7.13 Elliptic Integrals
=======================

The functions described in this section are declared in the header file
`gsl_sf_ellint.h'.  Further information about the elliptic integrals
can be found in Abramowitz & Stegun, Chapter 17.

* Menu:

* Definition of Legendre Forms::
* Definition of Carlson Forms::
* Legendre Form of Complete Elliptic Integrals::
* Legendre Form of Incomplete Elliptic Integrals::
* Carlson Forms::

\1f
File: gsl-ref.info,  Node: Definition of Legendre Forms,  Next: Definition of Carlson Forms,  Up: Elliptic Integrals

7.13.1 Definition of Legendre Forms
-----------------------------------

The Legendre forms of elliptic integrals F(\phi,k), E(\phi,k) and
\Pi(\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)))

     Pi(\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).

   The notation used here is based on Carlson, `Numerische Mathematik'
33 (1979) 1 and differs slightly from that used by Abramowitz & Stegun,
where the functions are given in terms of the parameter m = k^2 and n
is replaced by -n.

\1f
File: gsl-ref.info,  Node: Definition of Carlson Forms,  Next: Legendre Form of Complete Elliptic Integrals,  Prev: Definition of Legendre Forms,  Up: Elliptic Integrals

7.13.2 Definition of Carlson Forms
----------------------------------

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)

\1f
File: gsl-ref.info,  Node: Legendre Form of Complete Elliptic Integrals,  Next: Legendre Form of Incomplete Elliptic Integrals,  Prev: Definition of Carlson Forms,  Up: Elliptic Integrals

7.13.3 Legendre Form of Complete Elliptic Integrals
---------------------------------------------------

 -- Function: double gsl_sf_ellint_Kcomp (double K, gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_Kcomp_e (double K, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the complete elliptic integral K(k) to the
     accuracy specified by the mode variable MODE.  Note that
     Abramowitz & Stegun define this function in terms of the parameter
     m = k^2.

 -- Function: double gsl_sf_ellint_Ecomp (double K, gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_Ecomp_e (double K, gsl_mode_t MODE,
          gsl_sf_result * RESULT)
     These routines compute the complete elliptic integral E(k) to the
     accuracy specified by the mode variable MODE.  Note that
     Abramowitz & Stegun define this function in terms of the parameter
     m = k^2.

 -- Function: double gsl_sf_ellint_Pcomp (double K, double N,
          gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_Pcomp_e (double K, double N, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     These routines compute the complete elliptic integral \Pi(k,n) to
     the accuracy specified by the mode variable MODE.  Note that
     Abramowitz & Stegun define this function in terms of the
     parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of
     sign n \to -n.

\1f
File: gsl-ref.info,  Node: Legendre Form of Incomplete Elliptic Integrals,  Next: Carlson Forms,  Prev: Legendre Form of Complete Elliptic Integrals,  Up: Elliptic Integrals

7.13.4 Legendre Form of Incomplete Elliptic Integrals
-----------------------------------------------------

 -- Function: double gsl_sf_ellint_F (double PHI, double K, gsl_mode_t
          MODE)
 -- Function: int gsl_sf_ellint_F_e (double PHI, double K, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral F(\phi,k)
     to the accuracy specified by the mode variable MODE.  Note that
     Abramowitz & Stegun define this function in terms of the parameter
     m = k^2.

 -- Function: double gsl_sf_ellint_E (double PHI, double K, gsl_mode_t
          MODE)
 -- Function: int gsl_sf_ellint_E_e (double PHI, double K, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral E(\phi,k)
     to the accuracy specified by the mode variable MODE.  Note that
     Abramowitz & Stegun define this function in terms of the parameter
     m = k^2.

 -- Function: double gsl_sf_ellint_P (double PHI, double K, double N,
          gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_P_e (double PHI, double K, double N,
          gsl_mode_t MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral
     \Pi(\phi,k,n) to the accuracy specified by the mode variable MODE.
     Note that Abramowitz & Stegun define this function in terms of the
     parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of
     sign n \to -n.

 -- Function: double gsl_sf_ellint_D (double PHI, double K, double N,
          gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_D_e (double PHI, double K, double N,
          gsl_mode_t MODE, gsl_sf_result * RESULT)
     These functions compute the incomplete elliptic integral D(\phi,k)
     which is defined through the Carlson form RD(x,y,z) by the
     following relation,

          D(\phi,k,n) = (1/3)(\sin(\phi))^3 RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).
     The argument N is not used and will be removed in a future release.


\1f
File: gsl-ref.info,  Node: Carlson Forms,  Prev: Legendre Form of Incomplete Elliptic Integrals,  Up: Elliptic Integrals

7.13.5 Carlson Forms
--------------------

 -- Function: double gsl_sf_ellint_RC (double X, double Y, gsl_mode_t
          MODE)
 -- Function: int gsl_sf_ellint_RC_e (double X, double Y, gsl_mode_t
          MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral RC(x,y) to
     the accuracy specified by the mode variable MODE.

 -- Function: double gsl_sf_ellint_RD (double X, double Y, double Z,
          gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_RD_e (double X, double Y, double Z,
          gsl_mode_t MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral RD(x,y,z)
     to the accuracy specified by the mode variable MODE.

 -- Function: double gsl_sf_ellint_RF (double X, double Y, double Z,
          gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_RF_e (double X, double Y, double Z,
          gsl_mode_t MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral RF(x,y,z)
     to the accuracy specified by the mode variable MODE.

 -- Function: double gsl_sf_ellint_RJ (double X, double Y, double Z,
          double P, gsl_mode_t MODE)
 -- Function: int gsl_sf_ellint_RJ_e (double X, double Y, double Z,
          double P, gsl_mode_t MODE, gsl_sf_result * RESULT)
     These routines compute the incomplete elliptic integral RJ(x,y,z,p)
     to the accuracy specified by the mode variable MODE.

\1f
File: gsl-ref.info,  Node: Elliptic Functions (Jacobi),  Next: Error Functions,  Prev: Elliptic Integrals,  Up: Special Functions

7.14 Elliptic Functions (Jacobi)
================================

The Jacobian Elliptic functions are defined in Abramowitz & Stegun,
Chapter 16.  The functions are declared in the header file
`gsl_sf_elljac.h'.

 -- Function: int gsl_sf_elljac_e (double U, double M, double * SN,
          double * CN, double * DN)
     This function computes the Jacobian elliptic functions sn(u|m),
     cn(u|m), dn(u|m) by descending Landen transformations.

\1f
File: gsl-ref.info,  Node: Error Functions,  Next: Exponential Functions,  Prev: Elliptic Functions (Jacobi),  Up: Special Functions

7.15 Error Functions
====================

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'.

* Menu:

* Error Function::
* Complementary Error Function::
* Log Complementary Error Function::
* Probability functions::

\1f
File: gsl-ref.info,  Node: Error Function,  Next: Complementary Error Function,  Up: Error Functions

7.15.1 Error Function
---------------------

 -- Function: double gsl_sf_erf (double X)
 -- Function: int gsl_sf_erf_e (double X, gsl_sf_result * RESULT)
     These routines compute the error function erf(x), where erf(x) =
     (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2).

\1f
File: gsl-ref.info,  Node: Complementary Error Function,  Next: Log Complementary Error Function,  Prev: Error Function,  Up: Error Functions

7.15.2 Complementary Error Function
-----------------------------------

 -- Function: double gsl_sf_erfc (double X)
 -- Function: int gsl_sf_erfc_e (double X, gsl_sf_result * RESULT)
     These routines compute the complementary error function erfc(x) =
     1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2).

\1f
File: gsl-ref.info,  Node: Log Complementary Error Function,  Next: Probability functions,  Prev: Complementary Error Function,  Up: Error Functions

7.15.3 Log Complementary Error Function
---------------------------------------

 -- Function: double gsl_sf_log_erfc (double X)
 -- Function: int gsl_sf_log_erfc_e (double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of the complementary error
     function \log(\erfc(x)).

\1f
File: gsl-ref.info,  Node: Probability functions,  Prev: Log Complementary Error Function,  Up: Error Functions

7.15.4 Probability functions
----------------------------

The probability functions for the Normal or Gaussian distribution are
described in Abramowitz & Stegun, Section 26.2.

 -- Function: double gsl_sf_erf_Z (double X)
 -- Function: int gsl_sf_erf_Z_e (double X, gsl_sf_result * RESULT)
     These routines compute the Gaussian probability density function
     Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2).

 -- Function: double gsl_sf_erf_Q (double X)
 -- Function: int gsl_sf_erf_Q_e (double X, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_hazard (double X)
 -- Function: int gsl_sf_hazard_e (double X, gsl_sf_result * RESULT)
     These routines compute the hazard function for the normal
     distribution.

\1f
File: gsl-ref.info,  Node: Exponential Functions,  Next: Exponential Integrals,  Prev: Error Functions,  Up: Special Functions

7.16 Exponential Functions
==========================

The functions described in this section are declared in the header file
`gsl_sf_exp.h'.

* Menu:

* Exponential Function::
* Relative Exponential Functions::
* Exponentiation With Error Estimate::

\1f
File: gsl-ref.info,  Node: Exponential Function,  Next: Relative Exponential Functions,  Up: Exponential Functions

7.16.1 Exponential Function
---------------------------

 -- Function: double gsl_sf_exp (double X)
 -- Function: int gsl_sf_exp_e (double X, gsl_sf_result * RESULT)
     These routines provide an exponential function \exp(x) using GSL
     semantics and error checking.

 -- Function: int gsl_sf_exp_e10_e (double X, gsl_sf_result_e10 *
          RESULT)
     This function computes the exponential \exp(x) using the
     `gsl_sf_result_e10' type to return a result with extended range.
     This function may be useful if the value of \exp(x) would overflow
     the  numeric range of `double'.

 -- Function: double gsl_sf_exp_mult (double X, double Y)
 -- Function: int gsl_sf_exp_mult_e (double X, double Y, gsl_sf_result
          * RESULT)
     These routines exponentiate X and multiply by the factor Y to
     return the product y \exp(x).

 -- Function: int gsl_sf_exp_mult_e10_e (const double X, const double
          Y, gsl_sf_result_e10 * RESULT)
     This function computes the product y \exp(x) using the
     `gsl_sf_result_e10' type to return a result with extended numeric
     range.

\1f
File: gsl-ref.info,  Node: Relative Exponential Functions,  Next: Exponentiation With Error Estimate,  Prev: Exponential Function,  Up: Exponential Functions

7.16.2 Relative Exponential Functions
-------------------------------------

 -- Function: double gsl_sf_expm1 (double X)
 -- Function: int gsl_sf_expm1_e (double X, gsl_sf_result * RESULT)
     These routines compute the quantity \exp(x)-1 using an algorithm
     that is accurate for small x.

 -- Function: double gsl_sf_exprel (double X)
 -- Function: int gsl_sf_exprel_e (double X, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_exprel_2 (double X)
 -- Function: int gsl_sf_exprel_2_e (double X, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_exprel_n (int N, double X)
 -- Function: int gsl_sf_exprel_n_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the N-relative exponential, which is the
     N-th generalization of the functions `gsl_sf_exprel' and
     `gsl_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)

\1f
File: gsl-ref.info,  Node: Exponentiation With Error Estimate,  Prev: Relative Exponential Functions,  Up: Exponential Functions

7.16.3 Exponentiation With Error Estimate
-----------------------------------------

 -- Function: int gsl_sf_exp_err_e (double X, double DX, gsl_sf_result
          * RESULT)
     This function exponentiates X with an associated absolute error DX.

 -- Function: int gsl_sf_exp_err_e10_e (double X, double DX,
          gsl_sf_result_e10 * RESULT)
     This function exponentiates a quantity X with an associated
     absolute error DX using the `gsl_sf_result_e10' type to return a
     result with extended range.

 -- Function: int gsl_sf_exp_mult_err_e (double X, double DX, double Y,
          double DY, gsl_sf_result * RESULT)
     This routine computes the product y \exp(x) for the quantities X,
     Y with associated absolute errors DX, DY.

 -- Function: int gsl_sf_exp_mult_err_e10_e (double X, double DX,
          double Y, double DY, gsl_sf_result_e10 * RESULT)
     This routine computes the product y \exp(x) for the quantities X,
     Y with associated absolute errors DX, DY using the
     `gsl_sf_result_e10' type to return a result with extended range.

\1f
File: gsl-ref.info,  Node: Exponential Integrals,  Next: Fermi-Dirac Function,  Prev: Exponential Functions,  Up: Special Functions

7.17 Exponential Integrals
==========================

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'.

* Menu:

* Exponential Integral::
* Ei(x)::
* Hyperbolic Integrals::
* Ei_3(x)::
* Trigonometric Integrals::
* Arctangent Integral::

\1f
File: gsl-ref.info,  Node: Exponential Integral,  Next: Ei(x),  Up: Exponential Integrals

7.17.1 Exponential Integral
---------------------------

 -- Function: double gsl_sf_expint_E1 (double X)
 -- Function: int gsl_sf_expint_E1_e (double X, gsl_sf_result * RESULT)
     These routines compute the exponential integral E_1(x),

          E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t.



 -- Function: double gsl_sf_expint_E2 (double X)
 -- Function: int gsl_sf_expint_E2_e (double X, gsl_sf_result * RESULT)
     These routines compute the second-order exponential integral
     E_2(x),

          E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.



 -- Function: double gsl_sf_expint_En (int N, double X)
 -- Function: int gsl_sf_expint_En_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the exponential integral E_n(x) of order n,

          E_n(x) := \Re \int_1^\infty dt \exp(-xt)/t^n.



\1f
File: gsl-ref.info,  Node: Ei(x),  Next: Hyperbolic Integrals,  Prev: Exponential Integral,  Up: Exponential Integrals

7.17.2 Ei(x)
------------

 -- Function: double gsl_sf_expint_Ei (double X)
 -- Function: int gsl_sf_expint_Ei_e (double X, gsl_sf_result * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Hyperbolic Integrals,  Next: Ei_3(x),  Prev: Ei(x),  Up: Exponential Integrals

7.17.3 Hyperbolic Integrals
---------------------------

 -- Function: double gsl_sf_Shi (double X)
 -- Function: int gsl_sf_Shi_e (double X, gsl_sf_result * RESULT)
     These routines compute the integral Shi(x) = \int_0^x dt
     \sinh(t)/t.

 -- Function: double gsl_sf_Chi (double X)
 -- Function: int gsl_sf_Chi_e (double X, gsl_sf_result * RESULT)
     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').

\1f
File: gsl-ref.info,  Node: Ei_3(x),  Next: Trigonometric Integrals,  Prev: Hyperbolic Integrals,  Up: Exponential Integrals

7.17.4 Ei_3(x)
--------------

 -- Function: double gsl_sf_expint_3 (double X)
 -- Function: int gsl_sf_expint_3_e (double X, gsl_sf_result * RESULT)
     These routines compute the third-order exponential integral
     Ei_3(x) = \int_0^xdt \exp(-t^3) for x >= 0.

\1f
File: gsl-ref.info,  Node: Trigonometric Integrals,  Next: Arctangent Integral,  Prev: Ei_3(x),  Up: Exponential Integrals

7.17.5 Trigonometric Integrals
------------------------------

 -- Function: double gsl_sf_Si (const double X)
 -- Function: int gsl_sf_Si_e (double X, gsl_sf_result * RESULT)
     These routines compute the Sine integral Si(x) = \int_0^x dt
     \sin(t)/t.

 -- Function: double gsl_sf_Ci (const double X)
 -- Function: int gsl_sf_Ci_e (double X, gsl_sf_result * RESULT)
     These routines compute the Cosine integral Ci(x) = -\int_x^\infty
     dt \cos(t)/t for x > 0.

\1f
File: gsl-ref.info,  Node: Arctangent Integral,  Prev: Trigonometric Integrals,  Up: Exponential Integrals

7.17.6 Arctangent Integral
--------------------------

 -- Function: double gsl_sf_atanint (double X)
 -- Function: int gsl_sf_atanint_e (double X, gsl_sf_result * RESULT)
     These routines compute the Arctangent integral, which is defined
     as AtanInt(x) = \int_0^x dt \arctan(t)/t.

\1f
File: gsl-ref.info,  Node: Fermi-Dirac Function,  Next: Gamma and Beta Functions,  Prev: Exponential Integrals,  Up: Special Functions

7.18 Fermi-Dirac Function
=========================

The functions described in this section are declared in the header file
`gsl_sf_fermi_dirac.h'.

* Menu:

* Complete Fermi-Dirac Integrals::
* Incomplete Fermi-Dirac Integrals::

\1f
File: gsl-ref.info,  Node: Complete Fermi-Dirac Integrals,  Next: Incomplete Fermi-Dirac Integrals,  Up: Fermi-Dirac Function

7.18.1 Complete Fermi-Dirac Integrals
-------------------------------------

The complete Fermi-Dirac integral F_j(x) is given by,

     F_j(x)   := (1/\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1))
   Note that the Fermi-Dirac integral is sometimes defined without the
normalisation factor in other texts.

 -- Function: double gsl_sf_fermi_dirac_m1 (double X)
 -- Function: int gsl_sf_fermi_dirac_m1_e (double X, gsl_sf_result *
          RESULT)
     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).

 -- Function: double gsl_sf_fermi_dirac_0 (double X)
 -- Function: int gsl_sf_fermi_dirac_0_e (double X, gsl_sf_result *
          RESULT)
     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).

 -- Function: double gsl_sf_fermi_dirac_1 (double X)
 -- Function: int gsl_sf_fermi_dirac_1_e (double X, gsl_sf_result *
          RESULT)
     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)).

 -- Function: double gsl_sf_fermi_dirac_2 (double X)
 -- Function: int gsl_sf_fermi_dirac_2_e (double X, gsl_sf_result *
          RESULT)
     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)).

 -- Function: double gsl_sf_fermi_dirac_int (int J, double X)
 -- Function: int gsl_sf_fermi_dirac_int_e (int J, double X,
          gsl_sf_result * RESULT)
     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)).

 -- Function: double gsl_sf_fermi_dirac_mhalf (double X)
 -- Function: int gsl_sf_fermi_dirac_mhalf_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the complete Fermi-Dirac integral
     F_{-1/2}(x).

 -- Function: double gsl_sf_fermi_dirac_half (double X)
 -- Function: int gsl_sf_fermi_dirac_half_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the complete Fermi-Dirac integral
     F_{1/2}(x).

 -- Function: double gsl_sf_fermi_dirac_3half (double X)
 -- Function: int gsl_sf_fermi_dirac_3half_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the complete Fermi-Dirac integral
     F_{3/2}(x).

\1f
File: gsl-ref.info,  Node: Incomplete Fermi-Dirac Integrals,  Prev: Complete Fermi-Dirac Integrals,  Up: Fermi-Dirac Function

7.18.2 Incomplete Fermi-Dirac Integrals
---------------------------------------

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))

 -- Function: double gsl_sf_fermi_dirac_inc_0 (double X, double B)
 -- Function: int gsl_sf_fermi_dirac_inc_0_e (double X, double B,
          gsl_sf_result * RESULT)
     These routines compute the incomplete Fermi-Dirac integral with an
     index of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x).

\1f
File: gsl-ref.info,  Node: Gamma and Beta Functions,  Next: Gegenbauer Functions,  Prev: Fermi-Dirac Function,  Up: Special Functions

7.19 Gamma and Beta Functions
=============================

The  functions described in this section are declared in the header
file `gsl_sf_gamma.h'.

* Menu:

* Gamma Functions::
* Factorials::
* Pochhammer Symbol::
* Incomplete Gamma Functions::
* Beta Functions::
* Incomplete Beta Function::

\1f
File: gsl-ref.info,  Node: Gamma Functions,  Next: Factorials,  Up: Gamma and Beta Functions

7.19.1 Gamma Functions
----------------------

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'.

 -- Function: double gsl_sf_gamma (double X)
 -- Function: int gsl_sf_gamma_e (double X, gsl_sf_result * RESULT)
     These routines compute the Gamma function \Gamma(x), subject to x
     not being a negative integer or zero.  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_XMAX' and is 171.0.

 -- Function: double gsl_sf_lngamma (double X)
 -- Function: int gsl_sf_lngamma_e (double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of the Gamma function,
     \log(\Gamma(x)), subject to x not being a negative integer or
     zero.  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.

 -- Function: int gsl_sf_lngamma_sgn_e (double X, gsl_sf_result *
          RESULT_LG, double * SGN)
     This routine computes the sign of the gamma function and the
     logarithm of its magnitude, subject to x not being a negative
     integer or zero.  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).

 -- Function: double gsl_sf_gammastar (double X)
 -- Function: int gsl_sf_gammastar_e (double X, gsl_sf_result * RESULT)
     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 \infty
     and is a useful suggestion of Temme.

 -- Function: double gsl_sf_gammainv (double X)
 -- Function: int gsl_sf_gammainv_e (double X, gsl_sf_result * RESULT)
     These routines compute the reciprocal of the gamma function,
     1/\Gamma(x) using the real Lanczos method.

 -- Function: int gsl_sf_lngamma_complex_e (double ZR, double ZI,
          gsl_sf_result * LNR, gsl_sf_result * ARG)
     This routine computes \log(\Gamma(z)) for complex z=z_r+i z_i and
     z not a negative integer or zero, 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_ELOSS' error when it occurs.  The absolute value part (LNR),
     however, never suffers from loss of precision.

\1f
File: gsl-ref.info,  Node: Factorials,  Next: Pochhammer Symbol,  Prev: Gamma Functions,  Up: Gamma and Beta Functions

7.19.2 Factorials
-----------------

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.

 -- Function: double gsl_sf_fact (unsigned int N)
 -- Function: int gsl_sf_fact_e (unsigned int N, gsl_sf_result * RESULT)
     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_NMAX' and is 170.

 -- Function: double gsl_sf_doublefact (unsigned int N)
 -- Function: int gsl_sf_doublefact_e (unsigned int N, gsl_sf_result *
          RESULT)
     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_NMAX' and is 297.

 -- Function: double gsl_sf_lnfact (unsigned int N)
 -- Function: int gsl_sf_lnfact_e (unsigned int N, gsl_sf_result *
          RESULT)
     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_lngamma' for n < 170, but defers for larger N.

 -- Function: double gsl_sf_lndoublefact (unsigned int N)
 -- Function: int gsl_sf_lndoublefact_e (unsigned int N, gsl_sf_result
          * RESULT)
     These routines compute the logarithm of the double factorial of N,
     \log(n!!).

 -- Function: double gsl_sf_choose (unsigned int N, unsigned int M)
 -- Function: int gsl_sf_choose_e (unsigned int N, unsigned int M,
          gsl_sf_result * RESULT)
     These routines compute the combinatorial factor `n choose m' =
     n!/(m!(n-m)!)

 -- Function: double gsl_sf_lnchoose (unsigned int N, unsigned int M)
 -- Function: int gsl_sf_lnchoose_e (unsigned int N, unsigned int M,
          gsl_sf_result * RESULT)
     These routines compute the logarithm of `n choose m'.  This is
     equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).

 -- Function: double gsl_sf_taylorcoeff (int N, double X)
 -- Function: int gsl_sf_taylorcoeff_e (int N, double X, gsl_sf_result
          * RESULT)
     These routines compute the Taylor coefficient x^n / n! for x >= 0,
     n >= 0.

\1f
File: gsl-ref.info,  Node: Pochhammer Symbol,  Next: Incomplete Gamma Functions,  Prev: Factorials,  Up: Gamma and Beta Functions

7.19.3 Pochhammer Symbol
------------------------

 -- Function: double gsl_sf_poch (double A, double X)
 -- Function: int gsl_sf_poch_e (double A, double X, gsl_sf_result *
          RESULT)
     These routines compute the Pochhammer symbol (a)_x = \Gamma(a +
     x)/\Gamma(a), subject to a and a+x not being negative integers or
     zero. The Pochhammer symbol is also known as the Apell symbol and
     sometimes written as (a,x).

 -- Function: double gsl_sf_lnpoch (double A, double X)
 -- Function: int gsl_sf_lnpoch_e (double A, double X, gsl_sf_result *
          RESULT)
     These routines compute the logarithm of the Pochhammer symbol,
     \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)) for a > 0, a+x > 0.

 -- Function: int gsl_sf_lnpoch_sgn_e (double A, double X,
          gsl_sf_result * RESULT, double * SGN)
     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 or
     zero.

 -- Function: double gsl_sf_pochrel (double A, double X)
 -- Function: int gsl_sf_pochrel_e (double A, double X, gsl_sf_result *
          RESULT)
     These routines compute the relative Pochhammer symbol ((a)_x -
     1)/x where (a)_x = \Gamma(a + x)/\Gamma(a).

\1f
File: gsl-ref.info,  Node: Incomplete Gamma Functions,  Next: Beta Functions,  Prev: Pochhammer Symbol,  Up: Gamma and Beta Functions

7.19.4 Incomplete Gamma Functions
---------------------------------

 -- Function: double gsl_sf_gamma_inc (double A, double X)
 -- Function: int gsl_sf_gamma_inc_e (double A, double X, gsl_sf_result
          * RESULT)
     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.

 -- Function: double gsl_sf_gamma_inc_Q (double A, double X)
 -- Function: int gsl_sf_gamma_inc_Q_e (double A, double X,
          gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_gamma_inc_P (double A, double X)
 -- Function: int gsl_sf_gamma_inc_P_e (double A, double X,
          gsl_sf_result * RESULT)
     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).

\1f
File: gsl-ref.info,  Node: Beta Functions,  Next: Incomplete Beta Function,  Prev: Incomplete Gamma Functions,  Up: Gamma and Beta Functions

7.19.5 Beta Functions
---------------------

 -- Function: double gsl_sf_beta (double A, double B)
 -- Function: int gsl_sf_beta_e (double A, double B, gsl_sf_result *
          RESULT)
     These routines compute the Beta Function, B(a,b) =
     \Gamma(a)\Gamma(b)/\Gamma(a+b) subject to a and b not being
     negative integers.

 -- Function: double gsl_sf_lnbeta (double A, double B)
 -- Function: int gsl_sf_lnbeta_e (double A, double B, gsl_sf_result *
          RESULT)
     These routines compute the logarithm of the Beta Function,
     \log(B(a,b)) subject to a and b not being negative integers.

\1f
File: gsl-ref.info,  Node: Incomplete Beta Function,  Prev: Beta Functions,  Up: Gamma and Beta Functions

7.19.6 Incomplete Beta Function
-------------------------------

 -- Function: double gsl_sf_beta_inc (double A, double B, double X)
 -- Function: int gsl_sf_beta_inc_e (double A, double B, double X,
          gsl_sf_result * RESULT)
     These routines compute the normalized incomplete Beta function
     I_x(a,b)=B_x(a,b)/B(a,b) where B_x(a,b) = \int_0^x t^{a-1}
     (1-t)^{b-1} dt for 0 <= x <= 1.   For a > 0, b > 0 the value is
     computed using a continued fraction expansion.  For all other
     values it is computed using the relation I_x(a,b,x) = (1/a) x^a
     2F1(a,1-b,a+1,x)/B(a,b).

\1f
File: gsl-ref.info,  Node: Gegenbauer Functions,  Next: Hypergeometric Functions,  Prev: Gamma and Beta Functions,  Up: Special Functions

7.20 Gegenbauer Functions
=========================

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'.

 -- Function: double gsl_sf_gegenpoly_1 (double LAMBDA, double X)
 -- Function: double gsl_sf_gegenpoly_2 (double LAMBDA, double X)
 -- Function: double gsl_sf_gegenpoly_3 (double LAMBDA, double X)
 -- Function: int gsl_sf_gegenpoly_1_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
 -- Function: int gsl_sf_gegenpoly_2_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
 -- Function: int gsl_sf_gegenpoly_3_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
     These functions evaluate the Gegenbauer polynomials
     C^{(\lambda)}_n(x) using explicit representations for n =1, 2, 3.

 -- Function: double gsl_sf_gegenpoly_n (int N, double LAMBDA, double X)
 -- Function: int gsl_sf_gegenpoly_n_e (int N, double LAMBDA, double X,
          gsl_sf_result * RESULT)
     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.

 -- Function: int gsl_sf_gegenpoly_array (int NMAX, double LAMBDA,
          double X, double RESULT_ARRAY[])
     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.

\1f
File: gsl-ref.info,  Node: Hypergeometric Functions,  Next: Laguerre Functions,  Prev: Gegenbauer Functions,  Up: Special Functions

7.21 Hypergeometric Functions
=============================

Hypergeometric functions are described in Abramowitz & Stegun, Chapters
13 and 15.  These functions are declared in the header file
`gsl_sf_hyperg.h'.

 -- Function: double gsl_sf_hyperg_0F1 (double C, double X)
 -- Function: int gsl_sf_hyperg_0F1_e (double C, double X,
          gsl_sf_result * RESULT)
     These routines compute the hypergeometric function 0F1(c,x).

 -- Function: double gsl_sf_hyperg_1F1_int (int M, int N, double X)
 -- Function: int gsl_sf_hyperg_1F1_int_e (int M, int N, double X,
          gsl_sf_result * RESULT)
     These routines compute the confluent hypergeometric function
     1F1(m,n,x) = M(m,n,x) for integer parameters M, N.

 -- Function: double gsl_sf_hyperg_1F1 (double A, double B, double X)
 -- Function: int gsl_sf_hyperg_1F1_e (double A, double B, double X,
          gsl_sf_result * RESULT)
     These routines compute the confluent hypergeometric function
     1F1(a,b,x) = M(a,b,x) for general parameters A, B.

 -- Function: double gsl_sf_hyperg_U_int (int M, int N, double X)
 -- Function: int gsl_sf_hyperg_U_int_e (int M, int N, double X,
          gsl_sf_result * RESULT)
     These routines compute the confluent hypergeometric function
     U(m,n,x) for integer parameters M, N.

 -- Function: int gsl_sf_hyperg_U_int_e10_e (int M, int N, double X,
          gsl_sf_result_e10 * RESULT)
     This routine computes the confluent hypergeometric function
     U(m,n,x) for integer parameters M, N using the `gsl_sf_result_e10'
     type to return a result with extended range.

 -- Function: double gsl_sf_hyperg_U (double A, double B, double X)
 -- Function: int gsl_sf_hyperg_U_e (double A, double B, double X,
          gsl_sf_result * RESULT)
     These routines compute the confluent hypergeometric function
     U(a,b,x).

 -- Function: int gsl_sf_hyperg_U_e10_e (double A, double B, double X,
          gsl_sf_result_e10 * RESULT)
     This routine computes the confluent hypergeometric function
     U(a,b,x) using the `gsl_sf_result_e10' type to return a result
     with extended range.

 -- Function: double gsl_sf_hyperg_2F1 (double A, double B, double C,
          double X)
 -- Function: int gsl_sf_hyperg_2F1_e (double A, double B, double C,
          double X, gsl_sf_result * RESULT)
     These routines compute the Gauss hypergeometric function
     2F1(a,b,c,x) = F(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_EMAXITER' when the series
     approximation converges too slowly.  This occurs in the region of
     x=1, c - a - b = m for integer m.

 -- Function: double gsl_sf_hyperg_2F1_conj (double AR, double AI,
          double C, double X)
 -- Function: int gsl_sf_hyperg_2F1_conj_e (double AR, double AI,
          double C, double X, gsl_sf_result * RESULT)
     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:

 -- Function: double gsl_sf_hyperg_2F1_renorm (double A, double B,
          double C, double X)
 -- Function: int gsl_sf_hyperg_2F1_renorm_e (double A, double B,
          double C, double X, gsl_sf_result * RESULT)
     These routines compute the renormalized Gauss hypergeometric
     function 2F1(a,b,c,x) / \Gamma(c) for |x| < 1.

 -- Function: double gsl_sf_hyperg_2F1_conj_renorm (double AR, double
          AI, double C, double X)
 -- Function: int gsl_sf_hyperg_2F1_conj_renorm_e (double AR, double
          AI, double C, double X, gsl_sf_result * RESULT)
     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.

 -- Function: double gsl_sf_hyperg_2F0 (double A, double B, double X)
 -- Function: int gsl_sf_hyperg_2F0_e (double A, double B, double X,
          gsl_sf_result * RESULT)
     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)

\1f
File: gsl-ref.info,  Node: Laguerre Functions,  Next: Lambert W Functions,  Prev: Hypergeometric Functions,  Up: Special Functions

7.22 Laguerre Functions
=======================

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'.

 -- Function: double gsl_sf_laguerre_1 (double A, double X)
 -- Function: double gsl_sf_laguerre_2 (double A, double X)
 -- Function: double gsl_sf_laguerre_3 (double A, double X)
 -- Function: int gsl_sf_laguerre_1_e (double A, double X,
          gsl_sf_result * RESULT)
 -- Function: int gsl_sf_laguerre_2_e (double A, double X,
          gsl_sf_result * RESULT)
 -- Function: int gsl_sf_laguerre_3_e (double A, double X,
          gsl_sf_result * RESULT)
     These routines evaluate the generalized Laguerre polynomials
     L^a_1(x), L^a_2(x), L^a_3(x) using explicit representations.

 -- Function: double gsl_sf_laguerre_n (const int N, const double A,
          const double X)
 -- Function: int gsl_sf_laguerre_n_e (int N, double A, double X,
          gsl_sf_result * RESULT)
     These routines evaluate the generalized Laguerre polynomials
     L^a_n(x) for a > -1, n >= 0.


\1f
File: gsl-ref.info,  Node: Lambert W Functions,  Next: Legendre Functions and Spherical Harmonics,  Prev: Laguerre Functions,  Up: Special Functions

7.23 Lambert W Functions
========================

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'.

 -- Function: double gsl_sf_lambert_W0 (double X)
 -- Function: int gsl_sf_lambert_W0_e (double X, gsl_sf_result * RESULT)
     These compute the principal branch of the Lambert W function,
     W_0(x).

 -- Function: double gsl_sf_lambert_Wm1 (double X)
 -- Function: int gsl_sf_lambert_Wm1_e (double X, gsl_sf_result *
          RESULT)
     These compute the secondary real-valued branch of the Lambert W
     function, W_{-1}(x).

\1f
File: gsl-ref.info,  Node: Legendre Functions and Spherical Harmonics,  Next: Logarithm and Related Functions,  Prev: Lambert W Functions,  Up: Special Functions

7.24 Legendre Functions and Spherical Harmonics
===============================================

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'.

* Menu:

* Legendre Polynomials::
* Associated Legendre Polynomials and Spherical Harmonics::
* Conical Functions::
* Radial Functions for Hyperbolic Space::

\1f
File: gsl-ref.info,  Node: Legendre Polynomials,  Next: Associated Legendre Polynomials and Spherical Harmonics,  Up: Legendre Functions and Spherical Harmonics

7.24.1 Legendre Polynomials
---------------------------

 -- Function: double gsl_sf_legendre_P1 (double X)
 -- Function: double gsl_sf_legendre_P2 (double X)
 -- Function: double gsl_sf_legendre_P3 (double X)
 -- Function: int gsl_sf_legendre_P1_e (double X, gsl_sf_result *
          RESULT)
 -- Function: int gsl_sf_legendre_P2_e (double X, gsl_sf_result *
          RESULT)
 -- Function: int gsl_sf_legendre_P3_e (double X, gsl_sf_result *
          RESULT)
     These functions evaluate the Legendre polynomials P_l(x) using
     explicit representations for l=1, 2, 3.

 -- Function: double gsl_sf_legendre_Pl (int L, double X)
 -- Function: int gsl_sf_legendre_Pl_e (int L, double X, gsl_sf_result
          * RESULT)
     These functions evaluate the Legendre polynomial P_l(x) for a
     specific value of L, X subject to l >= 0, |x| <= 1

 -- Function: int gsl_sf_legendre_Pl_array (int LMAX, double X, double
          RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_Pl_deriv_array (int LMAX, double X,
          double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     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

 -- Function: double gsl_sf_legendre_Q0 (double X)
 -- Function: int gsl_sf_legendre_Q0_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the Legendre function Q_0(x) for x > -1, x
     != 1.

 -- Function: double gsl_sf_legendre_Q1 (double X)
 -- Function: int gsl_sf_legendre_Q1_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the Legendre function Q_1(x) for x > -1, x
     != 1.

 -- Function: double gsl_sf_legendre_Ql (int L, double X)
 -- Function: int gsl_sf_legendre_Ql_e (int L, double X, gsl_sf_result
          * RESULT)
     These routines compute the Legendre function Q_l(x) for x > -1, x
     != 1 and l >= 0.

\1f
File: gsl-ref.info,  Node: Associated Legendre Polynomials and Spherical Harmonics,  Next: Conical Functions,  Prev: Legendre Polynomials,  Up: Legendre Functions and Spherical Harmonics

7.24.2 Associated Legendre Polynomials and Spherical Harmonics
--------------------------------------------------------------

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.

 -- Function: double gsl_sf_legendre_Plm (int L, int M, double X)
 -- Function: int gsl_sf_legendre_Plm_e (int L, int M, double X,
          gsl_sf_result * RESULT)
     These routines compute the associated Legendre polynomial P_l^m(x)
     for m >= 0, l >= m, |x| <= 1.

 -- Function: int gsl_sf_legendre_Plm_array (int LMAX, int M, double X,
          double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_Plm_deriv_array (int LMAX, int M,
          double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     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.

 -- Function: double gsl_sf_legendre_sphPlm (int L, int M, double X)
 -- Function: int gsl_sf_legendre_sphPlm_e (int L, int M, double X,
          gsl_sf_result * RESULT)
     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).

 -- Function: int gsl_sf_legendre_sphPlm_array (int LMAX, int M, double
          X, double RESULT_ARRAY[])
 -- Function: int gsl_sf_legendre_sphPlm_deriv_array (int LMAX, int M,
          double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
     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

 -- Function: int gsl_sf_legendre_array_size (const int LMAX, const int
          M)
     This function returns the size of RESULT_ARRAY[] needed for the
     array versions of P_l^m(x), LMAX - M + 1.  An inline version of
     this function is used when `HAVE_INLINE' is defined.

\1f
File: gsl-ref.info,  Node: Conical Functions,  Next: Radial Functions for Hyperbolic Space,  Prev: Associated Legendre Polynomials and Spherical Harmonics,  Up: Legendre Functions and Spherical Harmonics

7.24.3 Conical Functions
------------------------

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.

 -- Function: double gsl_sf_conicalP_half (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_half_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
     These routines compute the irregular Spherical Conical Function
     P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_mhalf (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_mhalf_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
     These routines compute the regular Spherical Conical Function
     P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_0 (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_0_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
     These routines compute the conical function P^0_{-1/2 + i
     \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_1 (double LAMBDA, double X)
 -- Function: int gsl_sf_conicalP_1_e (double LAMBDA, double X,
          gsl_sf_result * RESULT)
     These routines compute the conical function P^1_{-1/2 + i
     \lambda}(x) for x > -1.

 -- Function: double gsl_sf_conicalP_sph_reg (int L, double LAMBDA,
          double X)
 -- Function: int gsl_sf_conicalP_sph_reg_e (int L, double LAMBDA,
          double X, gsl_sf_result * RESULT)
     These routines compute the Regular Spherical Conical Function
     P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.

 -- Function: double gsl_sf_conicalP_cyl_reg (int M, double LAMBDA,
          double X)
 -- Function: int gsl_sf_conicalP_cyl_reg_e (int M, double LAMBDA,
          double X, gsl_sf_result * RESULT)
     These routines compute the Regular Cylindrical Conical Function
     P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1.

\1f
File: gsl-ref.info,  Node: Radial Functions for Hyperbolic Space,  Prev: Conical Functions,  Up: Legendre Functions and Spherical Harmonics

7.24.4 Radial Functions for Hyperbolic Space
--------------------------------------------

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.

 -- Function: double gsl_sf_legendre_H3d_0 (double LAMBDA, double ETA)
 -- Function: int gsl_sf_legendre_H3d_0_e (double LAMBDA, double ETA,
          gsl_sf_result * RESULT)
     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).

 -- Function: double gsl_sf_legendre_H3d_1 (double LAMBDA, double ETA)
 -- Function: int gsl_sf_legendre_H3d_1_e (double LAMBDA, double ETA,
          gsl_sf_result * RESULT)
     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).

 -- Function: double gsl_sf_legendre_H3d (int L, double LAMBDA, double
          ETA)
 -- Function: int gsl_sf_legendre_H3d_e (int L, double LAMBDA, double
          ETA, gsl_sf_result * RESULT)
     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).

 -- Function: int gsl_sf_legendre_H3d_array (int LMAX, double LAMBDA,
          double ETA, double RESULT_ARRAY[])
     This function computes an array of radial eigenfunctions
     L^{H3d}_l(\lambda, \eta) for 0 <= l <= lmax.

\1f
File: gsl-ref.info,  Node: Logarithm and Related Functions,  Next: Mathieu Functions,  Prev: Legendre Functions and Spherical Harmonics,  Up: Special Functions

7.25 Logarithm and Related Functions
====================================

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'.

 -- Function: double gsl_sf_log (double X)
 -- Function: int gsl_sf_log_e (double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of X, \log(x), for x > 0.

 -- Function: double gsl_sf_log_abs (double X)
 -- Function: int gsl_sf_log_abs_e (double X, gsl_sf_result * RESULT)
     These routines compute the logarithm of the magnitude of X,
     \log(|x|), for x \ne 0.

 -- Function: int gsl_sf_complex_log_e (double ZR, double ZI,
          gsl_sf_result * LNR, gsl_sf_result * THETA)
     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].

 -- Function: double gsl_sf_log_1plusx (double X)
 -- Function: int gsl_sf_log_1plusx_e (double X, gsl_sf_result * RESULT)
     These routines compute \log(1 + x) for x > -1 using an algorithm
     that is accurate for small x.

 -- Function: double gsl_sf_log_1plusx_mx (double X)
 -- Function: int gsl_sf_log_1plusx_mx_e (double X, gsl_sf_result *
          RESULT)
     These routines compute \log(1 + x) - x for x > -1 using an
     algorithm that is accurate for small x.

\1f
File: gsl-ref.info,  Node: Mathieu Functions,  Next: Power Function,  Prev: Logarithm and Related Functions,  Up: Special Functions

7.26 Mathieu Functions
======================

The routines described in this section compute the angular and radial
Mathieu functions, and their characteristic values.  Mathieu functions
are the solutions of the following two differential equations:

     d^2y/dv^2 + (a - 2q\cos 2v)y = 0
     d^2f/du^2 - (a - 2q\cosh 2u)f = 0

The angular Mathieu functions ce_r(x,q), se_r(x,q) are the even and odd
periodic solutions of the first equation, which is known as Mathieu's
equation. These exist only for the discrete sequence of  characteristic
values a=a_r(q) (even-periodic) and a=b_r(q) (odd-periodic).

   The radial Mathieu functions Mc^{(j)}_{r}(z,q), Ms^{(j)}_{r}(z,q)
are the solutions of the second equation, which is referred to as
Mathieu's modified equation.  The radial Mathieu functions of the
first, second, third and fourth kind are denoted by the parameter j,
which takes the value 1, 2, 3 or 4.

   For more information on the Mathieu functions, see Abramowitz and
Stegun, Chapter 20.  These functions are defined in the header file
`gsl_sf_mathieu.h'.

* Menu:

* Mathieu Function Workspace::
* Mathieu Function Characteristic Values::
* Angular Mathieu Functions::
* Radial Mathieu Functions::

\1f
File: gsl-ref.info,  Node: Mathieu Function Workspace,  Next: Mathieu Function Characteristic Values,  Up: Mathieu Functions

7.26.1 Mathieu Function Workspace
---------------------------------

The Mathieu functions can be computed for a single order or for
multiple orders, using array-based routines.  The array-based routines
require a preallocated workspace.

 -- Function: gsl_sf_mathieu_workspace * gsl_sf_mathieu_alloc (size_t
          N, double QMAX)
     This function returns a workspace for the array versions of the
     Mathieu routines.  The arguments N and QMAX specify the maximum
     order and q-value of Mathieu functions which can be computed with
     this workspace.


 -- Function: void gsl_sf_mathieu_free (gsl_sf_mathieu_workspace *WORK)
     This function frees the workspace WORK.

\1f
File: gsl-ref.info,  Node: Mathieu Function Characteristic Values,  Next: Angular Mathieu Functions,  Prev: Mathieu Function Workspace,  Up: Mathieu Functions

7.26.2 Mathieu Function Characteristic Values
---------------------------------------------

 -- Function: int gsl_sf_mathieu_a (int N, double Q, gsl_sf_result
          *RESULT)
 -- Function: int gsl_sf_mathieu_b (int N, double Q, gsl_sf_result
          *RESULT)
     These routines compute the characteristic values a_n(q), b_n(q) of
     the Mathieu functions ce_n(q,x) and se_n(q,x), respectively.

 -- Function: int gsl_sf_mathieu_a_array (int ORDER_MIN, int ORDER_MAX,
          double Q, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_b_array (int ORDER_MIN, int ORDER_MAX,
          double Q, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
     These routines compute a series of Mathieu characteristic values
     a_n(q), b_n(q) for n from ORDER_MIN to ORDER_MAX inclusive,
     storing the results in the array RESULT_ARRAY.

\1f
File: gsl-ref.info,  Node: Angular Mathieu Functions,  Next: Radial Mathieu Functions,  Prev: Mathieu Function Characteristic Values,  Up: Mathieu Functions

7.26.3 Angular Mathieu Functions
--------------------------------

 -- Function: int gsl_sf_mathieu_ce (int N, double Q, double X,
          gsl_sf_result *RESULT)
 -- Function: int gsl_sf_mathieu_se (int N, double Q, double X,
          gsl_sf_result *RESULT)
     These routines compute the angular Mathieu functions ce_n(q,x) and
     se_n(q,x), respectively.

 -- Function: int gsl_sf_mathieu_ce_array (int NMIN, int NMAX, double
          Q, double X, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_se_array (int NMIN, int NMAX, double
          Q, double X, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
     These routines compute a series of the angular Mathieu functions
     ce_n(q,x) and se_n(q,x) of order n from NMIN to NMAX inclusive,
     storing the results in the array RESULT_ARRAY.

\1f
File: gsl-ref.info,  Node: Radial Mathieu Functions,  Prev: Angular Mathieu Functions,  Up: Mathieu Functions

7.26.4 Radial Mathieu Functions
-------------------------------

 -- Function: int gsl_sf_mathieu_Mc (int J, int N, double Q, double X,
          gsl_sf_result *RESULT)
 -- Function: int gsl_sf_mathieu_Ms (int J, int N, double Q, double X,
          gsl_sf_result *RESULT)
     These routines compute the radial J-th kind Mathieu functions
     Mc_n^{(j)}(q,x) and Ms_n^{(j)}(q,x) of order N.

     The allowed values of J are 1 and 2.  The functions for j = 3,4
     can be computed as M_n^{(3)} = M_n^{(1)} + iM_n^{(2)} and
     M_n^{(4)} = M_n^{(1)} - iM_n^{(2)}, where M_n^{(j)} = Mc_n^{(j)} or
     Ms_n^{(j)}.

 -- Function: int gsl_sf_mathieu_Mc_array (int J, int NMIN, int NMAX,
          double Q, double X, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
 -- Function: int gsl_sf_mathieu_Ms_array (int J, int NMIN, int NMAX,
          double Q, double X, gsl_sf_mathieu_workspace *WORK, double
          RESULT_ARRAY[])
     These routines compute a series of the radial Mathieu functions of
     kind J, with order from NMIN to NMAX inclusive, storing the
     results in the array RESULT_ARRAY.

\1f
File: gsl-ref.info,  Node: Power Function,  Next: Psi (Digamma) Function,  Prev: Mathieu Functions,  Up: Special Functions

7.27 Power Function
===================

The following functions are equivalent to the function `gsl_pow_int'
(*note Small integer powers::) with an error estimate.  These functions
are declared in the header file `gsl_sf_pow_int.h'.

 -- Function: double gsl_sf_pow_int (double X, int N)
 -- Function: int gsl_sf_pow_int_e (double X, int N, gsl_sf_result *
          RESULT)
     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);

\1f
File: gsl-ref.info,  Node: Psi (Digamma) Function,  Next: Synchrotron Functions,  Prev: Power Function,  Up: Special Functions

7.28 Psi (Digamma) Function
===========================

The polygamma functions of order n are defined by

     \psi^{(n)}(x) = (d/dx)^n \psi(x) = (d/dx)^{n+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'.

* Menu:

* Digamma Function::
* Trigamma Function::
* Polygamma Function::

\1f
File: gsl-ref.info,  Node: Digamma Function,  Next: Trigamma Function,  Up: Psi (Digamma) Function

7.28.1 Digamma Function
-----------------------

 -- Function: double gsl_sf_psi_int (int N)
 -- Function: int gsl_sf_psi_int_e (int N, gsl_sf_result * RESULT)
     These routines compute the digamma function \psi(n) for positive
     integer N.  The digamma function is also called the Psi function.

 -- Function: double gsl_sf_psi (double X)
 -- Function: int gsl_sf_psi_e (double X, gsl_sf_result * RESULT)
     These routines compute the digamma function \psi(x) for general x,
     x \ne 0.

 -- Function: double gsl_sf_psi_1piy (double Y)
 -- Function: int gsl_sf_psi_1piy_e (double Y, gsl_sf_result * RESULT)
     These routines compute the real part of the digamma function on
     the line 1+i y, \Re[\psi(1 + i y)].

\1f
File: gsl-ref.info,  Node: Trigamma Function,  Next: Polygamma Function,  Prev: Digamma Function,  Up: Psi (Digamma) Function

7.28.2 Trigamma Function
------------------------

 -- Function: double gsl_sf_psi_1_int (int N)
 -- Function: int gsl_sf_psi_1_int_e (int N, gsl_sf_result * RESULT)
     These routines compute the Trigamma function \psi'(n) for positive
     integer n.

 -- Function: double gsl_sf_psi_1 (double X)
 -- Function: int gsl_sf_psi_1_e (double X, gsl_sf_result * RESULT)
     These routines compute the Trigamma function \psi'(x) for general
     x.

\1f
File: gsl-ref.info,  Node: Polygamma Function,  Prev: Trigamma Function,  Up: Psi (Digamma) Function

7.28.3 Polygamma Function
-------------------------

 -- Function: double gsl_sf_psi_n (int N, double X)
 -- Function: int gsl_sf_psi_n_e (int N, double X, gsl_sf_result *
          RESULT)
     These routines compute the polygamma function \psi^{(n)}(x) for n
     >= 0, x > 0.

\1f
File: gsl-ref.info,  Node: Synchrotron Functions,  Next: Transport Functions,  Prev: Psi (Digamma) Function,  Up: Special Functions

7.29 Synchrotron Functions
==========================

The functions described in this section are declared in the header file
`gsl_sf_synchrotron.h'.

 -- Function: double gsl_sf_synchrotron_1 (double X)
 -- Function: int gsl_sf_synchrotron_1_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the first synchrotron function x
     \int_x^\infty dt K_{5/3}(t) for x >= 0.

 -- Function: double gsl_sf_synchrotron_2 (double X)
 -- Function: int gsl_sf_synchrotron_2_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the second synchrotron function x
     K_{2/3}(x) for x >= 0.

\1f
File: gsl-ref.info,  Node: Transport Functions,  Next: Trigonometric Functions,  Prev: Synchrotron Functions,  Up: Special Functions

7.30 Transport Functions
========================

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'.

 -- Function: double gsl_sf_transport_2 (double X)
 -- Function: int gsl_sf_transport_2_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the transport function J(2,x).

 -- Function: double gsl_sf_transport_3 (double X)
 -- Function: int gsl_sf_transport_3_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the transport function J(3,x).

 -- Function: double gsl_sf_transport_4 (double X)
 -- Function: int gsl_sf_transport_4_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the transport function J(4,x).

 -- Function: double gsl_sf_transport_5 (double X)
 -- Function: int gsl_sf_transport_5_e (double X, gsl_sf_result *
          RESULT)
     These routines compute the transport function J(5,x).

\1f
File: gsl-ref.info,  Node: Trigonometric Functions,  Next: Zeta Functions,  Prev: Transport Functions,  Up: Special Functions

7.31 Trigonometric Functions
============================

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'.

* Menu:

* Circular Trigonometric Functions::
* Trigonometric Functions for Complex Arguments::
* Hyperbolic Trigonometric Functions::
* Conversion Functions::
* Restriction Functions::
* Trigonometric Functions With Error Estimates::

\1f
File: gsl-ref.info,  Node: Circular Trigonometric Functions,  Next: Trigonometric Functions for Complex Arguments,  Up: Trigonometric Functions

7.31.1 Circular Trigonometric Functions
---------------------------------------

 -- Function: double gsl_sf_sin (double X)
 -- Function: int gsl_sf_sin_e (double X, gsl_sf_result * RESULT)
     These routines compute the sine function \sin(x).

 -- Function: double gsl_sf_cos (double X)
 -- Function: int gsl_sf_cos_e (double X, gsl_sf_result * RESULT)
     These routines compute the cosine function \cos(x).

 -- Function: double gsl_sf_hypot (double X, double Y)
 -- Function: int gsl_sf_hypot_e (double X, double Y, gsl_sf_result *
          RESULT)
     These routines compute the hypotenuse function \sqrt{x^2 + y^2}
     avoiding overflow and underflow.

 -- Function: double gsl_sf_sinc (double X)
 -- Function: int gsl_sf_sinc_e (double X, gsl_sf_result * RESULT)
     These routines compute \sinc(x) = \sin(\pi x) / (\pi x) for any
     value of X.

\1f
File: gsl-ref.info,  Node: Trigonometric Functions for Complex Arguments,  Next: Hyperbolic Trigonometric Functions,  Prev: Circular Trigonometric Functions,  Up: Trigonometric Functions

7.31.2 Trigonometric Functions for Complex Arguments
----------------------------------------------------

 -- Function: int gsl_sf_complex_sin_e (double ZR, double ZI,
          gsl_sf_result * SZR, gsl_sf_result * SZI)
     This function computes the complex sine, \sin(z_r + i z_i) storing
     the real and imaginary parts in SZR, SZI.

 -- Function: int gsl_sf_complex_cos_e (double ZR, double ZI,
          gsl_sf_result * CZR, gsl_sf_result * CZI)
     This function computes the complex cosine, \cos(z_r + i z_i)
     storing the real and imaginary parts in SZR, SZI.

 -- Function: int gsl_sf_complex_logsin_e (double ZR, double ZI,
          gsl_sf_result * LSZR, gsl_sf_result * LSZI)
     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.

\1f
File: gsl-ref.info,  Node: Hyperbolic Trigonometric Functions,  Next: Conversion Functions,  Prev: Trigonometric Functions for Complex Arguments,  Up: Trigonometric Functions

7.31.3 Hyperbolic Trigonometric Functions
-----------------------------------------

 -- Function: double gsl_sf_lnsinh (double X)
 -- Function: int gsl_sf_lnsinh_e (double X, gsl_sf_result * RESULT)
     These routines compute \log(\sinh(x)) for x > 0.

 -- Function: double gsl_sf_lncosh (double X)
 -- Function: int gsl_sf_lncosh_e (double X, gsl_sf_result * RESULT)
     These routines compute \log(\cosh(x)) for any X.

\1f
File: gsl-ref.info,  Node: Conversion Functions,  Next: Restriction Functions,  Prev: Hyperbolic Trigonometric Functions,  Up: Trigonometric Functions

7.31.4 Conversion Functions
---------------------------

 -- Function: int gsl_sf_polar_to_rect (double R, double THETA,
          gsl_sf_result * X, gsl_sf_result * Y);
     This function converts the polar coordinates (R,THETA) to
     rectilinear coordinates (X,Y), x = r\cos(\theta), y =
     r\sin(\theta).

 -- Function: int gsl_sf_rect_to_polar (double X, double Y,
          gsl_sf_result * R, gsl_sf_result * 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].

\1f
File: gsl-ref.info,  Node: Restriction Functions,  Next: Trigonometric Functions With Error Estimates,  Prev: Conversion Functions,  Up: Trigonometric Functions

7.31.5 Restriction Functions
----------------------------

 -- Function: double gsl_sf_angle_restrict_symm (double THETA)
 -- Function: int gsl_sf_angle_restrict_symm_e (double * THETA)
     These routines force the angle THETA to lie in the range
     (-\pi,\pi].

     Note that the mathematical value of \pi is slightly greater than
     `M_PI', so the machine numbers `M_PI' and `-M_PI' are included in
     the range.

 -- Function: double gsl_sf_angle_restrict_pos (double THETA)
 -- Function: int gsl_sf_angle_restrict_pos_e (double * THETA)
     These routines force the angle THETA to lie in the range [0, 2\pi).

     Note that the mathematical value of 2\pi is slightly greater than
     `2*M_PI', so the machine number `2*M_PI' is included in the range.


\1f
File: gsl-ref.info,  Node: Trigonometric Functions With Error Estimates,  Prev: Restriction Functions,  Up: Trigonometric Functions

7.31.6 Trigonometric Functions With Error Estimates
---------------------------------------------------

 -- Function: int gsl_sf_sin_err_e (double X, double DX, gsl_sf_result
          * RESULT)
     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.

 -- Function: int gsl_sf_cos_err_e (double X, double DX, gsl_sf_result
          * RESULT)
     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.

\1f
File: gsl-ref.info,  Node: Zeta Functions,  Next: Special Functions Examples,  Prev: Trigonometric Functions,  Up: Special Functions

7.32 Zeta Functions
===================

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'.

* Menu:

* Riemann Zeta Function::
* Riemann Zeta Function Minus One::
* Hurwitz Zeta Function::
* Eta Function::

\1f
File: gsl-ref.info,  Node: Riemann Zeta Function,  Next: Riemann Zeta Function Minus One,  Up: Zeta Functions

7.32.1 Riemann Zeta Function
----------------------------

The Riemann zeta function is defined by the infinite sum \zeta(s) =
\sum_{k=1}^\infty k^{-s}.

 -- Function: double gsl_sf_zeta_int (int N)
 -- Function: int gsl_sf_zeta_int_e (int N, gsl_sf_result * RESULT)
     These routines compute the Riemann zeta function \zeta(n) for
     integer N, n \ne 1.

 -- Function: double gsl_sf_zeta (double S)
 -- Function: int gsl_sf_zeta_e (double S, gsl_sf_result * RESULT)
     These routines compute the Riemann zeta function \zeta(s) for
     arbitrary S, s \ne 1.

\1f
File: gsl-ref.info,  Node: Riemann Zeta Function Minus One,  Next: Hurwitz Zeta Function,  Prev: Riemann Zeta Function,  Up: Zeta Functions

7.32.2 Riemann Zeta Function Minus One
--------------------------------------

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.

 -- Function: double gsl_sf_zetam1_int (int N)
 -- Function: int gsl_sf_zetam1_int_e (int N, gsl_sf_result * RESULT)
     These routines compute \zeta(n) - 1 for integer N, n \ne 1.

 -- Function: double gsl_sf_zetam1 (double S)
 -- Function: int gsl_sf_zetam1_e (double S, gsl_sf_result * RESULT)
     These routines compute \zeta(s) - 1 for arbitrary S, s \ne 1.

\1f
File: gsl-ref.info,  Node: Hurwitz Zeta Function,  Next: Eta Function,  Prev: Riemann Zeta Function Minus One,  Up: Zeta Functions

7.32.3 Hurwitz Zeta Function
----------------------------

The Hurwitz zeta function is defined by \zeta(s,q) = \sum_0^\infty
(k+q)^{-s}.

 -- Function: double gsl_sf_hzeta (double S, double Q)
 -- Function: int gsl_sf_hzeta_e (double S, double Q, gsl_sf_result *
          RESULT)
     These routines compute the Hurwitz zeta function \zeta(s,q) for s
     > 1, q > 0.

\1f
File: gsl-ref.info,  Node: Eta Function,  Prev: Hurwitz Zeta Function,  Up: Zeta Functions

7.32.4 Eta Function
-------------------

The eta function is defined by \eta(s) = (1-2^{1-s}) \zeta(s).

 -- Function: double gsl_sf_eta_int (int N)
 -- Function: int gsl_sf_eta_int_e (int N, gsl_sf_result * RESULT)
     These routines compute the eta function \eta(n) for integer N.

 -- Function: double gsl_sf_eta (double S)
 -- Function: int gsl_sf_eta_e (double S, gsl_sf_result * RESULT)
     These routines compute the eta function \eta(s) for arbitrary S.

\1f
File: gsl-ref.info,  Node: Special Functions Examples,  Next: Special Functions References and Further Reading,  Prev: Zeta Functions,  Up: Special Functions

7.33 Examples
=============

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

\1f
File: gsl-ref.info,  Node: Special Functions References and Further Reading,  Prev: Special Functions Examples,  Up: Special Functions

7.34 References and Further Reading
===================================

The library follows the conventions of `Abramowitz & Stegun' where
possible,
     Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions'

The following papers contain information on the algorithms used to
compute the special functions, 
     MISCFUN: A software package to compute uncommon special functions.
     `ACM Trans. Math. Soft.', vol. 22, 1996, 288-301

     G.N. Watson, A Treatise on the Theory of Bessel Functions, 2nd
     Edition (Cambridge University Press, 1944).

     G. Nemeth, Mathematical Approximations of Special Functions, Nova
     Science Publishers, ISBN 1-56072-052-2

     B.C. Carlson, Special Functions of Applied Mathematics (1977)

     W.J. Thompson, Atlas for Computing Mathematical Functions, John
     Wiley & Sons, New York (1997).

     Y.Y. Luke, Algorithms for the Computation of Mathematical
     Functions, Academic Press, New York (1977).


\1f
File: gsl-ref.info,  Node: Vectors and Matrices,  Next: Permutations,  Prev: Special Functions,  Up: Top

8 Vectors and Matrices
**********************

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.

* Menu:

* Data types::
* Blocks::
* Vectors::
* Matrices::
* Vector and Matrix References and Further Reading::

\1f
File: gsl-ref.info,  Node: Data types,  Next: Blocks,  Up: Vectors and Matrices

8.1 Data types
==============

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.

\1f
File: gsl-ref.info,  Node: Blocks,  Next: Vectors,  Prev: Data types,  Up: Vectors and Matrices

8.2 Blocks
==========

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'

* Menu:

* Block allocation::
* Reading and writing blocks::
* Example programs for blocks::

\1f
File: gsl-ref.info,  Node: Block allocation,  Next: Reading and writing blocks,  Up: Blocks

8.2.1 Block allocation
----------------------

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'.

 -- Function: gsl_block * gsl_block_alloc (size_t N)
     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_calloc' if 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.

 -- Function: gsl_block * gsl_block_calloc (size_t N)
     This function allocates memory for a block and initializes all the
     elements of the block to zero.

 -- Function: void gsl_block_free (gsl_block * B)
     This function frees the memory used by a block B previously
     allocated with `gsl_block_alloc' or `gsl_block_calloc'.  The block
     B must be a valid block object (a null pointer is not allowed).

\1f
File: gsl-ref.info,  Node: Reading and writing blocks,  Next: Example programs for blocks,  Prev: Block allocation,  Up: Blocks

8.2.2 Reading and writing blocks
--------------------------------

The library provides functions for reading and writing blocks to a file
as binary data or formatted text.

 -- Function: int gsl_block_fwrite (FILE * STREAM, const gsl_block * B)
     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_EFAILED' if 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.

 -- Function: int gsl_block_fread (FILE * STREAM, gsl_block * B)
     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_EFAILED' if 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.

 -- Function: int gsl_block_fprintf (FILE * STREAM, const gsl_block *
          B, const char * FORMAT)
     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', `%e' or `%f' formats for floating point
     numbers and `%d' for integers.  The function returns 0 for success
     and `GSL_EFAILED' if there was a problem writing to the file.

 -- Function: int gsl_block_fscanf (FILE * STREAM, gsl_block * B)
     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_EFAILED' if there was a problem reading from the file.

\1f
File: gsl-ref.info,  Node: Example programs for blocks,  Prev: Reading and writing blocks,  Up: Blocks

8.2.3 Example programs for blocks
---------------------------------

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

\1f
File: gsl-ref.info,  Node: Vectors,  Next: Matrices,  Prev: Blocks,  Up: Vectors and Matrices

8.3 Vectors
===========

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'

* Menu:

* Vector allocation::
* Accessing vector elements::
* Initializing vector elements::
* Reading and writing vectors::
* Vector views::
* Copying vectors::
* Exchanging elements::
* Vector operations::
* Finding maximum and minimum elements of vectors::
* Vector properties::
* Example programs for vectors::

\1f
File: gsl-ref.info,  Node: Vector allocation,  Next: Accessing vector elements,  Up: Vectors

8.3.1 Vector allocation
-----------------------

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'.

 -- Function: gsl_vector * gsl_vector_alloc (size_t N)
     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.

 -- Function: gsl_vector * gsl_vector_calloc (size_t N)
     This function allocates memory for a vector of length N and
     initializes all the elements of the vector to zero.

 -- Function: void gsl_vector_free (gsl_vector * V)
     This function frees a previously allocated vector V.  If the
     vector was created using `gsl_vector_alloc' then 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.  The vector V must be a
     valid vector object (a null pointer is not allowed).

\1f
File: gsl-ref.info,  Node: Accessing vector elements,  Next: Initializing vector elements,  Prev: Vector allocation,  Up: Vectors

8.3.2 Accessing vector elements
-------------------------------

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.

 -- Function: double gsl_vector_get (const gsl_vector * V, size_t I)
     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.  An inline version of this function is
     used when `HAVE_INLINE' is defined.

 -- Function: void gsl_vector_set (gsl_vector * V, size_t I, double X)
     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.  An inline version of this function is used
     when `HAVE_INLINE' is defined.

 -- Function: double * gsl_vector_ptr (gsl_vector * V, size_t I)
 -- Function: const double * gsl_vector_const_ptr (const gsl_vector *
          V, size_t I)
     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.  Inline
     versions of these functions are used when `HAVE_INLINE' is defined.

\1f
File: gsl-ref.info,  Node: Initializing vector elements,  Next: Reading and writing vectors,  Prev: Accessing vector elements,  Up: Vectors

8.3.3 Initializing vector elements
----------------------------------

 -- Function: void gsl_vector_set_all (gsl_vector * V, double X)
     This function sets all the elements of the vector V to the value X.

 -- Function: void gsl_vector_set_zero (gsl_vector * V)
     This function sets all the elements of the vector V to zero.

 -- Function: int gsl_vector_set_basis (gsl_vector * V, size_t I)
     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.

\1f
File: gsl-ref.info,  Node: Reading and writing vectors,  Next: Vector views,  Prev: Initializing vector elements,  Up: Vectors

8.3.4 Reading and writing vectors
---------------------------------

The library provides functions for reading and writing vectors to a file
as binary data or formatted text.

 -- Function: int gsl_vector_fwrite (FILE * STREAM, const gsl_vector *
          V)
     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_EFAILED' if 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.

 -- Function: int gsl_vector_fread (FILE * STREAM, gsl_vector * V)
     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_EFAILED' if 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.

 -- Function: int gsl_vector_fprintf (FILE * STREAM, const gsl_vector *
          V, const char * FORMAT)
     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', `%e' or `%f' formats for floating point
     numbers and `%d' for integers.  The function returns 0 for success
     and `GSL_EFAILED' if there was a problem writing to the file.

 -- Function: int gsl_vector_fscanf (FILE * STREAM, gsl_vector * V)
     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_EFAILED' if there was a problem reading from the file.

\1f
File: gsl-ref.info,  Node: Vector views,  Next: Copying vectors,  Prev: Reading and writing vectors,  Up: Vectors

8.3.5 Vector views
------------------

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.

 -- Function: gsl_vector_view gsl_vector_subvector (gsl_vector * V,
          size_t OFFSET, size_t N)
 -- Function: gsl_vector_const_view gsl_vector_const_subvector (const
          gsl_vector * V, size_t OFFSET, size_t N)
     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 `data' pointer 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_subvector' is equivalent to
     `gsl_vector_subvector' but can be used for vectors which are
     declared `const'.

 -- Function: gsl_vector_view gsl_vector_subvector_with_stride
          (gsl_vector * V, size_t OFFSET, size_t STRIDE, size_t N)
 -- Function: gsl_vector_const_view
gsl_vector_const_subvector_with_stride (const gsl_vector * V, size_t
          OFFSET, size_t STRIDE, size_t N)
     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_subvector' but 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 `v' of length `n', 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 of `v' using 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_stride' is equivalent
     to `gsl_vector_subvector_with_stride' but can be used for vectors
     which are declared `const'.

 -- Function: gsl_vector_view gsl_vector_complex_real
          (gsl_vector_complex * V)
 -- Function: gsl_vector_const_view gsl_vector_complex_const_real
          (const gsl_vector_complex * V)
     These functions return a vector view of the real parts of the
     complex vector V.

     The function `gsl_vector_complex_const_real' is equivalent to
     `gsl_vector_complex_real' but can be used for vectors which are
     declared `const'.

 -- Function: gsl_vector_view gsl_vector_complex_imag
          (gsl_vector_complex * V)
 -- Function: gsl_vector_const_view gsl_vector_complex_const_imag
          (const gsl_vector_complex * V)
     These functions return a vector view of the imaginary parts of the
     complex vector V.

     The function `gsl_vector_complex_const_imag' is equivalent to
     `gsl_vector_complex_imag' but can be used for vectors which are
     declared `const'.

 -- Function: gsl_vector_view gsl_vector_view_array (double * BASE,
          size_t N)
 -- Function: gsl_vector_const_view gsl_vector_const_view_array (const
          double * BASE, size_t N)
     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_array' is equivalent to
     `gsl_vector_view_array' but can be used for arrays which are
     declared `const'.

 -- Function: gsl_vector_view gsl_vector_view_array_with_stride (double
          * BASE, size_t STRIDE, size_t N)
 -- Function: gsl_vector_const_view
gsl_vector_const_view_array_with_stride (const double * BASE, size_t
          STRIDE, size_t N)
     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_array' but 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_stride' is
     equivalent to `gsl_vector_view_array_with_stride' but can be used
     for arrays which are declared `const'.

\1f
File: gsl-ref.info,  Node: Copying vectors,  Next: Exchanging elements,  Prev: Vector views,  Up: Vectors

8.3.6 Copying vectors
---------------------

Common operations on vectors such as addition and multiplication are
available in the BLAS part of the library (*note 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.

 -- Function: int gsl_vector_memcpy (gsl_vector * DEST, const
          gsl_vector * SRC)
     This function copies the elements of the vector SRC into the
     vector DEST.  The two vectors must have the same length.

 -- Function: int gsl_vector_swap (gsl_vector * V, gsl_vector * W)
     This function exchanges the elements of the vectors V and W by
     copying.  The two vectors must have the same length.

\1f
File: gsl-ref.info,  Node: Exchanging elements,  Next: Vector operations,  Prev: Copying vectors,  Up: Vectors

8.3.7 Exchanging elements
-------------------------

The following function can be used to exchange, or permute, the elements
of a vector.

 -- Function: int gsl_vector_swap_elements (gsl_vector * V, size_t I,
          size_t J)
     This function exchanges the I-th and J-th elements of the vector V
     in-place.

 -- Function: int gsl_vector_reverse (gsl_vector * V)
     This function reverses the order of the elements of the vector V.

\1f
File: gsl-ref.info,  Node: Vector operations,  Next: Finding maximum and minimum elements of vectors,  Prev: Exchanging elements,  Up: Vectors

8.3.8 Vector operations
-----------------------

The following operations are only defined for real vectors.

 -- Function: int gsl_vector_add (gsl_vector * A, const gsl_vector * B)
     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.

 -- Function: int gsl_vector_sub (gsl_vector * A, const gsl_vector * B)
     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.

 -- Function: int gsl_vector_mul (gsl_vector * A, const gsl_vector * B)
     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.

 -- Function: int gsl_vector_div (gsl_vector * A, const gsl_vector * B)
     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.

 -- Function: int gsl_vector_scale (gsl_vector * A, const double X)
     This function multiplies the elements of vector A by the constant
     factor X, a'_i = x a_i.

 -- Function: int gsl_vector_add_constant (gsl_vector * A, const double
          X)
     This function adds the constant value X to the elements of the
     vector A, a'_i = a_i + x.

\1f
File: gsl-ref.info,  Node: Finding maximum and minimum elements of vectors,  Next: Vector properties,  Prev: Vector operations,  Up: Vectors

8.3.9 Finding maximum and minimum elements of vectors
-----------------------------------------------------

 -- Function: double gsl_vector_max (const gsl_vector * V)
     This function returns the maximum value in the vector V.

 -- Function: double gsl_vector_min (const gsl_vector * V)
     This function returns the minimum value in the vector V.

 -- Function: void gsl_vector_minmax (const gsl_vector * V, double *
          MIN_OUT, double * MAX_OUT)
     This function returns the minimum and maximum values in the vector
     V, storing them in MIN_OUT and MAX_OUT.

 -- Function: size_t gsl_vector_max_index (const gsl_vector * V)
     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.

 -- Function: size_t gsl_vector_min_index (const gsl_vector * V)
     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.

 -- Function: void gsl_vector_minmax_index (const gsl_vector * V,
          size_t * IMIN, size_t * IMAX)
     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.

\1f
File: gsl-ref.info,  Node: Vector properties,  Next: Example programs for vectors,  Prev: Finding maximum and minimum elements of vectors,  Up: Vectors

8.3.10 Vector properties
------------------------

 -- Function: int gsl_vector_isnull (const gsl_vector * V)
 -- Function: int gsl_vector_ispos (const gsl_vector * V)
 -- Function: int gsl_vector_isneg (const gsl_vector * V)
 -- Function: int gsl_vector_isnonneg (const gsl_vector * V)
     These functions return 1 if all the elements of the vector V are
     zero, strictly positive, strictly negative, or non-negative
     respectively, and 0 otherwise.

\1f
File: gsl-ref.info,  Node: Example programs for vectors,  Prev: Vector properties,  Up: Vectors

8.3.11 Example programs for vectors
-----------------------------------

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++) /* OUT OF RANGE ERROR */
         {
           printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
         }

       gsl_vector_free (v);
       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);
       }

       gsl_vector_free (v);
       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));
         }

       gsl_vector_free (v);
       return 0;
     }

\1f
File: gsl-ref.info,  Node: Matrices,  Next: Vector and Matrix References and Further Reading,  Prev: Vectors,  Up: Vectors and Matrices

8.4 Matrices
============

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'

* Menu:

* Matrix allocation::
* Accessing matrix elements::
* Initializing matrix elements::
* Reading and writing matrices::
* Matrix views::
* Creating row and column views::
* Copying matrices::
* Copying rows and columns::
* Exchanging rows and columns::
* Matrix operations::
* Finding maximum and minimum elements of matrices::
* Matrix properties::
* Example programs for matrices::

\1f
File: gsl-ref.info,  Node: Matrix allocation,  Next: Accessing matrix elements,  Up: Matrices

8.4.1 Matrix allocation
-----------------------

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 matrix 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'.

 -- Function: gsl_matrix * gsl_matrix_alloc (size_t N1, size_t N2)
     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.

 -- Function: gsl_matrix * gsl_matrix_calloc (size_t N1, size_t N2)
     This function allocates memory for a matrix of size N1 rows by N2
     columns and initializes all the elements of the matrix to zero.

 -- Function: void gsl_matrix_free (gsl_matrix * M)
     This function frees a previously allocated matrix M.  If the
     matrix was created using `gsl_matrix_alloc' then 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 matrix M must be a
     valid matrix object (a null pointer is not allowed).

\1f
File: gsl-ref.info,  Node: Accessing matrix elements,  Next: Initializing matrix elements,  Prev: Matrix allocation,  Up: Matrices

8.4.2 Accessing matrix elements
-------------------------------

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.

 -- Function: double gsl_matrix_get (const gsl_matrix * M, size_t I,
          size_t J)
     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.  An inline version
     of this function is used when `HAVE_INLINE' is defined.

 -- Function: void gsl_matrix_set (gsl_matrix * M, size_t I, size_t J,
          double X)
     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.  An inline version of
     this function is used when `HAVE_INLINE' is defined.

 -- Function: double * gsl_matrix_ptr (gsl_matrix * M, size_t I, size_t
          J)
 -- Function: const double * gsl_matrix_const_ptr (const gsl_matrix *
          M, size_t I, size_t J)
     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.  Inline versions of these functions are used when
     `HAVE_INLINE' is defined.

\1f
File: gsl-ref.info,  Node: Initializing matrix elements,  Next: Reading and writing matrices,  Prev: Accessing matrix elements,  Up: Matrices

8.4.3 Initializing matrix elements
----------------------------------

 -- Function: void gsl_matrix_set_all (gsl_matrix * M, double X)
     This function sets all the elements of the matrix M to the value X.

 -- Function: void gsl_matrix_set_zero (gsl_matrix * M)
     This function sets all the elements of the matrix M to zero.

 -- Function: void gsl_matrix_set_identity (gsl_matrix * M)
     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.

\1f
File: gsl-ref.info,  Node: Reading and writing matrices,  Next: Matrix views,  Prev: Initializing matrix elements,  Up: Matrices

8.4.4 Reading and writing matrices
----------------------------------

The library provides functions for reading and writing matrices to a
file as binary data or formatted text.

 -- Function: int gsl_matrix_fwrite (FILE * STREAM, const gsl_matrix *
          M)
     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_EFAILED' if 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.

 -- Function: int gsl_matrix_fread (FILE * STREAM, gsl_matrix * M)
     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_EFAILED' if 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.

 -- Function: int gsl_matrix_fprintf (FILE * STREAM, const gsl_matrix *
          M, const char * FORMAT)
     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', `%e' or `%f' formats for floating point
     numbers and `%d' for integers.  The function returns 0 for success
     and `GSL_EFAILED' if there was a problem writing to the file.

 -- Function: int gsl_matrix_fscanf (FILE * STREAM, gsl_matrix * M)
     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_EFAILED' if there was a problem reading from the file.

\1f
File: gsl-ref.info,  Node: Matrix views,  Next: Creating row and column views,  Prev: Reading and writing matrices,  Up: Matrices

8.4.5 Matrix views
------------------

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.

 -- Function: gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * M,
          size_t K1, size_t K2, size_t N1, size_t N2)
 -- Function: gsl_matrix_const_view gsl_matrix_const_submatrix (const
          gsl_matrix * M, size_t K1, size_t K2, size_t N1, size_t N2)
     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-1' and the index J runs from
     0 to `n2-1'.

     The `data' pointer 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_submatrix' is equivalent to
     `gsl_matrix_submatrix' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_matrix_view gsl_matrix_view_array (double * BASE,
          size_t N1, size_t N2)
 -- Function: gsl_matrix_const_view gsl_matrix_const_view_array (const
          double * BASE, size_t N1, size_t N2)
     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-1' and the index J runs from
     0 to `n2-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' is equivalent to
     `gsl_matrix_view_array' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_matrix_view gsl_matrix_view_array_with_tda (double *
          BASE, size_t N1, size_t N2, size_t TDA)
 -- Function: gsl_matrix_const_view
          gsl_matrix_const_view_array_with_tda (const double * BASE,
          size_t N1, size_t N2, size_t TDA)
     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-1' and the index J runs from
     0 to `n2-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_tda' is equivalent
     to `gsl_matrix_view_array_with_tda' but can be used for matrices
     which are declared `const'.

 -- Function: gsl_matrix_view gsl_matrix_view_vector (gsl_vector * V,
          size_t N1, size_t N2)
 -- Function: gsl_matrix_const_view gsl_matrix_const_view_vector (const
          gsl_vector * V, size_t N1, size_t N2)
     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-1' and the index J runs from
     0 to `n2-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' is equivalent to
     `gsl_matrix_view_vector' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_matrix_view gsl_matrix_view_vector_with_tda
          (gsl_vector * V, size_t N1, size_t N2, size_t TDA)
 -- Function: gsl_matrix_const_view
gsl_matrix_const_view_vector_with_tda (const gsl_vector * V, size_t N1,
          size_t N2, size_t TDA)
     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-1' and the index J runs from
     0 to `n2-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_tda' is equivalent
     to `gsl_matrix_view_vector_with_tda' but can be used for matrices
     which are declared `const'.

\1f
File: gsl-ref.info,  Node: Creating row and column views,  Next: Copying matrices,  Prev: Matrix views,  Up: Matrices

8.4.6 Creating row and column views
-----------------------------------

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.

 -- Function: gsl_vector_view gsl_matrix_row (gsl_matrix * M, size_t I)
 -- Function: gsl_vector_const_view gsl_matrix_const_row (const
          gsl_matrix * M, size_t I)
     These functions return a vector view of the I-th row of the matrix
     M.  The `data' pointer of the new vector is set to null if I is
     out of range.

     The function `gsl_vector_const_row' is equivalent to
     `gsl_matrix_row' but can be used for matrices which are declared
     `const'.

 -- Function: gsl_vector_view gsl_matrix_column (gsl_matrix * M, size_t
          J)
 -- Function: gsl_vector_const_view gsl_matrix_const_column (const
          gsl_matrix * M, size_t J)
     These functions return a vector view of the J-th column of the
     matrix M.  The `data' pointer of the new vector is set to null if
     J is out of range.

     The function `gsl_vector_const_column' is equivalent to
     `gsl_matrix_column' but can be used for matrices which are declared
     `const'.

 -- Function: gsl_vector_view gsl_matrix_subrow (gsl_matrix * M, size_t
          I, size_t OFFSET, size_t N)
 -- Function: gsl_vector_const_view gsl_matrix_const_subrow (const
          gsl_matrix * M, size_t I, size_t OFFSET, size_t N)
     These functions return a vector view of the I-th row of the matrix
     M beginning at OFFSET elements past the first column and
     containing N elements. The `data' pointer of the new vector is set
     to null if I, OFFSET, or N are out of range.

     The function `gsl_vector_const_subrow' is equivalent to
     `gsl_matrix_subrow' but can be used for matrices which are declared
     `const'.

 -- Function: gsl_vector_view gsl_matrix_subcolumn (gsl_matrix * M,
          size_t J, size_t OFFSET, size_t N)
 -- Function: gsl_vector_const_view gsl_matrix_const_subcolumn (const
          gsl_matrix * M, size_t J, size_t OFFSET, size_t N)
     These functions return a vector view of the J-th column of the
     matrix M beginning at OFFSET elements past the first row and
     containing N elements. The `data' pointer of the new vector is set
     to null if J, OFFSET, or N are out of range.

     The function `gsl_vector_const_subcolumn' is equivalent to
     `gsl_matrix_subcolumn' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_vector_view gsl_matrix_diagonal (gsl_matrix * M)
 -- Function: gsl_vector_const_view gsl_matrix_const_diagonal (const
          gsl_matrix * M)
     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_diagonal' is equivalent to
     `gsl_matrix_diagonal' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * M,
          size_t K)
 -- Function: gsl_vector_const_view gsl_matrix_const_subdiagonal (const
          gsl_matrix * M, size_t K)
     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_subdiagonal' is equivalent to
     `gsl_matrix_subdiagonal' but can be used for matrices which are
     declared `const'.

 -- Function: gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * M,
          size_t K)
 -- Function: gsl_vector_const_view gsl_matrix_const_superdiagonal
          (const gsl_matrix * M, size_t K)
     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_superdiagonal' is equivalent to
     `gsl_matrix_superdiagonal' but can be used for matrices which are
     declared `const'.

\1f
File: gsl-ref.info,  Node: Copying matrices,  Next: Copying rows and columns,  Prev: Creating row and column views,  Up: Matrices

8.4.7 Copying matrices
----------------------

 -- Function: int gsl_matrix_memcpy (gsl_matrix * DEST, const
          gsl_matrix * SRC)
     This function copies the elements of the matrix SRC into the
     matrix DEST.  The two matrices must have the same size.

 -- Function: int gsl_matrix_swap (gsl_matrix * M1, gsl_matrix * M2)
     This function exchanges the elements of the matrices M1 and M2 by
     copying.  The two matrices must have the same size.

\1f
File: gsl-ref.info,  Node: Copying rows and columns,  Next: Exchanging rows and columns,  Prev: Copying matrices,  Up: Matrices

8.4.8 Copying rows and columns
------------------------------

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.

 -- Function: int gsl_matrix_get_row (gsl_vector * V, const gsl_matrix
          * M, size_t I)
     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.

 -- Function: int gsl_matrix_get_col (gsl_vector * V, const gsl_matrix
          * M, size_t J)
     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.

 -- Function: int gsl_matrix_set_row (gsl_matrix * M, size_t I, const
          gsl_vector * V)
     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.

 -- Function: int gsl_matrix_set_col (gsl_matrix * M, size_t J, const
          gsl_vector * V)
     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.

\1f
File: gsl-ref.info,  Node: Exchanging rows and columns,  Next: Matrix operations,  Prev: Copying rows and columns,  Up: Matrices

8.4.9 Exchanging rows and columns
---------------------------------

The following functions can be used to exchange the rows and columns of
a matrix.

 -- Function: int gsl_matrix_swap_rows (gsl_matrix * M, size_t I,
          size_t J)
     This function exchanges the I-th and J-th rows of the matrix M
     in-place.

 -- Function: int gsl_matrix_swap_columns (gsl_matrix * M, size_t I,
          size_t J)
     This function exchanges the I-th and J-th columns of the matrix M
     in-place.

 -- Function: int gsl_matrix_swap_rowcol (gsl_matrix * M, size_t I,
          size_t J)
     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.

 -- Function: int gsl_matrix_transpose_memcpy (gsl_matrix * DEST, const
          gsl_matrix * SRC)
     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.

 -- Function: int gsl_matrix_transpose (gsl_matrix * M)
     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.

\1f
File: gsl-ref.info,  Node: Matrix operations,  Next: Finding maximum and minimum elements of matrices,  Prev: Exchanging rows and columns,  Up: Matrices

8.4.10 Matrix operations
------------------------

The following operations are defined for real and complex matrices.

 -- Function: int gsl_matrix_add (gsl_matrix * A, const gsl_matrix * B)
     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.

 -- Function: int gsl_matrix_sub (gsl_matrix * A, const gsl_matrix * B)
     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.

 -- Function: int gsl_matrix_mul_elements (gsl_matrix * A, const
          gsl_matrix * B)
     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.

 -- Function: int gsl_matrix_div_elements (gsl_matrix * A, const
          gsl_matrix * B)
     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.

 -- Function: int gsl_matrix_scale (gsl_matrix * A, const double X)
     This function multiplies the elements of matrix A by the constant
     factor X, a'(i,j) = x a(i,j).

 -- Function: int gsl_matrix_add_constant (gsl_matrix * A, const double
          X)
     This function adds the constant value X to the elements of the
     matrix A, a'(i,j) = a(i,j) + x.

\1f
File: gsl-ref.info,  Node: Finding maximum and minimum elements of matrices,  Next: Matrix properties,  Prev: Matrix operations,  Up: Matrices

8.4.11 Finding maximum and minimum elements of matrices
-------------------------------------------------------

The following operations are only defined for real matrices.

 -- Function: double gsl_matrix_max (const gsl_matrix * M)
     This function returns the maximum value in the matrix M.

 -- Function: double gsl_matrix_min (const gsl_matrix * M)
     This function returns the minimum value in the matrix M.

 -- Function: void gsl_matrix_minmax (const gsl_matrix * M, double *
          MIN_OUT, double * MAX_OUT)
     This function returns the minimum and maximum values in the matrix
     M, storing them in MIN_OUT and MAX_OUT.

 -- Function: void gsl_matrix_max_index (const gsl_matrix * M, size_t *
          IMAX, size_t * JMAX)
     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.

 -- Function: void gsl_matrix_min_index (const gsl_matrix * M, size_t *
          IMIN, size_t * JMIN)
     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.

 -- Function: void gsl_matrix_minmax_index (const gsl_matrix * M,
          size_t * IMIN, size_t * JMIN, size_t * IMAX, size_t * JMAX)
     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.

\1f
File: gsl-ref.info,  Node: Matrix properties,  Next: Example programs for matrices,  Prev: Finding maximum and minimum elements of matrices,  Up: Matrices

8.4.12 Matrix properties
------------------------

 -- Function: int gsl_matrix_isnull (const gsl_matrix * M)
 -- Function: int gsl_matrix_ispos (const gsl_matrix * M)
 -- Function: int gsl_matrix_isneg (const gsl_matrix * M)
 -- Function: int gsl_matrix_isnonneg (const gsl_matrix * M)
     These functions return 1 if all the elements of the matrix M are
     zero, strictly positive, strictly negative, or non-negative
     respectively, and 0 otherwise. To test whether a matrix is
     positive-definite, use the Cholesky decomposition (*note Cholesky
     Decomposition::).

\1f
File: gsl-ref.info,  Node: Example programs for matrices,  Prev: Matrix properties,  Up: Matrices

8.4.13 Example programs for matrices
------------------------------------

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++)  /* OUT OF RANGE ERROR */
         for (j = 0; j < 3; j++)
           printf ("m(%d,%d) = %g\n", i, j,
                   gsl_matrix_get (m, i, j));

       gsl_matrix_free (m);

       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++;
           }

       gsl_matrix_free (m);
       gsl_matrix_free (a);

       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

\1f
File: gsl-ref.info,  Node: Vector and Matrix References and Further Reading,  Prev: Matrices,  Up: Vectors and Matrices

8.5 References and Further Reading
==================================

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,

     B. Stroustrup, `The C++ Programming Language' (3rd Ed), Section
     22.4 Vector Arithmetic.  Addison-Wesley 1997, ISBN 0-201-88954-4.

\1f
File: gsl-ref.info,  Node: Permutations,  Next: Combinations,  Prev: Vectors and Matrices,  Up: Top

9 Permutations
**************

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'.

* Menu:

* The Permutation struct::
* Permutation allocation::
* Accessing permutation elements::
* Permutation properties::
* Permutation functions::
* Applying Permutations::
* Reading and writing permutations::
* Permutations in cyclic form::
* Permutation Examples::
* Permutation References and Further Reading::

\1f
File: gsl-ref.info,  Node: The Permutation struct,  Next: Permutation allocation,  Up: Permutations

9.1 The Permutation struct
==========================

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;


\1f
File: gsl-ref.info,  Node: Permutation allocation,  Next: Accessing permutation elements,  Prev: The Permutation struct,  Up: Permutations

9.2 Permutation allocation
==========================

 -- Function: gsl_permutation * gsl_permutation_alloc (size_t N)
     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_calloc' if 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.

 -- Function: gsl_permutation * gsl_permutation_calloc (size_t N)
     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.

 -- Function: void gsl_permutation_init (gsl_permutation * P)
     This function initializes the permutation P to the identity, i.e.
     (0,1,2,...,n-1).

 -- Function: void gsl_permutation_free (gsl_permutation * P)
     This function frees all the memory used by the permutation P.

 -- Function: int gsl_permutation_memcpy (gsl_permutation * DEST, const
          gsl_permutation * SRC)
     This function copies the elements of the permutation SRC into the
     permutation DEST.  The two permutations must have the same size.

\1f
File: gsl-ref.info,  Node: Accessing permutation elements,  Next: Permutation properties,  Prev: Permutation allocation,  Up: Permutations

9.3 Accessing permutation elements
==================================

The following functions can be used to access and manipulate
permutations.

 -- Function: size_t gsl_permutation_get (const gsl_permutation * P,
          const size_t I)
     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.  An inline
     version of this function is used when `HAVE_INLINE' is defined.

 -- Function: int gsl_permutation_swap (gsl_permutation * P, const
          size_t I, const size_t J)
     This function exchanges the I-th and J-th elements of the
     permutation P.

\1f
File: gsl-ref.info,  Node: Permutation properties,  Next: Permutation functions,  Prev: Accessing permutation elements,  Up: Permutations

9.4 Permutation properties
==========================

 -- Function: size_t gsl_permutation_size (const gsl_permutation * P)
     This function returns the size of the permutation P.

 -- Function: size_t * gsl_permutation_data (const gsl_permutation * P)
     This function returns a pointer to the array of elements in the
     permutation P.

 -- Function: int gsl_permutation_valid (const gsl_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.

\1f
File: gsl-ref.info,  Node: Permutation functions,  Next: Applying Permutations,  Prev: Permutation properties,  Up: Permutations

9.5 Permutation functions
=========================

 -- Function: void gsl_permutation_reverse (gsl_permutation * P)
     This function reverses the elements of the permutation P.

 -- Function: int gsl_permutation_inverse (gsl_permutation * INV, const
          gsl_permutation * P)
     This function computes the inverse of the permutation P, storing
     the result in INV.

 -- Function: int gsl_permutation_next (gsl_permutation * P)
     This function advances the permutation P to the next permutation
     in lexicographic order and returns `GSL_SUCCESS'.  If no further
     permutations are available it returns `GSL_FAILURE' and leaves P
     unmodified.  Starting with the identity permutation and repeatedly
     applying this function will iterate through all possible
     permutations of a given order.

 -- Function: int gsl_permutation_prev (gsl_permutation * P)
     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 returns
     `GSL_FAILURE' and leaves P unmodified.

\1f
File: gsl-ref.info,  Node: Applying Permutations,  Next: Reading and writing permutations,  Prev: Permutation functions,  Up: Permutations

9.6 Applying Permutations
=========================

 -- Function: int gsl_permute (const size_t * P, double * DATA, size_t
          STRIDE, size_t N)
     This function applies the permutation P to the array DATA of size
     N with stride STRIDE.

 -- Function: int gsl_permute_inverse (const size_t * P, double * DATA,
          size_t STRIDE, size_t N)
     This function applies the inverse of the permutation P to the
     array DATA of size N with stride STRIDE.

 -- Function: int gsl_permute_vector (const gsl_permutation * P,
          gsl_vector * V)
     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.

 -- Function: int gsl_permute_vector_inverse (const gsl_permutation *
          P, gsl_vector * V)
     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.

 -- Function: int gsl_permutation_mul (gsl_permutation * P, const
          gsl_permutation * PA, const gsl_permutation * PB)
     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.

\1f
File: gsl-ref.info,  Node: Reading and writing permutations,  Next: Permutations in cyclic form,  Prev: Applying Permutations,  Up: Permutations

9.7 Reading and writing permutations
====================================

The library provides functions for reading and writing permutations to a
file as binary data or formatted text.

 -- Function: int gsl_permutation_fwrite (FILE * STREAM, const
          gsl_permutation * P)
     This function writes the elements of the permutation P to the
     stream STREAM in binary format.  The function returns
     `GSL_EFAILED' if 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.

 -- Function: int gsl_permutation_fread (FILE * STREAM, gsl_permutation
          * P)
     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_EFAILED' if 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.

 -- Function: int gsl_permutation_fprintf (FILE * STREAM, const
          gsl_permutation * P, const char * FORMAT)
     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.  In ISO C99
     the type modifier `z' represents `size_t', so `"%zu\n"' is a
     suitable format.(1) The function returns `GSL_EFAILED' if there
     was a problem writing to the file.

 -- Function: int gsl_permutation_fscanf (FILE * STREAM,
          gsl_permutation * P)
     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_EFAILED' if
     there was a problem reading from the file.

   ---------- Footnotes ----------

   (1) In versions of the GNU C library prior to the ISO C99 standard,
the type modifier `Z' was used instead.

\1f
File: gsl-ref.info,  Node: Permutations in cyclic form,  Next: Permutation Examples,  Prev: Reading and writing permutations,  Up: Permutations

9.8 Permutations in cyclic form
===============================

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,

  1. Write all singleton cycles explicitly

  2. Within each cycle, put the smallest number first

  3. Order the cycles in decreasing order of the first number in the
     cycle.

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.

 -- Function: int gsl_permutation_linear_to_canonical (gsl_permutation
          * Q, const gsl_permutation * P)
     This function computes the canonical form of the permutation P and
     stores it in the output argument Q.

 -- Function: int gsl_permutation_canonical_to_linear (gsl_permutation
          * P, const gsl_permutation * Q)
     This function converts a permutation Q in canonical form back into
     linear form storing it in the output argument P.

 -- Function: size_t gsl_permutation_inversions (const gsl_permutation
          * 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.

 -- Function: size_t gsl_permutation_linear_cycles (const
          gsl_permutation * P)
     This function counts the number of cycles in the permutation P,
     given in linear form.

 -- Function: size_t gsl_permutation_canonical_cycles (const
          gsl_permutation * Q)
     This function counts the number of cycles in the permutation Q,
     given in canonical form.

\1f
File: gsl-ref.info,  Node: Permutation Examples,  Next: Permutation References and Further Reading,  Prev: Permutations in cyclic form,  Up: Permutations

9.9 Examples
============

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");

       gsl_permutation_free (p);
       gsl_permutation_free (q);
       gsl_rng_free (r);

       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);

       gsl_permutation_free (p);

       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'.

\1f
File: gsl-ref.info,  Node: Permutation References and Further Reading,  Prev: Permutation Examples,  Up: Permutations

9.10 References and Further Reading
===================================

The subject of permutations is covered extensively in Knuth's `Sorting
and Searching',

     Donald E. Knuth, `The Art of Computer Programming: Sorting and
     Searching' (Vol 3, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850.

For the definition of the "canonical form" see,

     Donald E. Knuth, `The Art of Computer Programming: Fundamental
     Algorithms' (Vol 1, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850.
     Section 1.3.3, `An Unusual Correspondence', p.178-179.

\1f
File: gsl-ref.info,  Node: Combinations,  Next: Sorting,  Prev: Permutations,  Up: Top

10 Combinations
***************

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'.

* Menu:

* The Combination struct::
* Combination allocation::
* Accessing combination elements::
* Combination properties::
* Combination functions::
* Reading and writing combinations::
* Combination Examples::
* Combination References and Further Reading::

\1f
File: gsl-ref.info,  Node: The Combination struct,  Next: Combination allocation,  Up: Combinations

10.1 The Combination struct
===========================

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;


\1f
File: gsl-ref.info,  Node: Combination allocation,  Next: Accessing combination elements,  Prev: The Combination struct,  Up: Combinations

10.2 Combination allocation
===========================

 -- Function: gsl_combination * gsl_combination_alloc (size_t N, size_t
          K)
     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_calloc'
     if 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.

 -- Function: gsl_combination * gsl_combination_calloc (size_t N,
          size_t K)
     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.

 -- Function: void gsl_combination_init_first (gsl_combination * C)
     This function initializes the combination C to the
     lexicographically first combination, i.e.  (0,1,2,...,k-1).

 -- Function: void gsl_combination_init_last (gsl_combination * C)
     This function initializes the combination C to the
     lexicographically last combination, i.e.  (n-k,n-k+1,...,n-1).

 -- Function: void gsl_combination_free (gsl_combination * C)
     This function frees all the memory used by the combination C.

 -- Function: int gsl_combination_memcpy (gsl_combination * DEST, const
          gsl_combination * SRC)
     This function copies the elements of the combination SRC into the
     combination DEST.  The two combinations must have the same size.

\1f
File: gsl-ref.info,  Node: Accessing combination elements,  Next: Combination properties,  Prev: Combination allocation,  Up: Combinations

10.3 Accessing combination elements
===================================

The following function can be used to access the elements of a
combination.

 -- Function: size_t gsl_combination_get (const gsl_combination * C,
          const size_t I)
     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.  An inline
     version of this function is used when `HAVE_INLINE' is defined.

\1f
File: gsl-ref.info,  Node: Combination properties,  Next: Combination functions,  Prev: Accessing combination elements,  Up: Combinations

10.4 Combination properties
===========================

 -- Function: size_t gsl_combination_n (const gsl_combination * C)
     This function returns the range (n) of the combination C.

 -- Function: size_t gsl_combination_k (const gsl_combination * C)
     This function returns the number of elements (k) in the
     combination C.

 -- Function: size_t * gsl_combination_data (const gsl_combination * C)
     This function returns a pointer to the array of elements in the
     combination C.

 -- Function: int gsl_combination_valid (gsl_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.

\1f
File: gsl-ref.info,  Node: Combination functions,  Next: Reading and writing combinations,  Prev: Combination properties,  Up: Combinations

10.5 Combination functions
==========================

 -- Function: int gsl_combination_next (gsl_combination * C)
     This function advances the combination C to the next combination
     in lexicographic order and returns `GSL_SUCCESS'.  If no further
     combinations are available it returns `GSL_FAILURE' and leaves C
     unmodified.  Starting with the first combination and repeatedly
     applying this function will iterate through all possible
     combinations of a given order.

 -- Function: int gsl_combination_prev (gsl_combination * C)
     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 returns
     `GSL_FAILURE' and leaves C unmodified.

\1f
File: gsl-ref.info,  Node: Reading and writing combinations,  Next: Combination Examples,  Prev: Combination functions,  Up: Combinations

10.6 Reading and writing combinations
=====================================

The library provides functions for reading and writing combinations to a
file as binary data or formatted text.

 -- Function: int gsl_combination_fwrite (FILE * STREAM, const
          gsl_combination * C)
     This function writes the elements of the combination C to the
     stream STREAM in binary format.  The function returns
     `GSL_EFAILED' if 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.

 -- Function: int gsl_combination_fread (FILE * STREAM, gsl_combination
          * C)
     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_EFAILED' if 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.

 -- Function: int gsl_combination_fprintf (FILE * STREAM, const
          gsl_combination * C, const char * FORMAT)
     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.  In ISO C99
     the type modifier `z' represents `size_t', so `"%zu\n"' is a
     suitable format.(1)  The function returns `GSL_EFAILED' if there
     was a problem writing to the file.

 -- Function: int gsl_combination_fscanf (FILE * STREAM,
          gsl_combination * C)
     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_EFAILED' if there was a problem reading from the file.

   ---------- Footnotes ----------

   (1) In versions of the GNU C library prior to the ISO C99 standard,
the type modifier `Z' was used instead.

\1f
File: gsl-ref.info,  Node: Combination Examples,  Next: Combination References and Further Reading,  Prev: Reading and writing combinations,  Up: Combinations

10.7 Examples
=============

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.

\1f
File: gsl-ref.info,  Node: Combination References and Further Reading,  Prev: Combination Examples,  Up: Combinations

10.8 References and Further Reading
===================================

Further information on combinations can be found in,

     Donald L. Kreher, Douglas R. Stinson, `Combinatorial Algorithms:
     Generation, Enumeration and Search', 1998, CRC Press LLC, ISBN
     084933988X


\1f
File: gsl-ref.info,  Node: Sorting,  Next: BLAS Support,  Prev: Combinations,  Up: Top

11 Sorting
**********

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.

* Menu:

* Sorting objects::
* Sorting vectors::
* Selecting the k smallest or largest elements::
* Computing the rank::
* Sorting Examples::
* Sorting References and Further Reading::

\1f
File: gsl-ref.info,  Node: Sorting objects,  Next: Sorting vectors,  Up: Sorting

11.1 Sorting objects
====================

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'.

 -- Function: void gsl_heapsort (void * ARRAY, size_t COUNT, size_t
          SIZE, gsl_comparison_fn_t COMPARE)
     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, `0' if 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 `qsort' the 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.

 -- Function: int gsl_heapsort_index (size_t * P, const void * ARRAY,
          size_t COUNT, size_t SIZE, gsl_comparison_fn_t COMPARE)
     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.