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 . 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." 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. 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. 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/'). 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. 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::). 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 . 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. 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/'". 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'). 2 Using the library ******************* This chapter describes how to compile programs that use GSL, and introduces its conventions. 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 #include 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. 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 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. 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. 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::. 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 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_'. 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::. 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. 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 ' 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. 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. 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 All types #include double #include long double #include float #include long #include unsigned long #include int #include unsigned int #include short #include unsigned short #include char #include unsigned char 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'. 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. 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. 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. 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. 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'. 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. 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'. 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'. 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); } 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 #include #include ... 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. 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'. 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 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. 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)'. 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 double y = gsl_pow_4 (3.141) /* compute 3.141**4 */ 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). 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. 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. 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. 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'. ---------- Footnotes ---------- (1) Note that the first edition uses different definitions. 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. 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. 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). 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). 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). 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). 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). 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). 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 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'. 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. 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. 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. 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. 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). 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 #include 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). 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. 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'. 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'. 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; 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. 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'. 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. 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. 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). 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). 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'. 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. 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. 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. 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. 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. 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. 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. 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. 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). 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). 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. 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. 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. 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). 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. 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}. 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. 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. 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'. 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. 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. 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. 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. 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)). 7.11 Dilogarithm ================ The functions described in this section are declared in the header file `gsl_sf_dilog.h'. 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). 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. 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. 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. 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. 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) 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. 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. 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. 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. 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'. 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). 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). 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)). 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. 7.16 Exponential Functions ========================== The functions described in this section are declared in the header file `gsl_sf_exp.h'. 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. 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) 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. 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'. 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. 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. 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'). 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. 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. 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. 7.18 Fermi-Dirac Function ========================= The functions described in this section are declared in the header file `gsl_sf_fermi_dirac.h'. 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). 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). 7.19 Gamma and Beta Functions ============================= The functions described in this section are declared in the header file `gsl_sf_gamma.h'. 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