diff options
Diffstat (limited to 'gsl-1.9/doc/gsl-ref.info-1')
| -rw-r--r-- | gsl-1.9/doc/gsl-ref.info-1 | 7392 |
1 files changed, 7392 insertions, 0 deletions
diff --git a/gsl-1.9/doc/gsl-ref.info-1 b/gsl-1.9/doc/gsl-ref.info-1 new file mode 100644 index 0000000..2d4e483 --- /dev/null +++ b/gsl-1.9/doc/gsl-ref.info-1 @@ -0,0 +1,7392 @@ +This is gsl-ref.info, produced by makeinfo version 4.8 from +gsl-ref.texi. + +INFO-DIR-SECTION Scientific software +START-INFO-DIR-ENTRY +* gsl-ref: (gsl-ref). GNU Scientific Library - Reference +END-INFO-DIR-ENTRY + + Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, +2005, 2006, 2007 The GSL Team. + + Permission is granted to copy, distribute and/or modify this document +under the terms of the GNU Free Documentation License, Version 1.2 or +any later version published by the Free Software Foundation; with the +Invariant Sections being "GNU General Public License" and "Free Software +Needs Free Documentation", the Front-Cover text being "A GNU Manual", +and with the Back-Cover Text being (a) (see below). A copy of the +license is included in the section entitled "GNU Free Documentation +License". + + (a) The Back-Cover Text is: "You have freedom to copy and modify this +GNU Manual, like GNU software." + + +File: gsl-ref.info, Node: Top, Next: Introduction, Prev: (dir), Up: (dir) + +GSL +*** + +This file documents the GNU Scientific Library (GSL), a collection of +numerical routines for scientific computing. It corresponds to release +1.9 of the library. Please report any errors in this manual to +<bug-gsl@gnu.org>. + + More information about GSL can be found at the project homepage, +`http://www.gnu.org/software/gsl/'. + + Printed copies of this manual can be purchased from Network Theory +Ltd at `http://www.network-theory.co.uk/gsl/manual/'. The money raised +from sales of the manual helps support the development of GSL. + + A Japanese translation of this manual is available from the GSL +project homepage thanks to Daisuke Tominaga. + + Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, +2005, 2006, 2007 The GSL Team. + + Permission is granted to copy, distribute and/or modify this document +under the terms of the GNU Free Documentation License, Version 1.2 or +any later version published by the Free Software Foundation; with the +Invariant Sections being "GNU General Public License" and "Free Software +Needs Free Documentation", the Front-Cover text being "A GNU Manual", +and with the Back-Cover Text being (a) (see below). A copy of the +license is included in the section entitled "GNU Free Documentation +License". + + (a) The Back-Cover Text is: "You have freedom to copy and modify this +GNU Manual, like GNU software." + + +* Menu: + +* Introduction:: +* Using the library:: +* Error Handling:: +* Mathematical Functions:: +* Complex Numbers:: +* Polynomials:: +* Special Functions:: +* Vectors and Matrices:: +* Permutations:: +* Combinations:: +* Sorting:: +* BLAS Support:: +* Linear Algebra:: +* Eigensystems:: +* Fast Fourier Transforms:: +* Numerical Integration:: +* Random Number Generation:: +* Quasi-Random Sequences:: +* Random Number Distributions:: +* Statistics:: +* Histograms:: +* N-tuples:: +* Monte Carlo Integration:: +* Simulated Annealing:: +* Ordinary Differential Equations:: +* Interpolation:: +* Numerical Differentiation:: +* Chebyshev Approximations:: +* Series Acceleration:: +* Wavelet Transforms:: +* Discrete Hankel Transforms:: +* One dimensional Root-Finding:: +* One dimensional Minimization:: +* Multidimensional Root-Finding:: +* Multidimensional Minimization:: +* Least-Squares Fitting:: +* Nonlinear Least-Squares Fitting:: +* Basis Splines:: +* Physical Constants:: +* IEEE floating-point arithmetic:: +* Debugging Numerical Programs:: +* Contributors to GSL:: +* Autoconf Macros:: +* GSL CBLAS Library:: +* Free Software Needs Free Documentation:: +* GNU General Public License:: +* GNU Free Documentation License:: +* Function Index:: +* Variable Index:: +* Type Index:: +* Concept Index:: + + +File: gsl-ref.info, Node: Introduction, Next: Using the library, Prev: Top, Up: Top + +1 Introduction +************** + +The GNU Scientific Library (GSL) is a collection of routines for +numerical computing. The routines have been written from scratch in C, +and present a modern Applications Programming Interface (API) for C +programmers, allowing wrappers to be written for very high level +languages. The source code is distributed under the GNU General Public +License. + +* Menu: + +* Routines available in GSL:: +* GSL is Free Software:: +* Obtaining GSL:: +* No Warranty:: +* Reporting Bugs:: +* Further Information:: +* Conventions used in this manual:: + + +File: gsl-ref.info, Node: Routines available in GSL, Next: GSL is Free Software, Up: Introduction + +1.1 Routines available in GSL +============================= + +The library covers a wide range of topics in numerical computing. +Routines are available for the following areas, + + Complex Numbers Roots of Polynomials + Special Functions Vectors and Matrices + Permutations Combinations + Sorting BLAS Support + Linear Algebra CBLAS Library + Fast Fourier Transforms Eigensystems + Random Numbers Quadrature + Random Distributions Quasi-Random Sequences + Histograms Statistics + Monte Carlo Integration N-Tuples + Differential Equations Simulated Annealing + Numerical Differentiation Interpolation + Series Acceleration Chebyshev Approximations + Root-Finding Discrete Hankel Transforms + Least-Squares Fitting Minimization + IEEE Floating-Point Physical Constants + 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. + + +File: gsl-ref.info, Node: GSL is Free Software, Next: Obtaining GSL, Prev: Routines available in GSL, Up: Introduction + +1.2 GSL is Free Software +======================== + +The subroutines in the GNU Scientific Library are "free software"; this +means that everyone is free to use them, and to redistribute them in +other free programs. The library is not in the public domain; it is +copyrighted and there are conditions on its distribution. These +conditions are designed to permit everything that a good cooperating +citizen would want to do. What is not allowed is to try to prevent +others from further sharing any version of the software that they might +get from you. + + Specifically, we want to make sure that you have the right to share +copies of programs that you are given which use the GNU Scientific +Library, that you receive their source code or else can get it if you +want it, that you can change these programs or use pieces of them in new +free programs, and that you know you can do these things. + + To make sure that everyone has such rights, we have to forbid you to +deprive anyone else of these rights. For example, if you distribute +copies of any code which uses the GNU Scientific Library, you must give +the recipients all the rights that you have received. You must make +sure that they, too, receive or can get the source code, both to the +library and the code which uses it. And you must tell them their +rights. This means that the library should not be redistributed in +proprietary programs. + + Also, for our own protection, we must make certain that everyone +finds out that there is no warranty for the GNU Scientific Library. If +these programs are modified by someone else and passed on, we want their +recipients to know that what they have is not what we distributed, so +that any problems introduced by others will not reflect on our +reputation. + + The precise conditions for the distribution of software related to +the GNU Scientific Library are found in the GNU General Public License +(*note GNU General Public License::). Further information about this +license is available from the GNU Project webpage `Frequently Asked +Questions about the GNU GPL', + + `http://www.gnu.org/copyleft/gpl-faq.html' + +The Free Software Foundation also operates a license consulting service +for commercial users (contact details available from +`http://www.fsf.org/'). + + +File: gsl-ref.info, Node: Obtaining GSL, Next: No Warranty, Prev: GSL is Free Software, Up: Introduction + +1.3 Obtaining GSL +================= + +The source code for the library can be obtained in different ways, by +copying it from a friend, purchasing it on CDROM or downloading it from +the internet. A list of public ftp servers which carry the source code +can be found on the GNU website, + + `http://www.gnu.org/software/gsl/' + +The preferred platform for the library is a GNU system, which allows it +to take advantage of additional features in the GNU C compiler and GNU C +library. However, the library is fully portable and should compile on +most systems with a C compiler. Precompiled versions of the library +can be purchased from commercial redistributors listed on the website +above. + + Announcements of new releases, updates and other relevant events are +made on the `info-gsl@gnu.org' mailing list. To subscribe to this +low-volume list, send an email of the following form: + + To: info-gsl-request@gnu.org + Subject: subscribe + +You will receive a response asking you to reply in order to confirm +your subscription. + + +File: gsl-ref.info, Node: No Warranty, Next: Reporting Bugs, Prev: Obtaining GSL, Up: Introduction + +1.4 No Warranty +=============== + +The software described in this manual has no warranty, it is provided +"as is". It is your responsibility to validate the behavior of the +routines and their accuracy using the source code provided, or to +purchase support and warranties from commercial redistributors. Consult +the GNU General Public license for further details (*note GNU General +Public License::). + + +File: gsl-ref.info, Node: Reporting Bugs, Next: Further Information, Prev: No Warranty, Up: Introduction + +1.5 Reporting Bugs +================== + +A list of known bugs can be found in the `BUGS' file included in the +GSL distribution. Details of compilation problems can be found in the +`INSTALL' file. + + If you find a bug which is not listed in these files, please report +it to <bug-gsl@gnu.org>. + + All bug reports should include: + + * The version number of GSL + + * The hardware and operating system + + * The compiler used, including version number and compilation options + + * A description of the bug behavior + + * A short program which exercises the bug + +It is useful if you can check whether the same problem occurs when the +library is compiled without optimization. Thank you. + + Any errors or omissions in this manual can also be reported to the +same address. + + +File: gsl-ref.info, Node: Further Information, Next: Conventions used in this manual, Prev: Reporting Bugs, Up: Introduction + +1.6 Further Information +======================= + +Additional information, including online copies of this manual, links to +related projects, and mailing list archives are available from the +website mentioned above. + + Any questions about the use and installation of the library can be +asked on the mailing list `help-gsl@gnu.org'. To subscribe to this +list, send an email of the following form: + + To: help-gsl-request@gnu.org + Subject: subscribe + +This mailing list can be used to ask questions not covered by this +manual, and to contact the developers of the library. + + If you would like to refer to the GNU Scientific Library in a journal +article, the recommended way is to cite this reference manual, e.g. `M. +Galassi et al, GNU Scientific Library Reference Manual (2nd Ed.), ISBN +0954161734'. + + If you want to give a url, use "`http://www.gnu.org/software/gsl/'". + + +File: gsl-ref.info, Node: Conventions used in this manual, Prev: Further Information, Up: Introduction + +1.7 Conventions used in this manual +=================================== + +This manual contains many examples which can be typed at the keyboard. +A command entered at the terminal is shown like this, + + $ command + +The first character on the line is the terminal prompt, and should not +be typed. The dollar sign `$' is used as the standard prompt in this +manual, although some systems may use a different character. + + The examples assume the use of the GNU operating system. There may +be minor differences in the output on other systems. The commands for +setting environment variables use the Bourne shell syntax of the +standard GNU shell (`bash'). + + +File: gsl-ref.info, Node: Using the library, Next: Error Handling, Prev: Introduction, Up: Top + +2 Using the library +******************* + +This chapter describes how to compile programs that use GSL, and +introduces its conventions. + +* Menu: + +* An Example Program:: +* Compiling and Linking:: +* Shared Libraries:: +* ANSI C Compliance:: +* Inline functions:: +* Long double:: +* Portability functions:: +* Alternative optimized functions:: +* Support for different numeric types:: +* Compatibility with C++:: +* Aliasing of arrays:: +* Thread-safety:: +* Deprecated Functions:: +* Code Reuse:: + + +File: gsl-ref.info, Node: An Example Program, Next: Compiling and Linking, Up: Using the library + +2.1 An Example Program +====================== + +The following short program demonstrates the use of the library by +computing the value of the Bessel function J_0(x) for x=5, + + #include <stdio.h> + #include <gsl/gsl_sf_bessel.h> + + int + main (void) + { + double x = 5.0; + double y = gsl_sf_bessel_J0 (x); + printf ("J0(%g) = %.18e\n", x, y); + return 0; + } + +The output is shown below, and should be correct to double-precision +accuracy, + + J0(5) = -1.775967713143382920e-01 + +The steps needed to compile this program are described in the following +sections. + + +File: gsl-ref.info, Node: Compiling and Linking, Next: Shared Libraries, Prev: An Example Program, Up: Using the library + +2.2 Compiling and Linking +========================= + +The library header files are installed in their own `gsl' directory. +You should write any preprocessor include statements with a `gsl/' +directory prefix thus, + + #include <gsl/gsl_math.h> + +If the directory is not installed on the standard search path of your +compiler you will also need to provide its location to the preprocessor +as a command line flag. The default location of the `gsl' directory is +`/usr/local/include/gsl'. A typical compilation command for a source +file `example.c' with the GNU C compiler `gcc' is, + + $ gcc -Wall -I/usr/local/include -c example.c + +This results in an object file `example.o'. The default include path +for `gcc' searches `/usr/local/include' automatically so the `-I' +option can actually be omitted when GSL is installed in its default +location. + +* Menu: + +* Linking programs with the library:: +* Linking with an alternative BLAS library:: + + +File: gsl-ref.info, Node: Linking programs with the library, Next: Linking with an alternative BLAS library, Up: Compiling and Linking + +2.2.1 Linking programs with the library +--------------------------------------- + +The library is installed as a single file, `libgsl.a'. A shared +version of the library `libgsl.so' is also installed on systems that +support shared libraries. The default location of these files is +`/usr/local/lib'. If this directory is not on the standard search path +of your linker you will also need to provide its location as a command +line flag. + + To link against the library you need to specify both the main +library and a supporting CBLAS library, which provides standard basic +linear algebra subroutines. A suitable CBLAS implementation is +provided in the library `libgslcblas.a' if your system does not provide +one. The following example shows how to link an application with the +library, + + $ gcc -L/usr/local/lib example.o -lgsl -lgslcblas -lm + +The default library path for `gcc' searches `/usr/local/lib' +automatically so the `-L' option can be omitted when GSL is installed +in its default location. + + +File: gsl-ref.info, Node: Linking with an alternative BLAS library, Prev: Linking programs with the library, Up: Compiling and Linking + +2.2.2 Linking with an alternative BLAS library +---------------------------------------------- + +The following command line shows how you would link the same application +with an alternative CBLAS library called `libcblas', + + $ gcc example.o -lgsl -lcblas -lm + +For the best performance an optimized platform-specific CBLAS library +should be used for `-lcblas'. The library must conform to the CBLAS +standard. The ATLAS package provides a portable high-performance BLAS +library with a CBLAS interface. It is free software and should be +installed for any work requiring fast vector and matrix operations. +The following command line will link with the ATLAS library and its +CBLAS interface, + + $ gcc example.o -lgsl -lcblas -latlas -lm + +For more information see *Note BLAS Support::. + + +File: gsl-ref.info, Node: Shared Libraries, Next: ANSI C Compliance, Prev: Compiling and Linking, Up: Using the library + +2.3 Shared Libraries +==================== + +To run a program linked with the shared version of the library the +operating system must be able to locate the corresponding `.so' file at +runtime. If the library cannot be found, the following error will +occur: + + $ ./a.out + ./a.out: error while loading shared libraries: + libgsl.so.0: cannot open shared object file: No such + file or directory + +To avoid this error, define the shell variable `LD_LIBRARY_PATH' to +include the directory where the library is installed. + + For example, in the Bourne shell (`/bin/sh' or `/bin/bash'), the +library search path can be set with the following commands: + + $ LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH + $ export LD_LIBRARY_PATH + $ ./example + +In the C-shell (`/bin/csh' or `/bin/tcsh') the equivalent command is, + + % setenv LD_LIBRARY_PATH /usr/local/lib:$LD_LIBRARY_PATH + +The standard prompt for the C-shell in the example above is the percent +character `%', and should not be typed as part of the command. + + To save retyping these commands each session they should be placed +in an individual or system-wide login file. + + To compile a statically linked version of the program, use the +`-static' flag in `gcc', + + $ gcc -static example.o -lgsl -lgslcblas -lm + + +File: gsl-ref.info, Node: ANSI C Compliance, Next: Inline functions, Prev: Shared Libraries, Up: Using the library + +2.4 ANSI C Compliance +===================== + +The library is written in ANSI C and is intended to conform to the ANSI +C standard (C89). It should be portable to any system with a working +ANSI C compiler. + + The library does not rely on any non-ANSI extensions in the +interface it exports to the user. Programs you write using GSL can be +ANSI compliant. Extensions which can be used in a way compatible with +pure ANSI C are supported, however, via conditional compilation. This +allows the library to take advantage of compiler extensions on those +platforms which support them. + + When an ANSI C feature is known to be broken on a particular system +the library will exclude any related functions at compile-time. This +should make it impossible to link a program that would use these +functions and give incorrect results. + + To avoid namespace conflicts all exported function names and +variables have the prefix `gsl_', while exported macros have the prefix +`GSL_'. + + +File: gsl-ref.info, Node: Inline functions, Next: Long double, Prev: ANSI C Compliance, Up: Using the library + +2.5 Inline functions +==================== + +The `inline' keyword is not part of the original ANSI C standard (C89) +and the library does not export any inline function definitions by +default. However, the library provides optional inline versions of +performance-critical functions by conditional compilation. The inline +versions of these functions can be included by defining the macro +`HAVE_INLINE' when compiling an application, + + $ gcc -Wall -c -DHAVE_INLINE example.c + +If you use `autoconf' this macro can be defined automatically. If you +do not define the macro `HAVE_INLINE' then the slower non-inlined +versions of the functions will be used instead. + + Note that the actual usage of the inline keyword is `extern inline', +which eliminates unnecessary function definitions in GCC. If the form +`extern inline' causes problems with other compilers a stricter +autoconf test can be used, see *Note Autoconf Macros::. + + +File: gsl-ref.info, Node: Long double, Next: Portability functions, Prev: Inline functions, Up: Using the library + +2.6 Long double +=============== + +The extended numerical type `long double' is part of the ANSI C +standard and should be available in every modern compiler. However, the +precision of `long double' is platform dependent, and this should be +considered when using it. The IEEE standard only specifies the minimum +precision of extended precision numbers, while the precision of +`double' is the same on all platforms. + + In some system libraries the `stdio.h' formatted input/output +functions `printf' and `scanf' are not implemented correctly for `long +double'. Undefined or incorrect results are avoided by testing these +functions during the `configure' stage of library compilation and +eliminating certain GSL functions which depend on them if necessary. +The corresponding line in the `configure' output looks like this, + + checking whether printf works with long double... no + +Consequently when `long double' formatted input/output does not work on +a given system it should be impossible to link a program which uses GSL +functions dependent on this. + + If it is necessary to work on a system which does not support +formatted `long double' input/output then the options are to use binary +formats or to convert `long double' results into `double' for reading +and writing. + + +File: gsl-ref.info, Node: Portability functions, Next: Alternative optimized functions, Prev: Long double, Up: Using the library + +2.7 Portability functions +========================= + +To help in writing portable applications GSL provides some +implementations of functions that are found in other libraries, such as +the BSD math library. You can write your application to use the native +versions of these functions, and substitute the GSL versions via a +preprocessor macro if they are unavailable on another platform. + + For example, after determining whether the BSD function `hypot' is +available you can include the following macro definitions in a file +`config.h' with your application, + + /* Substitute gsl_hypot for missing system hypot */ + + #ifndef HAVE_HYPOT + #define hypot gsl_hypot + #endif + +The application source files can then use the include command `#include +<config.h>' to replace each occurrence of `hypot' by `gsl_hypot' when +`hypot' is not available. This substitution can be made automatically +if you use `autoconf', see *Note Autoconf Macros::. + + In most circumstances the best strategy is to use the native +versions of these functions when available, and fall back to GSL +versions otherwise, since this allows your application to take +advantage of any platform-specific optimizations in the system library. +This is the strategy used within GSL itself. + + +File: gsl-ref.info, Node: Alternative optimized functions, Next: Support for different numeric types, Prev: Portability functions, Up: Using the library + +2.8 Alternative optimized functions +=================================== + +The main implementation of some functions in the library will not be +optimal on all architectures. For example, there are several ways to +compute a Gaussian random variate and their relative speeds are +platform-dependent. In cases like this the library provides alternative +implementations of these functions with the same interface. If you +write your application using calls to the standard implementation you +can select an alternative version later via a preprocessor definition. +It is also possible to introduce your own optimized functions this way +while retaining portability. The following lines demonstrate the use of +a platform-dependent choice of methods for sampling from the Gaussian +distribution, + + #ifdef SPARC + #define gsl_ran_gaussian gsl_ran_gaussian_ratio_method + #endif + #ifdef INTEL + #define gsl_ran_gaussian my_gaussian + #endif + +These lines would be placed in the configuration header file `config.h' +of the application, which should then be included by all the source +files. Note that the alternative implementations will not produce +bit-for-bit identical results, and in the case of random number +distributions will produce an entirely different stream of random +variates. + + +File: gsl-ref.info, Node: Support for different numeric types, Next: Compatibility with C++, Prev: Alternative optimized functions, Up: Using the library + +2.9 Support for different numeric types +======================================= + +Many functions in the library are defined for different numeric types. +This feature is implemented by varying the name of the function with a +type-related modifier--a primitive form of C++ templates. The modifier +is inserted into the function name after the initial module prefix. +The following table shows the function names defined for all the +numeric types of an imaginary module `gsl_foo' with function `fn', + + gsl_foo_fn double + gsl_foo_long_double_fn long double + gsl_foo_float_fn float + gsl_foo_long_fn long + gsl_foo_ulong_fn unsigned long + gsl_foo_int_fn int + gsl_foo_uint_fn unsigned int + gsl_foo_short_fn short + gsl_foo_ushort_fn unsigned short + gsl_foo_char_fn char + gsl_foo_uchar_fn unsigned char + +The normal numeric precision `double' is considered the default and +does not require a suffix. For example, the function `gsl_stats_mean' +computes the mean of double precision numbers, while the function +`gsl_stats_int_mean' computes the mean of integers. + + A corresponding scheme is used for library defined types, such as +`gsl_vector' and `gsl_matrix'. In this case the modifier is appended +to the type name. For example, if a module defines a new +type-dependent struct or typedef `gsl_foo' it is modified for other +types in the following way, + + gsl_foo double + gsl_foo_long_double long double + gsl_foo_float float + gsl_foo_long long + gsl_foo_ulong unsigned long + gsl_foo_int int + gsl_foo_uint unsigned int + gsl_foo_short short + gsl_foo_ushort unsigned short + gsl_foo_char char + gsl_foo_uchar unsigned char + +When a module contains type-dependent definitions the library provides +individual header files for each type. The filenames are modified as +shown in the below. For convenience the default header includes the +definitions for all the types. To include only the double precision +header file, or any other specific type, use its individual filename. + + #include <gsl/gsl_foo.h> All types + #include <gsl/gsl_foo_double.h> double + #include <gsl/gsl_foo_long_double.h> long double + #include <gsl/gsl_foo_float.h> float + #include <gsl/gsl_foo_long.h> long + #include <gsl/gsl_foo_ulong.h> unsigned long + #include <gsl/gsl_foo_int.h> int + #include <gsl/gsl_foo_uint.h> unsigned int + #include <gsl/gsl_foo_short.h> short + #include <gsl/gsl_foo_ushort.h> unsigned short + #include <gsl/gsl_foo_char.h> char + #include <gsl/gsl_foo_uchar.h> unsigned char + + +File: gsl-ref.info, Node: Compatibility with C++, Next: Aliasing of arrays, Prev: Support for different numeric types, Up: Using the library + +2.10 Compatibility with C++ +=========================== + +The library header files automatically define functions to have `extern +"C"' linkage when included in C++ programs. This allows the functions +to be called directly from C++. + + To use C++ exception handling within user-defined functions passed to +the library as parameters, the library must be built with the +additional `CFLAGS' compilation option `-fexceptions'. + + +File: gsl-ref.info, Node: Aliasing of arrays, Next: Thread-safety, Prev: Compatibility with C++, Up: Using the library + +2.11 Aliasing of arrays +======================= + +The library assumes that arrays, vectors and matrices passed as +modifiable arguments are not aliased and do not overlap with each other. +This removes the need for the library to handle overlapping memory +regions as a special case, and allows additional optimizations to be +used. If overlapping memory regions are passed as modifiable arguments +then the results of such functions will be undefined. If the arguments +will not be modified (for example, if a function prototype declares them +as `const' arguments) then overlapping or aliased memory regions can be +safely used. + + +File: gsl-ref.info, Node: Thread-safety, Next: Deprecated Functions, Prev: Aliasing of arrays, Up: Using the library + +2.12 Thread-safety +================== + +The library can be used in multi-threaded programs. All the functions +are thread-safe, in the sense that they do not use static variables. +Memory is always associated with objects and not with functions. For +functions which use "workspace" objects as temporary storage the +workspaces should be allocated on a per-thread basis. For functions +which use "table" objects as read-only memory the tables can be used by +multiple threads simultaneously. Table arguments are always declared +`const' in function prototypes, to indicate that they may be safely +accessed by different threads. + + There are a small number of static global variables which are used to +control the overall behavior of the library (e.g. whether to use +range-checking, the function to call on fatal error, etc). These +variables are set directly by the user, so they should be initialized +once at program startup and not modified by different threads. + + +File: gsl-ref.info, Node: Deprecated Functions, Next: Code Reuse, Prev: Thread-safety, Up: Using the library + +2.13 Deprecated Functions +========================= + +From time to time, it may be necessary for the definitions of some +functions to be altered or removed from the library. In these +circumstances the functions will first be declared "deprecated" and +then removed from subsequent versions of the library. Functions that +are deprecated can be disabled in the current release by setting the +preprocessor definition `GSL_DISABLE_DEPRECATED'. This allows existing +code to be tested for forwards compatibility. + + +File: gsl-ref.info, Node: Code Reuse, Prev: Deprecated Functions, Up: Using the library + +2.14 Code Reuse +=============== + +Where possible the routines in the library have been written to avoid +dependencies between modules and files. This should make it possible to +extract individual functions for use in your own applications, without +needing to have the whole library installed. You may need to define +certain macros such as `GSL_ERROR' and remove some `#include' +statements in order to compile the files as standalone units. Reuse of +the library code in this way is encouraged, subject to the terms of the +GNU General Public License. + + +File: gsl-ref.info, Node: Error Handling, Next: Mathematical Functions, Prev: Using the library, Up: Top + +3 Error Handling +**************** + +This chapter describes the way that GSL functions report and handle +errors. By examining the status information returned by every function +you can determine whether it succeeded or failed, and if it failed you +can find out what the precise cause of failure was. You can also define +your own error handling functions to modify the default behavior of the +library. + + The functions described in this section are declared in the header +file `gsl_errno.h'. + +* Menu: + +* Error Reporting:: +* Error Codes:: +* Error Handlers:: +* Using GSL error reporting in your own functions:: +* Error Reporting Examples:: + + +File: gsl-ref.info, Node: Error Reporting, Next: Error Codes, Up: Error Handling + +3.1 Error Reporting +=================== + +The library follows the thread-safe error reporting conventions of the +POSIX Threads library. Functions return a non-zero error code to +indicate an error and `0' to indicate success. + + int status = gsl_function (...) + + if (status) { /* an error occurred */ + ..... + /* status value specifies the type of error */ + } + + The routines report an error whenever they cannot perform the task +requested of them. For example, a root-finding function would return a +non-zero error code if could not converge to the requested accuracy, or +exceeded a limit on the number of iterations. Situations like this are +a normal occurrence when using any mathematical library and you should +check the return status of the functions that you call. + + Whenever a routine reports an error the return value specifies the +type of error. The return value is analogous to the value of the +variable `errno' in the C library. The caller can examine the return +code and decide what action to take, including ignoring the error if it +is not considered serious. + + In addition to reporting errors by return codes the library also has +an error handler function `gsl_error'. This function is called by +other library functions when they report an error, just before they +return to the caller. The default behavior of the error handler is to +print a message and abort the program, + + gsl: file.c:67: ERROR: invalid argument supplied by user + Default GSL error handler invoked. + Aborted + + The purpose of the `gsl_error' handler is to provide a function +where a breakpoint can be set that will catch library errors when +running under the debugger. It is not intended for use in production +programs, which should handle any errors using the return codes. + + +File: gsl-ref.info, Node: Error Codes, Next: Error Handlers, Prev: Error Reporting, Up: Error Handling + +3.2 Error Codes +=============== + +The error code numbers returned by library functions are defined in the +file `gsl_errno.h'. They all have the prefix `GSL_' and expand to +non-zero constant integer values. Error codes above 1024 are reserved +for applications, and are not used by the library. Many of the error +codes use the same base name as the corresponding error code in the C +library. Here are some of the most common error codes, + + -- Macro: int GSL_EDOM + Domain error; used by mathematical functions when an argument + value does not fall into the domain over which the function is + defined (like EDOM in the C library) + + -- Macro: int GSL_ERANGE + Range error; used by mathematical functions when the result value + is not representable because of overflow or underflow (like ERANGE + in the C library) + + -- Macro: int GSL_ENOMEM + No memory available. The system cannot allocate more virtual + memory because its capacity is full (like ENOMEM in the C + library). This error is reported when a GSL routine encounters + problems when trying to allocate memory with `malloc'. + + -- Macro: int GSL_EINVAL + Invalid argument. This is used to indicate various kinds of + problems with passing the wrong argument to a library function + (like EINVAL in the C library). + + The error codes can be converted into an error message using the +function `gsl_strerror'. + + -- Function: const char * gsl_strerror (const int GSL_ERRNO) + This function returns a pointer to a string describing the error + code GSL_ERRNO. For example, + + printf ("error: %s\n", gsl_strerror (status)); + + would print an error message like `error: output range error' for a + status value of `GSL_ERANGE'. + + +File: gsl-ref.info, Node: Error Handlers, Next: Using GSL error reporting in your own functions, Prev: Error Codes, Up: Error Handling + +3.3 Error Handlers +================== + +The default behavior of the GSL error handler is to print a short +message and call `abort'. When this default is in use programs will +stop with a core-dump whenever a library routine reports an error. +This is intended as a fail-safe default for programs which do not check +the return status of library routines (we don't encourage you to write +programs this way). + + If you turn off the default error handler it is your responsibility +to check the return values of routines and handle them yourself. You +can also customize the error behavior by providing a new error handler. +For example, an alternative error handler could log all errors to a +file, ignore certain error conditions (such as underflows), or start the +debugger and attach it to the current process when an error occurs. + + All GSL error handlers have the type `gsl_error_handler_t', which is +defined in `gsl_errno.h', + + -- Data Type: gsl_error_handler_t + This is the type of GSL error handler functions. An error handler + will be passed four arguments which specify the reason for the + error (a string), the name of the source file in which it occurred + (also a string), the line number in that file (an integer) and the + error number (an integer). The source file and line number are + set at compile time using the `__FILE__' and `__LINE__' directives + in the preprocessor. An error handler function returns type + `void'. Error handler functions should be defined like this, + + void handler (const char * reason, + const char * file, + int line, + int gsl_errno) + +To request the use of your own error handler you need to call the +function `gsl_set_error_handler' which is also declared in +`gsl_errno.h', + + -- Function: gsl_error_handler_t * gsl_set_error_handler + (gsl_error_handler_t NEW_HANDLER) + This function sets a new error handler, NEW_HANDLER, for the GSL + library routines. The previous handler is returned (so that you + can restore it later). Note that the pointer to a user defined + error handler function is stored in a static variable, so there + can be only one error handler per program. This function should + be not be used in multi-threaded programs except to set up a + program-wide error handler from a master thread. The following + example shows how to set and restore a new error handler, + + /* save original handler, install new handler */ + old_handler = gsl_set_error_handler (&my_handler); + + /* code uses new handler */ + ..... + + /* restore original handler */ + gsl_set_error_handler (old_handler); + + To use the default behavior (`abort' on error) set the error + handler to `NULL', + + old_handler = gsl_set_error_handler (NULL); + + -- Function: gsl_error_handler_t * gsl_set_error_handler_off () + This function turns off the error handler by defining an error + handler which does nothing. This will cause the program to + continue after any error, so the return values from any library + routines must be checked. This is the recommended behavior for + production programs. The previous handler is returned (so that + you can restore it later). + + The error behavior can be changed for specific applications by +recompiling the library with a customized definition of the `GSL_ERROR' +macro in the file `gsl_errno.h'. + + +File: gsl-ref.info, Node: Using GSL error reporting in your own functions, Next: Error Reporting Examples, Prev: Error Handlers, Up: Error Handling + +3.4 Using GSL error reporting in your own functions +=================================================== + +If you are writing numerical functions in a program which also uses GSL +code you may find it convenient to adopt the same error reporting +conventions as in the library. + + To report an error you need to call the function `gsl_error' with a +string describing the error and then return an appropriate error code +from `gsl_errno.h', or a special value, such as `NaN'. For convenience +the file `gsl_errno.h' defines two macros which carry out these steps: + + -- Macro: GSL_ERROR (REASON, GSL_ERRNO) + This macro reports an error using the GSL conventions and returns a + status value of `gsl_errno'. It expands to the following code + fragment, + + gsl_error (reason, __FILE__, __LINE__, gsl_errno); + return gsl_errno; + + The macro definition in `gsl_errno.h' actually wraps the code in a + `do { ... } while (0)' block to prevent possible parsing problems. + + Here is an example of how the macro could be used to report that a +routine did not achieve a requested tolerance. To report the error the +routine needs to return the error code `GSL_ETOL'. + + if (residual > tolerance) + { + GSL_ERROR("residual exceeds tolerance", GSL_ETOL); + } + + -- Macro: GSL_ERROR_VAL (REASON, GSL_ERRNO, VALUE) + This macro is the same as `GSL_ERROR' but returns a user-defined + value of VALUE instead of an error code. It can be used for + mathematical functions that return a floating point value. + + The following example shows how to return a `NaN' at a mathematical +singularity using the `GSL_ERROR_VAL' macro, + + if (x == 0) + { + GSL_ERROR_VAL("argument lies on singularity", + GSL_ERANGE, GSL_NAN); + } + + +File: gsl-ref.info, Node: Error Reporting Examples, Prev: Using GSL error reporting in your own functions, Up: Error Handling + +3.5 Examples +============ + +Here is an example of some code which checks the return value of a +function where an error might be reported, + + #include <stdio.h> + #include <gsl/gsl_errno.h> + #include <gsl/gsl_fft_complex.h> + + ... + int status; + size_t n = 37; + + gsl_set_error_handler_off(); + + status = gsl_fft_complex_radix2_forward (data, 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. + + +File: gsl-ref.info, Node: Mathematical Functions, Next: Complex Numbers, Prev: Error Handling, Up: Top + +4 Mathematical Functions +************************ + +This chapter describes basic mathematical functions. Some of these +functions are present in system libraries, but the alternative versions +given here can be used as a substitute when the system functions are not +available. + + The functions and macros described in this chapter are defined in the +header file `gsl_math.h'. + +* Menu: + +* Mathematical Constants:: +* Infinities and Not-a-number:: +* Elementary Functions:: +* Small integer powers:: +* Testing the Sign of Numbers:: +* Testing for Odd and Even Numbers:: +* Maximum and Minimum functions:: +* Approximate Comparison of Floating Point Numbers:: + + +File: gsl-ref.info, Node: Mathematical Constants, Next: Infinities and Not-a-number, Up: Mathematical Functions + +4.1 Mathematical Constants +========================== + +The library ensures that the standard BSD mathematical constants are +defined. For reference, here is a list of the constants: + +`M_E' + The base of exponentials, e + +`M_LOG2E' + The base-2 logarithm of e, \log_2 (e) + +`M_LOG10E' + The base-10 logarithm of e, \log_10 (e) + +`M_SQRT2' + The square root of two, \sqrt 2 + +`M_SQRT1_2' + The square root of one-half, \sqrt{1/2} + +`M_SQRT3' + The square root of three, \sqrt 3 + +`M_PI' + The constant pi, \pi + +`M_PI_2' + Pi divided by two, \pi/2 + +`M_PI_4' + Pi divided by four, \pi/4 + +`M_SQRTPI' + The square root of pi, \sqrt\pi + +`M_2_SQRTPI' + Two divided by the square root of pi, 2/\sqrt\pi + +`M_1_PI' + The reciprocal of pi, 1/\pi + +`M_2_PI' + Twice the reciprocal of pi, 2/\pi + +`M_LN10' + The natural logarithm of ten, \ln(10) + +`M_LN2' + The natural logarithm of two, \ln(2) + +`M_LNPI' + The natural logarithm of pi, \ln(\pi) + +`M_EULER' + Euler's constant, \gamma + + + +File: gsl-ref.info, Node: Infinities and Not-a-number, Next: Elementary Functions, Prev: Mathematical Constants, Up: Mathematical Functions + +4.2 Infinities and Not-a-number +=============================== + + -- Macro: GSL_POSINF + This macro contains the IEEE representation of positive infinity, + +\infty. It is computed from the expression `+1.0/0.0'. + + -- Macro: GSL_NEGINF + This macro contains the IEEE representation of negative infinity, + -\infty. It is computed from the expression `-1.0/0.0'. + + -- Macro: GSL_NAN + This macro contains the IEEE representation of the Not-a-Number + symbol, `NaN'. It is computed from the ratio `0.0/0.0'. + + -- Function: int gsl_isnan (const double X) + This function returns 1 if X is not-a-number. + + -- Function: int gsl_isinf (const double X) + This function returns +1 if X is positive infinity, -1 if X is + negative infinity and 0 otherwise. + + -- Function: int gsl_finite (const double X) + This function returns 1 if X is a real number, and 0 if it is + infinite or not-a-number. + + +File: gsl-ref.info, Node: Elementary Functions, Next: Small integer powers, Prev: Infinities and Not-a-number, Up: Mathematical Functions + +4.3 Elementary Functions +======================== + +The following routines provide portable implementations of functions +found in the BSD math library. When native versions are not available +the functions described here can be used instead. The substitution can +be made automatically if you use `autoconf' to compile your application +(*note Portability functions::). + + -- Function: double gsl_log1p (const double X) + This function computes the value of \log(1+x) in a way that is + accurate for small X. It provides an alternative to the BSD math + function `log1p(x)'. + + -- Function: double gsl_expm1 (const double X) + This function computes the value of \exp(x)-1 in a way that is + accurate for small X. It provides an alternative to the BSD math + function `expm1(x)'. + + -- Function: double gsl_hypot (const double X, const double Y) + This function computes the value of \sqrt{x^2 + y^2} in a way that + avoids overflow. It provides an alternative to the BSD math + function `hypot(x,y)'. + + -- Function: double gsl_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)'. + + +File: gsl-ref.info, Node: Small integer powers, Next: Testing the Sign of Numbers, Prev: Elementary Functions, Up: Mathematical Functions + +4.4 Small integer powers +======================== + +A common complaint about the standard C library is its lack of a +function for calculating (small) integer powers. GSL provides some +simple functions to fill this gap. For reasons of efficiency, these +functions do not check for overflow or underflow conditions. + + -- Function: double gsl_pow_int (double X, int N) + This routine computes the power x^n for integer N. The power is + computed efficiently--for example, x^8 is computed as ((x^2)^2)^2, + requiring only 3 multiplications. A version of this function + which also computes the numerical error in the result is available + as `gsl_sf_pow_int_e'. + + -- Function: double gsl_pow_2 (const double X) + -- Function: double gsl_pow_3 (const double X) + -- Function: double gsl_pow_4 (const double X) + -- Function: double gsl_pow_5 (const double X) + -- Function: double gsl_pow_6 (const double X) + -- Function: double gsl_pow_7 (const double X) + -- Function: double gsl_pow_8 (const double X) + -- Function: double gsl_pow_9 (const double X) + These functions can be used to compute small integer powers x^2, + x^3, etc. efficiently. The functions will be inlined when possible + so that use of these functions should be as efficient as + explicitly writing the corresponding product expression. + + #include <gsl/gsl_math.h> + double y = gsl_pow_4 (3.141) /* compute 3.141**4 */ + + +File: gsl-ref.info, Node: Testing the Sign of Numbers, Next: Testing for Odd and Even Numbers, Prev: Small integer powers, Up: Mathematical Functions + +4.5 Testing the Sign of Numbers +=============================== + + -- Macro: GSL_SIGN (x) + This macro returns the sign of X. It is defined as `((x) >= 0 ? 1 + : -1)'. Note that with this definition the sign of zero is positive + (regardless of its IEEE sign bit). + + +File: gsl-ref.info, Node: Testing for Odd and Even Numbers, Next: Maximum and Minimum functions, Prev: Testing the Sign of Numbers, Up: Mathematical Functions + +4.6 Testing for Odd and Even Numbers +==================================== + + -- Macro: GSL_IS_ODD (n) + This macro evaluates to 1 if N is odd and 0 if N is even. The + argument N must be of integer type. + + -- Macro: GSL_IS_EVEN (n) + This macro is the opposite of `GSL_IS_ODD(n)'. It evaluates to 1 if + N is even and 0 if N is odd. The argument N must be of integer + type. + + +File: gsl-ref.info, Node: Maximum and Minimum functions, Next: Approximate Comparison of Floating Point Numbers, Prev: Testing for Odd and Even Numbers, Up: Mathematical Functions + +4.7 Maximum and Minimum functions +================================= + + -- Macro: GSL_MAX (a, b) + This macro returns the maximum of A and B. It is defined as `((a) + > (b) ? (a):(b))'. + + -- Macro: GSL_MIN (a, b) + This macro returns the minimum of A and B. It is defined as `((a) + < (b) ? (a):(b))'. + + -- Function: extern inline double GSL_MAX_DBL (double A, double B) + This function returns the maximum of the double precision numbers + A and B using an inline function. The use of a function allows for + type checking of the arguments as an extra safety feature. On + platforms where inline functions are not available the macro + `GSL_MAX' will be automatically substituted. + + -- Function: extern inline double GSL_MIN_DBL (double A, double B) + This function returns the minimum of the double precision numbers + A and B using an inline function. The use of a function allows for + type checking of the arguments as an extra safety feature. On + platforms where inline functions are not available the macro + `GSL_MIN' will be automatically substituted. + + -- Function: extern inline int GSL_MAX_INT (int A, int B) + -- Function: extern inline int GSL_MIN_INT (int A, int B) + These functions return the maximum or minimum of the integers A + and B using an inline function. On platforms where inline + functions are not available the macros `GSL_MAX' or `GSL_MIN' will + be automatically substituted. + + -- Function: extern inline long double GSL_MAX_LDBL (long double A, + long double B) + -- Function: extern inline long double GSL_MIN_LDBL (long double A, + long double B) + These functions return the maximum or minimum of the long doubles A + and B using an inline function. On platforms where inline + functions are not available the macros `GSL_MAX' or `GSL_MIN' will + be automatically substituted. + + +File: gsl-ref.info, Node: Approximate Comparison of Floating Point Numbers, Prev: Maximum and Minimum functions, Up: Mathematical Functions + +4.8 Approximate Comparison of Floating Point Numbers +==================================================== + +It is sometimes useful to be able to compare two floating point numbers +approximately, to allow for rounding and truncation errors. The +following function implements the approximate floating-point comparison +algorithm proposed by D.E. Knuth in Section 4.2.2 of `Seminumerical +Algorithms' (3rd edition). + + -- Function: int gsl_fcmp (double X, double Y, double EPSILON) + This function determines whether x and y are approximately equal + to a relative accuracy EPSILON. + + The relative accuracy is measured using an interval of size 2 + \delta, where \delta = 2^k \epsilon and k is the maximum base-2 + exponent of x and y as computed by the function `frexp'. + + If x and y lie within this interval, they are considered + approximately equal and the function returns 0. Otherwise if x < + y, the function returns -1, or if x > y, the function returns +1. + + The implementation is based on the package `fcmp' by T.C. Belding. + + +File: gsl-ref.info, Node: Complex Numbers, Next: Polynomials, Prev: Mathematical Functions, Up: Top + +5 Complex Numbers +***************** + +The functions described in this chapter provide support for complex +numbers. The algorithms take care to avoid unnecessary intermediate +underflows and overflows, allowing the functions to be evaluated over +as much of the complex plane as possible. + + For multiple-valued functions the branch cuts have been chosen to +follow the conventions of Abramowitz and Stegun in the `Handbook of +Mathematical Functions'. The functions return principal values which are +the same as those in GNU Calc, which in turn are the same as those in +`Common Lisp, The Language (Second Edition)'(1) and the HP-28/48 series +of calculators. + + The complex types are defined in the header file `gsl_complex.h', +while the corresponding complex functions and arithmetic operations are +defined in `gsl_complex_math.h'. + +* Menu: + +* Complex numbers:: +* Properties of complex numbers:: +* Complex arithmetic operators:: +* Elementary Complex Functions:: +* Complex Trigonometric Functions:: +* Inverse Complex Trigonometric Functions:: +* Complex Hyperbolic Functions:: +* Inverse Complex Hyperbolic Functions:: +* Complex Number References and Further Reading:: + + ---------- Footnotes ---------- + + (1) Note that the first edition uses different definitions. + + +File: gsl-ref.info, Node: Complex numbers, Next: Properties of complex numbers, Up: Complex Numbers + +5.1 Complex numbers +=================== + +Complex numbers are represented using the type `gsl_complex'. The +internal representation of this type may vary across platforms and +should not be accessed directly. The functions and macros described +below allow complex numbers to be manipulated in a portable way. + + For reference, the default form of the `gsl_complex' type is given +by the following struct, + + typedef struct + { + double dat[2]; + } gsl_complex; + +The real and imaginary part are stored in contiguous elements of a two +element array. This eliminates any padding between the real and +imaginary parts, `dat[0]' and `dat[1]', allowing the struct to be +mapped correctly onto packed complex arrays. + + -- Function: gsl_complex gsl_complex_rect (double X, double Y) + This function uses the rectangular cartesian components (X,Y) to + return the complex number z = x + i y. + + -- 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. + + +File: gsl-ref.info, Node: Properties of complex numbers, Next: Complex arithmetic operators, Prev: Complex numbers, Up: Complex Numbers + +5.2 Properties of complex numbers +================================= + + -- Function: double gsl_complex_arg (gsl_complex Z) + This function returns the argument of the complex number Z, + \arg(z), where -\pi < \arg(z) <= \pi. + + -- Function: double gsl_complex_abs (gsl_complex Z) + This function returns the magnitude of the complex number Z, |z|. + + -- Function: double gsl_complex_abs2 (gsl_complex Z) + This function returns the squared magnitude of the complex number + Z, |z|^2. + + -- Function: double gsl_complex_logabs (gsl_complex Z) + This function returns the natural logarithm of the magnitude of the + complex number Z, \log|z|. It allows an accurate evaluation of + \log|z| when |z| is close to one. The direct evaluation of + `log(gsl_complex_abs(z))' would lead to a loss of precision in + this case. + + +File: gsl-ref.info, Node: Complex arithmetic operators, Next: Elementary Complex Functions, Prev: Properties of complex numbers, Up: Complex Numbers + +5.3 Complex arithmetic operators +================================ + + -- Function: gsl_complex gsl_complex_add (gsl_complex A, gsl_complex B) + This function returns the sum of the complex numbers A and B, + z=a+b. + + -- Function: gsl_complex gsl_complex_sub (gsl_complex A, gsl_complex B) + This function returns the difference of the complex numbers A and + B, z=a-b. + + -- Function: gsl_complex gsl_complex_mul (gsl_complex A, gsl_complex B) + This function returns the product of the complex numbers A and B, + z=ab. + + -- Function: gsl_complex gsl_complex_div (gsl_complex A, gsl_complex B) + This function returns the quotient of the complex numbers A and B, + z=a/b. + + -- Function: gsl_complex gsl_complex_add_real (gsl_complex A, double X) + This function returns the sum of the complex number A and the real + number X, z=a+x. + + -- Function: gsl_complex gsl_complex_sub_real (gsl_complex A, double X) + This function returns the difference of the complex number A and + the real number X, z=a-x. + + -- Function: gsl_complex gsl_complex_mul_real (gsl_complex A, double X) + This function returns the product of the complex number A and the + real number X, z=ax. + + -- Function: gsl_complex gsl_complex_div_real (gsl_complex A, double X) + This function returns the quotient of the complex number A and the + real number X, z=a/x. + + -- Function: gsl_complex gsl_complex_add_imag (gsl_complex A, double Y) + This function returns the sum of the complex number A and the + imaginary number iY, z=a+iy. + + -- Function: gsl_complex gsl_complex_sub_imag (gsl_complex A, double Y) + This function returns the difference of the complex number A and + the imaginary number iY, z=a-iy. + + -- Function: gsl_complex gsl_complex_mul_imag (gsl_complex A, double Y) + This function returns the product of the complex number A and the + imaginary number iY, z=a*(iy). + + -- Function: gsl_complex gsl_complex_div_imag (gsl_complex A, double Y) + This function returns the quotient of the complex number A and the + imaginary number iY, z=a/(iy). + + -- Function: gsl_complex gsl_complex_conjugate (gsl_complex Z) + This function returns the complex conjugate of the complex number + Z, z^* = x - i y. + + -- Function: gsl_complex gsl_complex_inverse (gsl_complex Z) + This function returns the inverse, or reciprocal, of the complex + number Z, 1/z = (x - i y)/(x^2 + y^2). + + -- Function: gsl_complex gsl_complex_negative (gsl_complex Z) + This function returns the negative of the complex number Z, -z = + (-x) + i(-y). + + +File: gsl-ref.info, Node: Elementary Complex Functions, Next: Complex Trigonometric Functions, Prev: Complex arithmetic operators, Up: Complex Numbers + +5.4 Elementary Complex Functions +================================ + + -- Function: gsl_complex gsl_complex_sqrt (gsl_complex Z) + This function returns the square root of the complex number Z, + \sqrt z. The branch cut is the negative real axis. The result + always lies in the right half of the complex plane. + + -- Function: gsl_complex gsl_complex_sqrt_real (double X) + This function returns the complex square root of the real number + X, where X may be negative. + + -- Function: gsl_complex gsl_complex_pow (gsl_complex Z, gsl_complex A) + The function returns the complex number Z raised to the complex + power A, z^a. This is computed as \exp(\log(z)*a) using complex + logarithms and complex exponentials. + + -- Function: gsl_complex gsl_complex_pow_real (gsl_complex Z, double X) + This function returns the complex number Z raised to the real + power X, z^x. + + -- Function: gsl_complex gsl_complex_exp (gsl_complex Z) + This function returns the complex exponential of the complex number + Z, \exp(z). + + -- Function: gsl_complex gsl_complex_log (gsl_complex Z) + This function returns the complex natural logarithm (base e) of + the complex number Z, \log(z). The branch cut is the negative + real axis. + + -- Function: gsl_complex gsl_complex_log10 (gsl_complex Z) + This function returns the complex base-10 logarithm of the complex + number Z, \log_10 (z). + + -- Function: gsl_complex gsl_complex_log_b (gsl_complex Z, gsl_complex + B) + This function returns the complex base-B logarithm of the complex + number Z, \log_b(z). This quantity is computed as the ratio + \log(z)/\log(b). + + +File: gsl-ref.info, Node: Complex Trigonometric Functions, Next: Inverse Complex Trigonometric Functions, Prev: Elementary Complex Functions, Up: Complex Numbers + +5.5 Complex Trigonometric Functions +=================================== + + -- Function: gsl_complex gsl_complex_sin (gsl_complex Z) + This function returns the complex sine of the complex number Z, + \sin(z) = (\exp(iz) - \exp(-iz))/(2i). + + -- Function: gsl_complex gsl_complex_cos (gsl_complex Z) + This function returns the complex cosine of the complex number Z, + \cos(z) = (\exp(iz) + \exp(-iz))/2. + + -- Function: gsl_complex gsl_complex_tan (gsl_complex Z) + This function returns the complex tangent of the complex number Z, + \tan(z) = \sin(z)/\cos(z). + + -- Function: gsl_complex gsl_complex_sec (gsl_complex Z) + This function returns the complex secant of the complex number Z, + \sec(z) = 1/\cos(z). + + -- Function: gsl_complex gsl_complex_csc (gsl_complex Z) + This function returns the complex cosecant of the complex number Z, + \csc(z) = 1/\sin(z). + + -- Function: gsl_complex gsl_complex_cot (gsl_complex Z) + This function returns the complex cotangent of the complex number + Z, \cot(z) = 1/\tan(z). + + +File: gsl-ref.info, Node: Inverse Complex Trigonometric Functions, Next: Complex Hyperbolic Functions, Prev: Complex Trigonometric Functions, Up: Complex Numbers + +5.6 Inverse Complex Trigonometric Functions +=========================================== + + -- Function: gsl_complex gsl_complex_arcsin (gsl_complex Z) + This function returns the complex arcsine of the complex number Z, + \arcsin(z). The branch cuts are on the real axis, less than -1 and + greater than 1. + + -- Function: gsl_complex gsl_complex_arcsin_real (double Z) + This function returns the complex arcsine of the real number Z, + \arcsin(z). For z between -1 and 1, the function returns a real + value in the range [-\pi/2,\pi/2]. For z less than -1 the result + has a real part of -\pi/2 and a positive imaginary part. For z + greater than 1 the result has a real part of \pi/2 and a negative + imaginary part. + + -- Function: gsl_complex gsl_complex_arccos (gsl_complex Z) + This function returns the complex arccosine of the complex number + Z, \arccos(z). The branch cuts are on the real axis, less than -1 + and greater than 1. + + -- Function: gsl_complex gsl_complex_arccos_real (double Z) + This function returns the complex arccosine of the real number Z, + \arccos(z). For z between -1 and 1, the function returns a real + value in the range [0,\pi]. For z less than -1 the result has a + real part of \pi and a negative imaginary part. For z greater + than 1 the result is purely imaginary and positive. + + -- Function: gsl_complex gsl_complex_arctan (gsl_complex Z) + This function returns the complex arctangent of the complex number + Z, \arctan(z). The branch cuts are on the imaginary axis, below -i + and above i. + + -- Function: gsl_complex gsl_complex_arcsec (gsl_complex Z) + This function returns the complex arcsecant of the complex number + Z, \arcsec(z) = \arccos(1/z). + + -- Function: gsl_complex gsl_complex_arcsec_real (double Z) + This function returns the complex arcsecant of the real number Z, + \arcsec(z) = \arccos(1/z). + + -- Function: gsl_complex gsl_complex_arccsc (gsl_complex Z) + This function returns the complex arccosecant of the complex + number Z, \arccsc(z) = \arcsin(1/z). + + -- Function: gsl_complex gsl_complex_arccsc_real (double Z) + This function returns the complex arccosecant of the real number Z, + \arccsc(z) = \arcsin(1/z). + + -- Function: gsl_complex gsl_complex_arccot (gsl_complex Z) + This function returns the complex arccotangent of the complex + number Z, \arccot(z) = \arctan(1/z). + + +File: gsl-ref.info, Node: Complex Hyperbolic Functions, Next: Inverse Complex Hyperbolic Functions, Prev: Inverse Complex Trigonometric Functions, Up: Complex Numbers + +5.7 Complex Hyperbolic Functions +================================ + + -- Function: gsl_complex gsl_complex_sinh (gsl_complex Z) + This function returns the complex hyperbolic sine of the complex + number Z, \sinh(z) = (\exp(z) - \exp(-z))/2. + + -- Function: gsl_complex gsl_complex_cosh (gsl_complex Z) + This function returns the complex hyperbolic cosine of the complex + number Z, \cosh(z) = (\exp(z) + \exp(-z))/2. + + -- Function: gsl_complex gsl_complex_tanh (gsl_complex Z) + This function returns the complex hyperbolic tangent of the + complex number Z, \tanh(z) = \sinh(z)/\cosh(z). + + -- Function: gsl_complex gsl_complex_sech (gsl_complex Z) + This function returns the complex hyperbolic secant of the complex + number Z, \sech(z) = 1/\cosh(z). + + -- Function: gsl_complex gsl_complex_csch (gsl_complex Z) + This function returns the complex hyperbolic cosecant of the + complex number Z, \csch(z) = 1/\sinh(z). + + -- Function: gsl_complex gsl_complex_coth (gsl_complex Z) + This function returns the complex hyperbolic cotangent of the + complex number Z, \coth(z) = 1/\tanh(z). + + +File: gsl-ref.info, Node: Inverse Complex Hyperbolic Functions, Next: Complex Number References and Further Reading, Prev: Complex Hyperbolic Functions, Up: Complex Numbers + +5.8 Inverse Complex Hyperbolic Functions +======================================== + + -- Function: gsl_complex gsl_complex_arcsinh (gsl_complex Z) + This function returns the complex hyperbolic arcsine of the + complex number Z, \arcsinh(z). The branch cuts are on the + imaginary axis, below -i and above i. + + -- Function: gsl_complex gsl_complex_arccosh (gsl_complex Z) + This function returns the complex hyperbolic arccosine of the + complex number Z, \arccosh(z). The branch cut is on the real + axis, less than 1. Note that in this case we use the negative + square root in formula 4.6.21 of Abramowitz & Stegun giving + \arccosh(z)=\log(z-\sqrt{z^2-1}). + + -- Function: gsl_complex gsl_complex_arccosh_real (double Z) + This function returns the complex hyperbolic arccosine of the real + number Z, \arccosh(z). + + -- Function: gsl_complex gsl_complex_arctanh (gsl_complex Z) + This function returns the complex hyperbolic arctangent of the + complex number Z, \arctanh(z). The branch cuts are on the real + axis, less than -1 and greater than 1. + + -- Function: gsl_complex gsl_complex_arctanh_real (double Z) + This function returns the complex hyperbolic arctangent of the real + number Z, \arctanh(z). + + -- Function: gsl_complex gsl_complex_arcsech (gsl_complex Z) + This function returns the complex hyperbolic arcsecant of the + complex number Z, \arcsech(z) = \arccosh(1/z). + + -- Function: gsl_complex gsl_complex_arccsch (gsl_complex Z) + This function returns the complex hyperbolic arccosecant of the + complex number Z, \arccsch(z) = \arcsin(1/z). + + -- Function: gsl_complex gsl_complex_arccoth (gsl_complex Z) + This function returns the complex hyperbolic arccotangent of the + complex number Z, \arccoth(z) = \arctanh(1/z). + + +File: gsl-ref.info, Node: Complex Number References and Further Reading, Prev: Inverse Complex Hyperbolic Functions, Up: Complex Numbers + +5.9 References and Further Reading +================================== + +The implementations of the elementary and trigonometric functions are +based on the following papers, + + T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang, + "Implementing Complex Elementary Functions Using Exception + Handling", `ACM Transactions on Mathematical Software', Volume 20 + (1994), pp 215-244, Corrigenda, p553 + + T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang, + "Implementing the complex arcsin and arccosine functions using + exception handling", `ACM Transactions on Mathematical Software', + Volume 23 (1997) pp 299-335 + +The general formulas and details of branch cuts can be found in the +following books, + + Abramowitz and Stegun, `Handbook of Mathematical Functions', + "Circular Functions in Terms of Real and Imaginary Parts", Formulas + 4.3.55-58, "Inverse Circular Functions in Terms of Real and + Imaginary Parts", Formulas 4.4.37-39, "Hyperbolic Functions in + Terms of Real and Imaginary Parts", Formulas 4.5.49-52, "Inverse + Hyperbolic Functions--relation to Inverse Circular Functions", + Formulas 4.6.14-19. + + Dave Gillespie, `Calc Manual', Free Software Foundation, ISBN + 1-882114-18-3 + + +File: gsl-ref.info, Node: Polynomials, Next: Special Functions, Prev: Complex Numbers, Up: Top + +6 Polynomials +************* + +This chapter describes functions for evaluating and solving polynomials. +There are routines for finding real and complex roots of quadratic and +cubic equations using analytic methods. An iterative polynomial solver +is also available for finding the roots of general polynomials with real +coefficients (of any order). The functions are declared in the header +file `gsl_poly.h'. + +* Menu: + +* Polynomial Evaluation:: +* Divided Difference Representation of Polynomials:: +* Quadratic Equations:: +* Cubic Equations:: +* General Polynomial Equations:: +* Roots of Polynomials Examples:: +* Roots of Polynomials References and Further Reading:: + + +File: gsl-ref.info, Node: Polynomial Evaluation, Next: Divided Difference Representation of Polynomials, Up: Polynomials + +6.1 Polynomial Evaluation +========================= + + -- Function: double gsl_poly_eval (const double C[], const int LEN, + const double X) + This function evaluates the polynomial c[0] + c[1] x + c[2] x^2 + + \dots + c[len-1] x^{len-1} using Horner's method for stability. + The function is inlined when possible. + + +File: gsl-ref.info, Node: Divided Difference Representation of Polynomials, Next: Quadratic Equations, Prev: Polynomial Evaluation, Up: Polynomials + +6.2 Divided Difference Representation of Polynomials +==================================================== + +The functions described here manipulate polynomials stored in Newton's +divided-difference representation. The use of divided-differences is +described in Abramowitz & Stegun sections 25.1.4 and 25.2.26. + + -- Function: int gsl_poly_dd_init (double DD[], const double XA[], + const double YA[], size_t SIZE) + This function computes a divided-difference representation of the + interpolating polynomial for the points (XA, YA) stored in the + arrays XA and YA of length SIZE. On output the + divided-differences of (XA,YA) are stored in the array DD, also of + length SIZE. + + -- Function: double gsl_poly_dd_eval (const double DD[], const double + XA[], const size_t SIZE, const double X) + This function evaluates the polynomial stored in + divided-difference form in the arrays DD and XA of length SIZE at + the point X. + + -- 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. + + +File: gsl-ref.info, Node: Quadratic Equations, Next: Cubic Equations, Prev: Divided Difference Representation of Polynomials, Up: Polynomials + +6.3 Quadratic Equations +======================= + + -- Function: int gsl_poly_solve_quadratic (double A, double B, double + C, double * X0, double * X1) + This function finds the real roots of the quadratic equation, + + a x^2 + b x + c = 0 + + The number of real roots (either zero, one or two) is returned, and + their locations are stored in X0 and X1. If no real roots are + found then X0 and X1 are not modified. If one real root is found + (i.e. if a=0) then it is stored in X0. When two real roots are + found they are stored in X0 and X1 in ascending order. The case + of coincident roots is not considered special. For example + (x-1)^2=0 will have two roots, which happen to have exactly equal + values. + + The number of roots found depends on the sign of the discriminant + b^2 - 4 a c. This will be subject to rounding and cancellation + errors when computed in double precision, and will also be subject + to errors if the coefficients of the polynomial are inexact. + These errors may cause a discrete change in the number of roots. + However, for polynomials with small integer coefficients the + discriminant can always be computed exactly. + + + -- Function: int gsl_poly_complex_solve_quadratic (double A, double B, + double C, gsl_complex * Z0, gsl_complex * Z1) + This function finds the complex roots of the quadratic equation, + + a z^2 + b z + c = 0 + + The number of complex roots is returned (either one or two) and the + locations of the roots are stored in Z0 and Z1. The roots are + returned in ascending order, sorted first by their real components + and then by their imaginary components. If only one real root is + found (i.e. if a=0) then it is stored in Z0. + + + +File: gsl-ref.info, Node: Cubic Equations, Next: General Polynomial Equations, Prev: Quadratic Equations, Up: Polynomials + +6.4 Cubic Equations +=================== + + -- Function: int gsl_poly_solve_cubic (double A, double B, double C, + double * X0, double * X1, double * X2) + This function finds the real roots of the cubic equation, + + x^3 + a x^2 + b x + c = 0 + + with a leading coefficient of unity. The number of real roots + (either one or three) is returned, and their locations are stored + in X0, X1 and X2. If one real root is found then only X0 is + modified. When three real roots are found they are stored in X0, + X1 and X2 in ascending order. The case of coincident roots is not + considered special. For example, the equation (x-1)^3=0 will have + three roots with exactly equal values. + + + -- 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. + + + +File: gsl-ref.info, Node: General Polynomial Equations, Next: Roots of Polynomials Examples, Prev: Cubic Equations, Up: Polynomials + +6.5 General Polynomial Equations +================================ + +The roots of polynomial equations cannot be found analytically beyond +the special cases of the quadratic, cubic and quartic equation. The +algorithm described in this section uses an iterative method to find the +approximate locations of roots of higher order polynomials. + + -- Function: gsl_poly_complex_workspace * +gsl_poly_complex_workspace_alloc (size_t N) + This function allocates space for a `gsl_poly_complex_workspace' + struct and a workspace suitable for solving a polynomial with N + coefficients using the routine `gsl_poly_complex_solve'. + + The function returns a pointer to the newly allocated + `gsl_poly_complex_workspace' if no errors were detected, and a null + pointer in the case of error. + + -- Function: void gsl_poly_complex_workspace_free + (gsl_poly_complex_workspace * W) + This function frees all the memory associated with the workspace W. + + -- Function: int gsl_poly_complex_solve (const double * A, size_t N, + gsl_poly_complex_workspace * W, gsl_complex_packed_ptr Z) + This function computes the roots of the general polynomial P(x) = + a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using balanced-QR + reduction of the companion matrix. The parameter N specifies the + length of the coefficient array. The coefficient of the highest + order term must be non-zero. The function requires a workspace W + of the appropriate size. The n-1 roots are returned in the packed + complex array Z of length 2(n-1), alternating real and imaginary + parts. + + The function returns `GSL_SUCCESS' if all the roots are found. If + the QR reduction does not converge, the error handler is invoked + with an error code of `GSL_EFAILED'. Note that due to finite + precision, roots of higher multiplicity are returned as a cluster + of simple roots with reduced accuracy. The solution of + polynomials with higher-order roots requires specialized + algorithms that take the multiplicity structure into account (see + e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical + Software, Volume 30, Issue 2 (2004), pp 218-236). + + +File: gsl-ref.info, Node: Roots of Polynomials Examples, Next: Roots of Polynomials References and Further Reading, Prev: General Polynomial Equations, Up: Polynomials + +6.6 Examples +============ + +To demonstrate the use of the general polynomial solver we will take the +polynomial P(x) = x^5 - 1 which has the following roots, + + 1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5} + +The following program will find these roots. + + #include <stdio.h> + #include <gsl/gsl_poly.h> + + int + main (void) + { + int i; + /* coefficients of P(x) = -1 + x^5 */ + double a[6] = { -1, 0, 0, 0, 0, 1 }; + double z[10]; + + gsl_poly_complex_workspace * w + = gsl_poly_complex_workspace_alloc (6); + + gsl_poly_complex_solve (a, 6, w, z); + + gsl_poly_complex_workspace_free (w); + + for (i = 0; i < 5; i++) + { + printf ("z%d = %+.18f %+.18f\n", + i, z[2*i], z[2*i+1]); + } + + return 0; + } + +The output of the program is, + + $ ./a.out + z0 = -0.809016994374947451 +0.587785252292473137 + z1 = -0.809016994374947451 -0.587785252292473137 + z2 = +0.309016994374947451 +0.951056516295153642 + z3 = +0.309016994374947451 -0.951056516295153642 + z4 = +1.000000000000000000 +0.000000000000000000 + +which agrees with the analytic result, z_n = \exp(2 \pi n i/5). + + +File: gsl-ref.info, Node: Roots of Polynomials References and Further Reading, Prev: Roots of Polynomials Examples, Up: Polynomials + +6.7 References and Further Reading +================================== + +The balanced-QR method and its error analysis are described in the +following papers, + + R.S. Martin, G. Peters and J.H. Wilkinson, "The QR Algorithm for + Real Hessenberg Matrices", `Numerische Mathematik', 14 (1970), + 219-231. + + B.N. Parlett and C. Reinsch, "Balancing a Matrix for Calculation of + Eigenvalues and Eigenvectors", `Numerische Mathematik', 13 (1969), + 293-304. + + A. Edelman and H. Murakami, "Polynomial roots from companion matrix + eigenvalues", `Mathematics of Computation', Vol. 64, No. 210 + (1995), 763-776. + +The formulas for divided differences are given in Abramowitz and Stegun, + + Abramowitz and Stegun, `Handbook of Mathematical Functions', + Sections 25.1.4 and 25.2.26. + + +File: gsl-ref.info, Node: Special Functions, Next: Vectors and Matrices, Prev: Polynomials, Up: Top + +7 Special Functions +******************* + +This chapter describes the GSL special function library. The library +includes routines for calculating the values of Airy functions, Bessel +functions, Clausen functions, Coulomb wave functions, Coupling +coefficients, the Dawson function, Debye functions, Dilogarithms, +Elliptic integrals, Jacobi elliptic functions, Error functions, +Exponential integrals, Fermi-Dirac functions, Gamma functions, +Gegenbauer functions, Hypergeometric functions, Laguerre functions, +Legendre functions and Spherical Harmonics, the Psi (Digamma) Function, +Synchrotron functions, Transport functions, Trigonometric functions and +Zeta functions. Each routine also computes an estimate of the numerical +error in the calculated value of the function. + + The functions in this chapter are declared in individual header +files, such as `gsl_sf_airy.h', `gsl_sf_bessel.h', etc. The complete +set of header files can be included using the file `gsl_sf.h'. + +* Menu: + +* Special Function Usage:: +* The gsl_sf_result struct:: +* Special Function Modes:: +* Airy Functions and Derivatives:: +* Bessel Functions:: +* Clausen Functions:: +* Coulomb Functions:: +* Coupling Coefficients:: +* Dawson Function:: +* Debye Functions:: +* Dilogarithm:: +* Elementary Operations:: +* Elliptic Integrals:: +* Elliptic Functions (Jacobi):: +* Error Functions:: +* Exponential Functions:: +* Exponential Integrals:: +* Fermi-Dirac Function:: +* Gamma and Beta Functions:: +* Gegenbauer Functions:: +* Hypergeometric Functions:: +* Laguerre Functions:: +* Lambert W Functions:: +* Legendre Functions and Spherical Harmonics:: +* Logarithm and Related Functions:: +* Mathieu Functions:: +* Power Function:: +* Psi (Digamma) Function:: +* Synchrotron Functions:: +* Transport Functions:: +* Trigonometric Functions:: +* Zeta Functions:: +* Special Functions Examples:: +* Special Functions References and Further Reading:: + + +File: gsl-ref.info, Node: Special Function Usage, Next: The gsl_sf_result struct, Up: Special Functions + +7.1 Usage +========= + +The special functions are available in two calling conventions, a +"natural form" which returns the numerical value of the function and an +"error-handling form" which returns an error code. The two types of +function provide alternative ways of accessing the same underlying code. + + The "natural form" returns only the value of the function and can be +used directly in mathematical expressions. For example, the following +function call will compute the value of the Bessel function J_0(x), + + double y = gsl_sf_bessel_J0 (x); + +There is no way to access an error code or to estimate the error using +this method. To allow access to this information the alternative +error-handling form stores the value and error in a modifiable argument, + + gsl_sf_result result; + int status = gsl_sf_bessel_J0_e (x, &result); + +The error-handling functions have the suffix `_e'. The returned status +value indicates error conditions such as overflow, underflow or loss of +precision. If there are no errors the error-handling functions return +`GSL_SUCCESS'. + + +File: gsl-ref.info, Node: The gsl_sf_result struct, Next: Special Function Modes, Prev: Special Function Usage, Up: Special Functions + +7.2 The gsl_sf_result struct +============================ + +The error handling form of the special functions always calculate an +error estimate along with the value of the result. Therefore, +structures are provided for amalgamating a value and error estimate. +These structures are declared in the header file `gsl_sf_result.h'. + + The `gsl_sf_result' struct contains value and error fields. + + typedef struct + { + double val; + double err; + } gsl_sf_result; + +The field VAL contains the value and the field ERR contains an estimate +of the absolute error in the value. + + In some cases, an overflow or underflow can be detected and handled +by a function. In this case, it may be possible to return a scaling +exponent as well as an error/value pair in order to save the result from +exceeding the dynamic range of the built-in types. The +`gsl_sf_result_e10' struct contains value and error fields as well as +an exponent field such that the actual result is obtained as `result * +10^(e10)'. + + typedef struct + { + double val; + double err; + int e10; + } gsl_sf_result_e10; + + +File: gsl-ref.info, Node: Special Function Modes, Next: Airy Functions and Derivatives, Prev: The gsl_sf_result struct, Up: Special Functions + +7.3 Modes +========= + +The goal of the library is to achieve double precision accuracy wherever +possible. However the cost of evaluating some special functions to +double precision can be significant, particularly where very high order +terms are required. In these cases a `mode' argument allows the +accuracy of the function to be reduced in order to improve performance. +The following precision levels are available for the mode argument, + +`GSL_PREC_DOUBLE' + Double-precision, a relative accuracy of approximately 2 * 10^-16. + +`GSL_PREC_SINGLE' + Single-precision, a relative accuracy of approximately 10^-7. + +`GSL_PREC_APPROX' + Approximate values, a relative accuracy of approximately 5 * 10^-4. + +The approximate mode provides the fastest evaluation at the lowest +accuracy. + + +File: gsl-ref.info, Node: Airy Functions and Derivatives, Next: Bessel Functions, Prev: Special Function Modes, Up: Special Functions + +7.4 Airy Functions and Derivatives +================================== + +The Airy functions Ai(x) and Bi(x) are defined by the integral +representations, + + Ai(x) = (1/\pi) \int_0^\infty \cos((1/3) t^3 + xt) dt + Bi(x) = (1/\pi) \int_0^\infty (e^(-(1/3) t^3) + \sin((1/3) t^3 + xt)) dt + +For further information see Abramowitz & Stegun, Section 10.4. The Airy +functions are defined in the header file `gsl_sf_airy.h'. + +* Menu: + +* Airy Functions:: +* Derivatives of Airy Functions:: +* Zeros of Airy Functions:: +* Zeros of Derivatives of Airy Functions:: + + +File: gsl-ref.info, Node: Airy Functions, Next: Derivatives of Airy Functions, Up: Airy Functions and Derivatives + +7.4.1 Airy Functions +-------------------- + + -- Function: double gsl_sf_airy_Ai (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Ai_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the Airy function Ai(x) with an accuracy + specified by MODE. + + -- Function: double gsl_sf_airy_Bi (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Bi_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the Airy function Bi(x) with an accuracy + specified by MODE. + + -- Function: double gsl_sf_airy_Ai_scaled (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Ai_scaled_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute a scaled version of the Airy function + S_A(x) Ai(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3) + x^(3/2)), and is 1 for x<0. + + -- Function: double gsl_sf_airy_Bi_scaled (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Bi_scaled_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute a scaled version of the Airy function + S_B(x) Bi(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) + x^(3/2)), and is 1 for x<0. + + +File: gsl-ref.info, Node: Derivatives of Airy Functions, Next: Zeros of Airy Functions, Prev: Airy Functions, Up: Airy Functions and Derivatives + +7.4.2 Derivatives of Airy Functions +----------------------------------- + + -- Function: double gsl_sf_airy_Ai_deriv (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Ai_deriv_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the Airy function derivative Ai'(x) with an + accuracy specified by MODE. + + -- Function: double gsl_sf_airy_Bi_deriv (double X, gsl_mode_t MODE) + -- Function: int gsl_sf_airy_Bi_deriv_e (double X, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the Airy function derivative Bi'(x) with an + accuracy specified by MODE. + + -- Function: double gsl_sf_airy_Ai_deriv_scaled (double X, gsl_mode_t + MODE) + -- Function: int gsl_sf_airy_Ai_deriv_scaled_e (double X, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the scaled Airy function derivative S_A(x) + Ai'(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3) + x^(3/2)), and is 1 for x<0. + + -- Function: double gsl_sf_airy_Bi_deriv_scaled (double X, gsl_mode_t + MODE) + -- Function: int gsl_sf_airy_Bi_deriv_scaled_e (double X, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the scaled Airy function derivative S_B(x) + Bi'(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)), + and is 1 for x<0. + + +File: gsl-ref.info, Node: Zeros of Airy Functions, Next: Zeros of Derivatives of Airy Functions, Prev: Derivatives of Airy Functions, Up: Airy Functions and Derivatives + +7.4.3 Zeros of Airy Functions +----------------------------- + + -- Function: double gsl_sf_airy_zero_Ai (unsigned int S) + -- Function: int gsl_sf_airy_zero_Ai_e (unsigned int S, gsl_sf_result + * RESULT) + These routines compute the location of the S-th zero of the Airy + function Ai(x). + + -- Function: double gsl_sf_airy_zero_Bi (unsigned int S) + -- Function: int gsl_sf_airy_zero_Bi_e (unsigned int S, gsl_sf_result + * RESULT) + These routines compute the location of the S-th zero of the Airy + function Bi(x). + + +File: gsl-ref.info, Node: Zeros of Derivatives of Airy Functions, Prev: Zeros of Airy Functions, Up: Airy Functions and Derivatives + +7.4.4 Zeros of Derivatives of Airy Functions +-------------------------------------------- + + -- Function: double gsl_sf_airy_zero_Ai_deriv (unsigned int S) + -- Function: int gsl_sf_airy_zero_Ai_deriv_e (unsigned int S, + gsl_sf_result * RESULT) + These routines compute the location of the S-th zero of the Airy + function derivative Ai'(x). + + -- Function: double gsl_sf_airy_zero_Bi_deriv (unsigned int S) + -- Function: int gsl_sf_airy_zero_Bi_deriv_e (unsigned int S, + gsl_sf_result * RESULT) + These routines compute the location of the S-th zero of the Airy + function derivative Bi'(x). + + +File: gsl-ref.info, Node: Bessel Functions, Next: Clausen Functions, Prev: Airy Functions and Derivatives, Up: Special Functions + +7.5 Bessel Functions +==================== + +The routines described in this section compute the Cylindrical Bessel +functions J_n(x), Y_n(x), Modified cylindrical Bessel functions I_n(x), +K_n(x), Spherical Bessel functions j_l(x), y_l(x), and Modified +Spherical Bessel functions i_l(x), k_l(x). For more information see +Abramowitz & Stegun, Chapters 9 and 10. The Bessel functions are +defined in the header file `gsl_sf_bessel.h'. + +* Menu: + +* Regular Cylindrical Bessel Functions:: +* Irregular Cylindrical Bessel Functions:: +* Regular Modified Cylindrical Bessel Functions:: +* Irregular Modified Cylindrical Bessel Functions:: +* Regular Spherical Bessel Functions:: +* Irregular Spherical Bessel Functions:: +* Regular Modified Spherical Bessel Functions:: +* Irregular Modified Spherical Bessel Functions:: +* Regular Bessel Function - Fractional Order:: +* Irregular Bessel Functions - Fractional Order:: +* Regular Modified Bessel Functions - Fractional Order:: +* Irregular Modified Bessel Functions - Fractional Order:: +* Zeros of Regular Bessel Functions:: + + +File: gsl-ref.info, Node: Regular Cylindrical Bessel Functions, Next: Irregular Cylindrical Bessel Functions, Up: Bessel Functions + +7.5.1 Regular Cylindrical Bessel Functions +------------------------------------------ + + -- Function: double gsl_sf_bessel_J0 (double X) + -- Function: int gsl_sf_bessel_J0_e (double X, gsl_sf_result * RESULT) + These routines compute the regular cylindrical Bessel function of + zeroth order, J_0(x). + + -- Function: double gsl_sf_bessel_J1 (double X) + -- Function: int gsl_sf_bessel_J1_e (double X, gsl_sf_result * RESULT) + These routines compute the regular cylindrical Bessel function of + first order, J_1(x). + + -- Function: double gsl_sf_bessel_Jn (int N, double X) + -- Function: int gsl_sf_bessel_Jn_e (int N, double X, gsl_sf_result * + RESULT) + These routines compute the regular cylindrical Bessel function of + order N, J_n(x). + + -- Function: int gsl_sf_bessel_Jn_array (int NMIN, int NMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the regular cylindrical Bessel + functions J_n(x) for n from NMIN to NMAX inclusive, storing the + results in the array RESULT_ARRAY. The values are computed using + recurrence relations for efficiency, and therefore may differ + slightly from the exact values. + + +File: gsl-ref.info, Node: Irregular Cylindrical Bessel Functions, Next: Regular Modified Cylindrical Bessel Functions, Prev: Regular Cylindrical Bessel Functions, Up: Bessel Functions + +7.5.2 Irregular Cylindrical Bessel Functions +-------------------------------------------- + + -- Function: double gsl_sf_bessel_Y0 (double X) + -- Function: int gsl_sf_bessel_Y0_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular cylindrical Bessel function + of zeroth order, Y_0(x), for x>0. + + -- Function: double gsl_sf_bessel_Y1 (double X) + -- Function: int gsl_sf_bessel_Y1_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular cylindrical Bessel function + of first order, Y_1(x), for x>0. + + -- Function: double gsl_sf_bessel_Yn (int N,double X) + -- Function: int gsl_sf_bessel_Yn_e (int N,double X, gsl_sf_result * + RESULT) + These routines compute the irregular cylindrical Bessel function of + order N, Y_n(x), for x>0. + + -- Function: int gsl_sf_bessel_Yn_array (int NMIN, int NMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the irregular cylindrical + Bessel functions Y_n(x) for n from NMIN to NMAX inclusive, storing + the results in the array RESULT_ARRAY. The domain of the function + is x>0. The values are computed using recurrence relations for + efficiency, and therefore may differ slightly from the exact + values. + + +File: gsl-ref.info, Node: Regular Modified Cylindrical Bessel Functions, Next: Irregular Modified Cylindrical Bessel Functions, Prev: Irregular Cylindrical Bessel Functions, Up: Bessel Functions + +7.5.3 Regular Modified Cylindrical Bessel Functions +--------------------------------------------------- + + -- Function: double gsl_sf_bessel_I0 (double X) + -- Function: int gsl_sf_bessel_I0_e (double X, gsl_sf_result * RESULT) + These routines compute the regular modified cylindrical Bessel + function of zeroth order, I_0(x). + + -- Function: double gsl_sf_bessel_I1 (double X) + -- Function: int gsl_sf_bessel_I1_e (double X, gsl_sf_result * RESULT) + These routines compute the regular modified cylindrical Bessel + function of first order, I_1(x). + + -- Function: double gsl_sf_bessel_In (int N, double X) + -- Function: int gsl_sf_bessel_In_e (int N, double X, gsl_sf_result * + RESULT) + These routines compute the regular modified cylindrical Bessel + function of order N, I_n(x). + + -- Function: int gsl_sf_bessel_In_array (int NMIN, int NMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the regular modified + cylindrical Bessel functions I_n(x) for n from NMIN to NMAX + inclusive, storing the results in the array RESULT_ARRAY. The + start of the range NMIN must be positive or zero. The values are + computed using recurrence relations for efficiency, and therefore + may differ slightly from the exact values. + + -- Function: double gsl_sf_bessel_I0_scaled (double X) + -- Function: int gsl_sf_bessel_I0_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled regular modified cylindrical + Bessel function of zeroth order \exp(-|x|) I_0(x). + + -- Function: double gsl_sf_bessel_I1_scaled (double X) + -- Function: int gsl_sf_bessel_I1_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled regular modified cylindrical + Bessel function of first order \exp(-|x|) I_1(x). + + -- Function: double gsl_sf_bessel_In_scaled (int N, double X) + -- Function: int gsl_sf_bessel_In_scaled_e (int N, double X, + gsl_sf_result * RESULT) + These routines compute the scaled regular modified cylindrical + Bessel function of order N, \exp(-|x|) I_n(x) + + -- Function: int gsl_sf_bessel_In_scaled_array (int NMIN, int NMAX, + double X, double RESULT_ARRAY[]) + This routine computes the values of the scaled regular cylindrical + Bessel functions \exp(-|x|) I_n(x) for n from NMIN to NMAX + inclusive, storing the results in the array RESULT_ARRAY. The + start of the range NMIN must be positive or zero. The values are + computed using recurrence relations for efficiency, and therefore + may differ slightly from the exact values. + + +File: gsl-ref.info, Node: Irregular Modified Cylindrical Bessel Functions, Next: Regular Spherical Bessel Functions, Prev: Regular Modified Cylindrical Bessel Functions, Up: Bessel Functions + +7.5.4 Irregular Modified Cylindrical Bessel Functions +----------------------------------------------------- + + -- Function: double gsl_sf_bessel_K0 (double X) + -- Function: int gsl_sf_bessel_K0_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular modified cylindrical Bessel + function of zeroth order, K_0(x), for x > 0. + + -- Function: double gsl_sf_bessel_K1 (double X) + -- Function: int gsl_sf_bessel_K1_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular modified cylindrical Bessel + function of first order, K_1(x), for x > 0. + + -- Function: double gsl_sf_bessel_Kn (int N, double X) + -- Function: int gsl_sf_bessel_Kn_e (int N, double X, gsl_sf_result * + RESULT) + These routines compute the irregular modified cylindrical Bessel + function of order N, K_n(x), for x > 0. + + -- Function: int gsl_sf_bessel_Kn_array (int NMIN, int NMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the irregular modified + cylindrical Bessel functions K_n(x) for n from NMIN to NMAX + inclusive, storing the results in the array RESULT_ARRAY. The + start of the range NMIN must be positive or zero. The domain of + the function is x>0. The values are computed using recurrence + relations for efficiency, and therefore may differ slightly from + the exact values. + + -- Function: double gsl_sf_bessel_K0_scaled (double X) + -- Function: int gsl_sf_bessel_K0_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled irregular modified cylindrical + Bessel function of zeroth order \exp(x) K_0(x) for x>0. + + -- Function: double gsl_sf_bessel_K1_scaled (double X) + -- Function: int gsl_sf_bessel_K1_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled irregular modified cylindrical + Bessel function of first order \exp(x) K_1(x) for x>0. + + -- Function: double gsl_sf_bessel_Kn_scaled (int N, double X) + -- Function: int gsl_sf_bessel_Kn_scaled_e (int N, double X, + gsl_sf_result * RESULT) + These routines compute the scaled irregular modified cylindrical + Bessel function of order N, \exp(x) K_n(x), for x>0. + + -- Function: int gsl_sf_bessel_Kn_scaled_array (int NMIN, int NMAX, + double X, double RESULT_ARRAY[]) + This routine computes the values of the scaled irregular + cylindrical Bessel functions \exp(x) K_n(x) for n from NMIN to + NMAX inclusive, storing the results in the array RESULT_ARRAY. The + start of the range NMIN must be positive or zero. The domain of + the function is x>0. The values are computed using recurrence + relations for efficiency, and therefore may differ slightly from + the exact values. + + +File: gsl-ref.info, Node: Regular Spherical Bessel Functions, Next: Irregular Spherical Bessel Functions, Prev: Irregular Modified Cylindrical Bessel Functions, Up: Bessel Functions + +7.5.5 Regular Spherical Bessel Functions +---------------------------------------- + + -- Function: double gsl_sf_bessel_j0 (double X) + -- Function: int gsl_sf_bessel_j0_e (double X, gsl_sf_result * RESULT) + These routines compute the regular spherical Bessel function of + zeroth order, j_0(x) = \sin(x)/x. + + -- Function: double gsl_sf_bessel_j1 (double X) + -- Function: int gsl_sf_bessel_j1_e (double X, gsl_sf_result * RESULT) + These routines compute the regular spherical Bessel function of + first order, j_1(x) = (\sin(x)/x - \cos(x))/x. + + -- Function: double gsl_sf_bessel_j2 (double X) + -- Function: int gsl_sf_bessel_j2_e (double X, gsl_sf_result * RESULT) + These routines compute the regular spherical Bessel function of + second order, j_2(x) = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x. + + -- Function: double gsl_sf_bessel_jl (int L, double X) + -- Function: int gsl_sf_bessel_jl_e (int L, double X, gsl_sf_result * + RESULT) + These routines compute the regular spherical Bessel function of + order L, j_l(x), for l >= 0 and x >= 0. + + -- Function: int gsl_sf_bessel_jl_array (int LMAX, double X, double + RESULT_ARRAY[]) + This routine computes the values of the regular spherical Bessel + functions j_l(x) for l from 0 to LMAX inclusive for lmax >= 0 and + x >= 0, storing the results in the array RESULT_ARRAY. The values + are computed using recurrence relations for efficiency, and + therefore may differ slightly from the exact values. + + -- Function: int gsl_sf_bessel_jl_steed_array (int LMAX, double X, + double * JL_X_ARRAY) + This routine uses Steed's method to compute the values of the + regular spherical Bessel functions j_l(x) for l from 0 to LMAX + inclusive for lmax >= 0 and x >= 0, storing the results in the + array RESULT_ARRAY. The Steed/Barnett algorithm is described in + `Comp. Phys. Comm.' 21, 297 (1981). Steed's method is more stable + than the recurrence used in the other functions but is also slower. + + +File: gsl-ref.info, Node: Irregular Spherical Bessel Functions, Next: Regular Modified Spherical Bessel Functions, Prev: Regular Spherical Bessel Functions, Up: Bessel Functions + +7.5.6 Irregular Spherical Bessel Functions +------------------------------------------ + + -- Function: double gsl_sf_bessel_y0 (double X) + -- Function: int gsl_sf_bessel_y0_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular spherical Bessel function of + zeroth order, y_0(x) = -\cos(x)/x. + + -- Function: double gsl_sf_bessel_y1 (double X) + -- Function: int gsl_sf_bessel_y1_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular spherical Bessel function of + first order, y_1(x) = -(\cos(x)/x + \sin(x))/x. + + -- Function: double gsl_sf_bessel_y2 (double X) + -- Function: int gsl_sf_bessel_y2_e (double X, gsl_sf_result * RESULT) + These routines compute the irregular spherical Bessel function of + second order, y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x). + + -- Function: double gsl_sf_bessel_yl (int L, double X) + -- Function: int gsl_sf_bessel_yl_e (int L, double X, gsl_sf_result * + RESULT) + These routines compute the irregular spherical Bessel function of + order L, y_l(x), for l >= 0. + + -- Function: int gsl_sf_bessel_yl_array (int LMAX, double X, double + RESULT_ARRAY[]) + This routine computes the values of the irregular spherical Bessel + functions y_l(x) for l from 0 to LMAX inclusive for lmax >= 0, + storing the results in the array RESULT_ARRAY. The values are + computed using recurrence relations for efficiency, and therefore + may differ slightly from the exact values. + + +File: gsl-ref.info, Node: Regular Modified Spherical Bessel Functions, Next: Irregular Modified Spherical Bessel Functions, Prev: Irregular Spherical Bessel Functions, Up: Bessel Functions + +7.5.7 Regular Modified Spherical Bessel Functions +------------------------------------------------- + +The regular modified spherical Bessel functions i_l(x) are related to +the modified Bessel functions of fractional order, i_l(x) = +\sqrt{\pi/(2x)} I_{l+1/2}(x) + + -- Function: double gsl_sf_bessel_i0_scaled (double X) + -- Function: int gsl_sf_bessel_i0_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled regular modified spherical Bessel + function of zeroth order, \exp(-|x|) i_0(x). + + -- Function: double gsl_sf_bessel_i1_scaled (double X) + -- Function: int gsl_sf_bessel_i1_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled regular modified spherical Bessel + function of first order, \exp(-|x|) i_1(x). + + -- Function: double gsl_sf_bessel_i2_scaled (double X) + -- Function: int gsl_sf_bessel_i2_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled regular modified spherical Bessel + function of second order, \exp(-|x|) i_2(x) + + -- Function: double gsl_sf_bessel_il_scaled (int L, double X) + -- Function: int gsl_sf_bessel_il_scaled_e (int L, double X, + gsl_sf_result * RESULT) + These routines compute the scaled regular modified spherical Bessel + function of order L, \exp(-|x|) i_l(x) + + -- Function: int gsl_sf_bessel_il_scaled_array (int LMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the scaled regular modified + cylindrical Bessel functions \exp(-|x|) i_l(x) for l from 0 to + LMAX inclusive for lmax >= 0, storing the results in the array + RESULT_ARRAY. The values are computed using recurrence relations + for efficiency, and therefore may differ slightly from the exact + values. + + +File: gsl-ref.info, Node: Irregular Modified Spherical Bessel Functions, Next: Regular Bessel Function - Fractional Order, Prev: Regular Modified Spherical Bessel Functions, Up: Bessel Functions + +7.5.8 Irregular Modified Spherical Bessel Functions +--------------------------------------------------- + +The irregular modified spherical Bessel functions k_l(x) are related to +the irregular modified Bessel functions of fractional order, k_l(x) = +\sqrt{\pi/(2x)} K_{l+1/2}(x). + + -- Function: double gsl_sf_bessel_k0_scaled (double X) + -- Function: int gsl_sf_bessel_k0_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled irregular modified spherical + Bessel function of zeroth order, \exp(x) k_0(x), for x>0. + + -- Function: double gsl_sf_bessel_k1_scaled (double X) + -- Function: int gsl_sf_bessel_k1_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled irregular modified spherical + Bessel function of first order, \exp(x) k_1(x), for x>0. + + -- Function: double gsl_sf_bessel_k2_scaled (double X) + -- Function: int gsl_sf_bessel_k2_scaled_e (double X, gsl_sf_result * + RESULT) + These routines compute the scaled irregular modified spherical + Bessel function of second order, \exp(x) k_2(x), for x>0. + + -- Function: double gsl_sf_bessel_kl_scaled (int L, double X) + -- Function: int gsl_sf_bessel_kl_scaled_e (int L, double X, + gsl_sf_result * RESULT) + These routines compute the scaled irregular modified spherical + Bessel function of order L, \exp(x) k_l(x), for x>0. + + -- Function: int gsl_sf_bessel_kl_scaled_array (int LMAX, double X, + double RESULT_ARRAY[]) + This routine computes the values of the scaled irregular modified + spherical Bessel functions \exp(x) k_l(x) for l from 0 to LMAX + inclusive for lmax >= 0 and x>0, storing the results in the array + RESULT_ARRAY. The values are computed using recurrence relations + for efficiency, and therefore may differ slightly from the exact + values. + + +File: gsl-ref.info, Node: Regular Bessel Function - Fractional Order, Next: Irregular Bessel Functions - Fractional Order, Prev: Irregular Modified Spherical Bessel Functions, Up: Bessel Functions + +7.5.9 Regular Bessel Function--Fractional Order +----------------------------------------------- + + -- Function: double gsl_sf_bessel_Jnu (double NU, double X) + -- Function: int gsl_sf_bessel_Jnu_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the regular cylindrical Bessel function of + fractional order \nu, J_\nu(x). + + -- Function: int gsl_sf_bessel_sequence_Jnu_e (double NU, gsl_mode_t + MODE, size_t SIZE, double V[]) + This function computes the regular cylindrical Bessel function of + fractional order \nu, J_\nu(x), evaluated at a series of x values. + The array V of length SIZE contains the x values. They are + assumed to be strictly ordered and positive. The array is + over-written with the values of J_\nu(x_i). + + +File: gsl-ref.info, Node: Irregular Bessel Functions - Fractional Order, Next: Regular Modified Bessel Functions - Fractional Order, Prev: Regular Bessel Function - Fractional Order, Up: Bessel Functions + +7.5.10 Irregular Bessel Functions--Fractional Order +--------------------------------------------------- + + -- Function: double gsl_sf_bessel_Ynu (double NU, double X) + -- Function: int gsl_sf_bessel_Ynu_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the irregular cylindrical Bessel function of + fractional order \nu, Y_\nu(x). + + +File: gsl-ref.info, Node: Regular Modified Bessel Functions - Fractional Order, Next: Irregular Modified Bessel Functions - Fractional Order, Prev: Irregular Bessel Functions - Fractional Order, Up: Bessel Functions + +7.5.11 Regular Modified Bessel Functions--Fractional Order +---------------------------------------------------------- + + -- Function: double gsl_sf_bessel_Inu (double NU, double X) + -- Function: int gsl_sf_bessel_Inu_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the regular modified Bessel function of + fractional order \nu, I_\nu(x) for x>0, \nu>0. + + -- Function: double gsl_sf_bessel_Inu_scaled (double NU, double X) + -- Function: int gsl_sf_bessel_Inu_scaled_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the scaled regular modified Bessel function + of fractional order \nu, \exp(-|x|)I_\nu(x) for x>0, \nu>0. + + +File: gsl-ref.info, Node: Irregular Modified Bessel Functions - Fractional Order, Next: Zeros of Regular Bessel Functions, Prev: Regular Modified Bessel Functions - Fractional Order, Up: Bessel Functions + +7.5.12 Irregular Modified Bessel Functions--Fractional Order +------------------------------------------------------------ + + -- Function: double gsl_sf_bessel_Knu (double NU, double X) + -- Function: int gsl_sf_bessel_Knu_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the irregular modified Bessel function of + fractional order \nu, K_\nu(x) for x>0, \nu>0. + + -- Function: double gsl_sf_bessel_lnKnu (double NU, double X) + -- Function: int gsl_sf_bessel_lnKnu_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the logarithm of the irregular modified + Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0, + \nu>0. + + -- Function: double gsl_sf_bessel_Knu_scaled (double NU, double X) + -- Function: int gsl_sf_bessel_Knu_scaled_e (double NU, double X, + gsl_sf_result * RESULT) + These routines compute the scaled irregular modified Bessel + function of fractional order \nu, \exp(+|x|) K_\nu(x) for x>0, + \nu>0. + + +File: gsl-ref.info, Node: Zeros of Regular Bessel Functions, Prev: Irregular Modified Bessel Functions - Fractional Order, Up: Bessel Functions + +7.5.13 Zeros of Regular Bessel Functions +---------------------------------------- + + -- Function: double gsl_sf_bessel_zero_J0 (unsigned int S) + -- Function: int gsl_sf_bessel_zero_J0_e (unsigned int S, + gsl_sf_result * RESULT) + These routines compute the location of the S-th positive zero of + the Bessel function J_0(x). + + -- Function: double gsl_sf_bessel_zero_J1 (unsigned int S) + -- Function: int gsl_sf_bessel_zero_J1_e (unsigned int S, + gsl_sf_result * RESULT) + These routines compute the location of the S-th positive zero of + the Bessel function J_1(x). + + -- Function: double gsl_sf_bessel_zero_Jnu (double NU, unsigned int S) + -- Function: int gsl_sf_bessel_zero_Jnu_e (double NU, unsigned int S, + gsl_sf_result * RESULT) + These routines compute the location of the S-th positive zero of + the Bessel function J_\nu(x). The current implementation does not + support negative values of NU. + + +File: gsl-ref.info, Node: Clausen Functions, Next: Coulomb Functions, Prev: Bessel Functions, Up: Special Functions + +7.6 Clausen Functions +===================== + +The Clausen function is defined by the following integral, + + Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2)) + +It is related to the dilogarithm by Cl_2(\theta) = \Im +Li_2(\exp(i\theta)). The Clausen functions are declared in the header +file `gsl_sf_clausen.h'. + + -- Function: double gsl_sf_clausen (double X) + -- Function: int gsl_sf_clausen_e (double X, gsl_sf_result * RESULT) + These routines compute the Clausen integral Cl_2(x). + + +File: gsl-ref.info, Node: Coulomb Functions, Next: Coupling Coefficients, Prev: Clausen Functions, Up: Special Functions + +7.7 Coulomb Functions +===================== + +The prototypes of the Coulomb functions are declared in the header file +`gsl_sf_coulomb.h'. Both bound state and scattering solutions are +available. + +* Menu: + +* Normalized Hydrogenic Bound States:: +* Coulomb Wave Functions:: +* Coulomb Wave Function Normalization Constant:: + + +File: gsl-ref.info, Node: Normalized Hydrogenic Bound States, Next: Coulomb Wave Functions, Up: Coulomb Functions + +7.7.1 Normalized Hydrogenic Bound States +---------------------------------------- + + -- Function: double gsl_sf_hydrogenicR_1 (double Z, double R) + -- Function: int gsl_sf_hydrogenicR_1_e (double Z, double R, + gsl_sf_result * RESULT) + These routines compute the lowest-order normalized hydrogenic bound + state radial wavefunction R_1 := 2Z \sqrt{Z} \exp(-Z r). + + -- Function: double gsl_sf_hydrogenicR (int N, int L, double Z, double + R) + -- Function: int gsl_sf_hydrogenicR_e (int N, int L, double Z, double + R, gsl_sf_result * RESULT) + These routines compute the N-th normalized hydrogenic bound state + radial wavefunction, + + R_n := 2 (Z^{3/2}/n^2) \sqrt{(n-l-1)!/(n+l)!} \exp(-Z r/n) (2Zr/n)^l + L^{2l+1}_{n-l-1}(2Zr/n). + + where L^a_b(x) is the generalized Laguerre polynomial (*note + Laguerre Functions::). The normalization is chosen such that the + wavefunction \psi is given by \psi(n,l,r) = R_n Y_{lm}. + + +File: gsl-ref.info, Node: Coulomb Wave Functions, Next: Coulomb Wave Function Normalization Constant, Prev: Normalized Hydrogenic Bound States, Up: Coulomb Functions + +7.7.2 Coulomb Wave Functions +---------------------------- + +The Coulomb wave functions F_L(\eta,x), G_L(\eta,x) are described in +Abramowitz & Stegun, Chapter 14. Because there can be a large dynamic +range of values for these functions, overflows are handled gracefully. +If an overflow occurs, `GSL_EOVRFLW' is signalled and exponent(s) are +returned through the modifiable parameters EXP_F, EXP_G. The full +solution can be reconstructed from the following relations, + + F_L(eta,x) = fc[k_L] * exp(exp_F) + G_L(eta,x) = gc[k_L] * exp(exp_G) + + F_L'(eta,x) = fcp[k_L] * exp(exp_F) + G_L'(eta,x) = gcp[k_L] * exp(exp_G) + + + -- Function: int gsl_sf_coulomb_wave_FG_e (double ETA, double X, + double L_F, int K, gsl_sf_result * F, gsl_sf_result * FP, + gsl_sf_result * G, gsl_sf_result * GP, double * EXP_F, double + * EXP_G) + This function computes the Coulomb wave functions F_L(\eta,x), + G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x), + G'_{L-k}(\eta,x) with respect to x. The parameters are restricted + to L, L-k > -1/2, x > 0 and integer k. Note that L itself is not + restricted to being an integer. The results are stored in the + parameters F, G for the function values and FP, GP for the + derivative values. If an overflow occurs, `GSL_EOVRFLW' is + returned and scaling exponents are stored in the modifiable + parameters EXP_F, EXP_G. + + -- Function: int gsl_sf_coulomb_wave_F_array (double L_MIN, int KMAX, + double ETA, double X, double FC_ARRAY[], double * F_EXPONENT) + This function computes the Coulomb wave function F_L(\eta,x) for L + = Lmin \dots Lmin + kmax, storing the results in FC_ARRAY. In the + case of overflow the exponent is stored in F_EXPONENT. + + -- Function: int gsl_sf_coulomb_wave_FG_array (double L_MIN, int KMAX, + double ETA, double X, double FC_ARRAY[], double GC_ARRAY[], + double * F_EXPONENT, double * G_EXPONENT) + This function computes the functions F_L(\eta,x), G_L(\eta,x) for + L = Lmin \dots Lmin + kmax storing the results in FC_ARRAY and + GC_ARRAY. In the case of overflow the exponents are stored in + F_EXPONENT and G_EXPONENT. + + -- Function: int gsl_sf_coulomb_wave_FGp_array (double L_MIN, int + KMAX, double ETA, double X, double FC_ARRAY[], double + FCP_ARRAY[], double GC_ARRAY[], double GCP_ARRAY[], double * + F_EXPONENT, double * G_EXPONENT) + This function computes the functions F_L(\eta,x), G_L(\eta,x) and + their derivatives F'_L(\eta,x), G'_L(\eta,x) for L = Lmin \dots + Lmin + kmax storing the results in FC_ARRAY, GC_ARRAY, FCP_ARRAY + and GCP_ARRAY. In the case of overflow the exponents are stored + in F_EXPONENT and G_EXPONENT. + + -- Function: int gsl_sf_coulomb_wave_sphF_array (double L_MIN, int + KMAX, double ETA, double X, double FC_ARRAY[], double + F_EXPONENT[]) + This function computes the Coulomb wave function divided by the + argument F_L(\eta, x)/x for L = Lmin \dots Lmin + kmax, storing the + results in FC_ARRAY. In the case of overflow the exponent is + stored in F_EXPONENT. This function reduces to spherical Bessel + functions in the limit \eta \to 0. + + +File: gsl-ref.info, Node: Coulomb Wave Function Normalization Constant, Prev: Coulomb Wave Functions, Up: Coulomb Functions + +7.7.3 Coulomb Wave Function Normalization Constant +-------------------------------------------------- + +The Coulomb wave function normalization constant is defined in +Abramowitz 14.1.7. + + -- Function: int gsl_sf_coulomb_CL_e (double L, double ETA, + gsl_sf_result * RESULT) + This function computes the Coulomb wave function normalization + constant C_L(\eta) for L > -1. + + -- Function: int gsl_sf_coulomb_CL_array (double LMIN, int KMAX, + double ETA, double CL[]) + This function computes the Coulomb wave function normalization + constant C_L(\eta) for L = Lmin \dots Lmin + kmax, Lmin > -1. + + +File: gsl-ref.info, Node: Coupling Coefficients, Next: Dawson Function, Prev: Coulomb Functions, Up: Special Functions + +7.8 Coupling Coefficients +========================= + +The Wigner 3-j, 6-j and 9-j symbols give the coupling coefficients for +combined angular momentum vectors. Since the arguments of the standard +coupling coefficient functions are integer or half-integer, the +arguments of the following functions are, by convention, integers equal +to twice the actual spin value. For information on the 3-j coefficients +see Abramowitz & Stegun, Section 27.9. The functions described in this +section are declared in the header file `gsl_sf_coupling.h'. + +* Menu: + +* 3-j Symbols:: +* 6-j Symbols:: +* 9-j Symbols:: + + +File: gsl-ref.info, Node: 3-j Symbols, Next: 6-j Symbols, Up: Coupling Coefficients + +7.8.1 3-j Symbols +----------------- + + -- Function: double gsl_sf_coupling_3j (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC) + -- Function: int gsl_sf_coupling_3j_e (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC, gsl_sf_result * + RESULT) + These routines compute the Wigner 3-j coefficient, + + (ja jb jc + ma mb mc) + + where the arguments are given in half-integer units, ja = + TWO_JA/2, ma = TWO_MA/2, etc. + + +File: gsl-ref.info, Node: 6-j Symbols, Next: 9-j Symbols, Prev: 3-j Symbols, Up: Coupling Coefficients + +7.8.2 6-j Symbols +----------------- + + -- Function: double gsl_sf_coupling_6j (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF) + -- Function: int gsl_sf_coupling_6j_e (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, gsl_sf_result * + RESULT) + These routines compute the Wigner 6-j coefficient, + + {ja jb jc + jd je jf} + + where the arguments are given in half-integer units, ja = + TWO_JA/2, ma = TWO_MA/2, etc. + + +File: gsl-ref.info, Node: 9-j Symbols, Prev: 6-j Symbols, Up: Coupling Coefficients + +7.8.3 9-j Symbols +----------------- + + -- Function: double gsl_sf_coupling_9j (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int + TWO_JH, int TWO_JI) + -- Function: int gsl_sf_coupling_9j_e (int TWO_JA, int TWO_JB, int + TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int + TWO_JH, int TWO_JI, gsl_sf_result * RESULT) + These routines compute the Wigner 9-j coefficient, + + {ja jb jc + jd je jf + jg jh ji} + + where the arguments are given in half-integer units, ja = + TWO_JA/2, ma = TWO_MA/2, etc. + + +File: gsl-ref.info, Node: Dawson Function, Next: Debye Functions, Prev: Coupling Coefficients, Up: Special Functions + +7.9 Dawson Function +=================== + +The Dawson integral is defined by \exp(-x^2) \int_0^x dt \exp(t^2). A +table of Dawson's integral can be found in Abramowitz & Stegun, Table +7.5. The Dawson functions are declared in the header file +`gsl_sf_dawson.h'. + + -- Function: double gsl_sf_dawson (double X) + -- Function: int gsl_sf_dawson_e (double X, gsl_sf_result * RESULT) + These routines compute the value of Dawson's integral for X. + + +File: gsl-ref.info, Node: Debye Functions, Next: Dilogarithm, Prev: Dawson Function, Up: Special Functions + +7.10 Debye Functions +==================== + +The Debye functions D_n(x) are defined by the following integral, + + D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1)) + +For further information see Abramowitz & Stegun, Section 27.1. The +Debye functions are declared in the header file `gsl_sf_debye.h'. + + -- Function: double gsl_sf_debye_1 (double X) + -- Function: int gsl_sf_debye_1_e (double X, gsl_sf_result * RESULT) + These routines compute the first-order Debye function D_1(x) = + (1/x) \int_0^x dt (t/(e^t - 1)). + + -- Function: double gsl_sf_debye_2 (double X) + -- Function: int gsl_sf_debye_2_e (double X, gsl_sf_result * RESULT) + These routines compute the second-order Debye function D_2(x) = + (2/x^2) \int_0^x dt (t^2/(e^t - 1)). + + -- Function: double gsl_sf_debye_3 (double X) + -- Function: int gsl_sf_debye_3_e (double X, gsl_sf_result * RESULT) + These routines compute the third-order Debye function D_3(x) = + (3/x^3) \int_0^x dt (t^3/(e^t - 1)). + + -- Function: double gsl_sf_debye_4 (double X) + -- Function: int gsl_sf_debye_4_e (double X, gsl_sf_result * RESULT) + These routines compute the fourth-order Debye function D_4(x) = + (4/x^4) \int_0^x dt (t^4/(e^t - 1)). + + -- Function: double gsl_sf_debye_5 (double X) + -- Function: int gsl_sf_debye_5_e (double X, gsl_sf_result * RESULT) + These routines compute the fifth-order Debye function D_5(x) = + (5/x^5) \int_0^x dt (t^5/(e^t - 1)). + + -- Function: double gsl_sf_debye_6 (double X) + -- Function: int gsl_sf_debye_6_e (double X, gsl_sf_result * RESULT) + These routines compute the sixth-order Debye function D_6(x) = + (6/x^6) \int_0^x dt (t^6/(e^t - 1)). + + +File: gsl-ref.info, Node: Dilogarithm, Next: Elementary Operations, Prev: Debye Functions, Up: Special Functions + +7.11 Dilogarithm +================ + +The functions described in this section are declared in the header file +`gsl_sf_dilog.h'. + +* Menu: + +* Real Argument:: +* Complex Argument:: + + +File: gsl-ref.info, Node: Real Argument, Next: Complex Argument, Up: Dilogarithm + +7.11.1 Real Argument +-------------------- + + -- Function: double gsl_sf_dilog (double X) + -- Function: int gsl_sf_dilog_e (double X, gsl_sf_result * RESULT) + These routines compute the dilogarithm for a real argument. In + Lewin's notation this is Li_2(x), the real part of the dilogarithm + of a real x. It is defined by the integral representation Li_2(x) + = - \Re \int_0^x ds \log(1-s) / s. Note that \Im(Li_2(x)) = 0 for + x <= 1, and -\pi\log(x) for x > 1. + + + +File: gsl-ref.info, Node: Complex Argument, Prev: Real Argument, Up: Dilogarithm + +7.11.2 Complex Argument +----------------------- + + -- Function: int gsl_sf_complex_dilog_e (double R, double THETA, + gsl_sf_result * RESULT_RE, gsl_sf_result * RESULT_IM) + This function computes the full complex-valued dilogarithm for the + complex argument z = r \exp(i \theta). The real and imaginary + parts of the result are returned in RESULT_RE, RESULT_IM. + + +File: gsl-ref.info, Node: Elementary Operations, Next: Elliptic Integrals, Prev: Dilogarithm, Up: Special Functions + +7.12 Elementary Operations +========================== + +The following functions allow for the propagation of errors when +combining quantities by multiplication. The functions are declared in +the header file `gsl_sf_elementary.h'. + + -- Function: int gsl_sf_multiply_e (double X, double Y, gsl_sf_result + * RESULT) + This function multiplies X and Y storing the product and its + associated error in RESULT. + + -- Function: int gsl_sf_multiply_err_e (double X, double DX, double Y, + double DY, gsl_sf_result * RESULT) + This function multiplies X and Y with associated absolute errors + DX and DY. The product xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2) is + stored in RESULT. + + +File: gsl-ref.info, Node: Elliptic Integrals, Next: Elliptic Functions (Jacobi), Prev: Elementary Operations, Up: Special Functions + +7.13 Elliptic Integrals +======================= + +The functions described in this section are declared in the header file +`gsl_sf_ellint.h'. Further information about the elliptic integrals +can be found in Abramowitz & Stegun, Chapter 17. + +* Menu: + +* Definition of Legendre Forms:: +* Definition of Carlson Forms:: +* Legendre Form of Complete Elliptic Integrals:: +* Legendre Form of Incomplete Elliptic Integrals:: +* Carlson Forms:: + + +File: gsl-ref.info, Node: Definition of Legendre Forms, Next: Definition of Carlson Forms, Up: Elliptic Integrals + +7.13.1 Definition of Legendre Forms +----------------------------------- + +The Legendre forms of elliptic integrals F(\phi,k), E(\phi,k) and +\Pi(\phi,k,n) are defined by, + + F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t))) + + E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t))) + + Pi(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t))) + +The complete Legendre forms are denoted by K(k) = F(\pi/2, k) and E(k) += E(\pi/2, k). + + The notation used here is based on Carlson, `Numerische Mathematik' +33 (1979) 1 and differs slightly from that used by Abramowitz & Stegun, +where the functions are given in terms of the parameter m = k^2 and n +is replaced by -n. + + +File: gsl-ref.info, Node: Definition of Carlson Forms, Next: Legendre Form of Complete Elliptic Integrals, Prev: Definition of Legendre Forms, Up: Elliptic Integrals + +7.13.2 Definition of Carlson Forms +---------------------------------- + +The Carlson symmetric forms of elliptical integrals RC(x,y), RD(x,y,z), +RF(x,y,z) and RJ(x,y,z,p) are defined by, + + RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1) + + RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2) + + RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) + + RJ(x,y,z,p) = 3/2 \int_0^\infty dt + (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1) + + +File: gsl-ref.info, Node: Legendre Form of Complete Elliptic Integrals, Next: Legendre Form of Incomplete Elliptic Integrals, Prev: Definition of Carlson Forms, Up: Elliptic Integrals + +7.13.3 Legendre Form of Complete Elliptic Integrals +--------------------------------------------------- + + -- Function: double gsl_sf_ellint_Kcomp (double K, gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_Kcomp_e (double K, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the complete elliptic integral K(k) to the + accuracy specified by the mode variable MODE. Note that + Abramowitz & Stegun define this function in terms of the parameter + m = k^2. + + -- Function: double gsl_sf_ellint_Ecomp (double K, gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_Ecomp_e (double K, gsl_mode_t MODE, + gsl_sf_result * RESULT) + These routines compute the complete elliptic integral E(k) to the + accuracy specified by the mode variable MODE. Note that + Abramowitz & Stegun define this function in terms of the parameter + m = k^2. + + -- Function: double gsl_sf_ellint_Pcomp (double K, double N, + gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_Pcomp_e (double K, double N, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the complete elliptic integral \Pi(k,n) to + the accuracy specified by the mode variable MODE. Note that + Abramowitz & Stegun define this function in terms of the + parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of + sign n \to -n. + + +File: gsl-ref.info, Node: Legendre Form of Incomplete Elliptic Integrals, Next: Carlson Forms, Prev: Legendre Form of Complete Elliptic Integrals, Up: Elliptic Integrals + +7.13.4 Legendre Form of Incomplete Elliptic Integrals +----------------------------------------------------- + + -- Function: double gsl_sf_ellint_F (double PHI, double K, gsl_mode_t + MODE) + -- Function: int gsl_sf_ellint_F_e (double PHI, double K, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral F(\phi,k) + to the accuracy specified by the mode variable MODE. Note that + Abramowitz & Stegun define this function in terms of the parameter + m = k^2. + + -- Function: double gsl_sf_ellint_E (double PHI, double K, gsl_mode_t + MODE) + -- Function: int gsl_sf_ellint_E_e (double PHI, double K, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral E(\phi,k) + to the accuracy specified by the mode variable MODE. Note that + Abramowitz & Stegun define this function in terms of the parameter + m = k^2. + + -- Function: double gsl_sf_ellint_P (double PHI, double K, double N, + gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_P_e (double PHI, double K, double N, + gsl_mode_t MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral + \Pi(\phi,k,n) to the accuracy specified by the mode variable MODE. + Note that Abramowitz & Stegun define this function in terms of the + parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of + sign n \to -n. + + -- Function: double gsl_sf_ellint_D (double PHI, double K, double N, + gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_D_e (double PHI, double K, double N, + gsl_mode_t MODE, gsl_sf_result * RESULT) + These functions compute the incomplete elliptic integral D(\phi,k) + which is defined through the Carlson form RD(x,y,z) by the + following relation, + + D(\phi,k,n) = (1/3)(\sin(\phi))^3 RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1). + The argument N is not used and will be removed in a future release. + + + +File: gsl-ref.info, Node: Carlson Forms, Prev: Legendre Form of Incomplete Elliptic Integrals, Up: Elliptic Integrals + +7.13.5 Carlson Forms +-------------------- + + -- Function: double gsl_sf_ellint_RC (double X, double Y, gsl_mode_t + MODE) + -- Function: int gsl_sf_ellint_RC_e (double X, double Y, gsl_mode_t + MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral RC(x,y) to + the accuracy specified by the mode variable MODE. + + -- Function: double gsl_sf_ellint_RD (double X, double Y, double Z, + gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_RD_e (double X, double Y, double Z, + gsl_mode_t MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral RD(x,y,z) + to the accuracy specified by the mode variable MODE. + + -- Function: double gsl_sf_ellint_RF (double X, double Y, double Z, + gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_RF_e (double X, double Y, double Z, + gsl_mode_t MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral RF(x,y,z) + to the accuracy specified by the mode variable MODE. + + -- Function: double gsl_sf_ellint_RJ (double X, double Y, double Z, + double P, gsl_mode_t MODE) + -- Function: int gsl_sf_ellint_RJ_e (double X, double Y, double Z, + double P, gsl_mode_t MODE, gsl_sf_result * RESULT) + These routines compute the incomplete elliptic integral RJ(x,y,z,p) + to the accuracy specified by the mode variable MODE. + + +File: gsl-ref.info, Node: Elliptic Functions (Jacobi), Next: Error Functions, Prev: Elliptic Integrals, Up: Special Functions + +7.14 Elliptic Functions (Jacobi) +================================ + +The Jacobian Elliptic functions are defined in Abramowitz & Stegun, +Chapter 16. The functions are declared in the header file +`gsl_sf_elljac.h'. + + -- Function: int gsl_sf_elljac_e (double U, double M, double * SN, + double * CN, double * DN) + This function computes the Jacobian elliptic functions sn(u|m), + cn(u|m), dn(u|m) by descending Landen transformations. + + +File: gsl-ref.info, Node: Error Functions, Next: Exponential Functions, Prev: Elliptic Functions (Jacobi), Up: Special Functions + +7.15 Error Functions +==================== + +The error function is described in Abramowitz & Stegun, Chapter 7. The +functions in this section are declared in the header file +`gsl_sf_erf.h'. + +* Menu: + +* Error Function:: +* Complementary Error Function:: +* Log Complementary Error Function:: +* Probability functions:: + + +File: gsl-ref.info, Node: Error Function, Next: Complementary Error Function, Up: Error Functions + +7.15.1 Error Function +--------------------- + + -- Function: double gsl_sf_erf (double X) + -- Function: int gsl_sf_erf_e (double X, gsl_sf_result * RESULT) + These routines compute the error function erf(x), where erf(x) = + (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2). + + +File: gsl-ref.info, Node: Complementary Error Function, Next: Log Complementary Error Function, Prev: Error Function, Up: Error Functions + +7.15.2 Complementary Error Function +----------------------------------- + + -- Function: double gsl_sf_erfc (double X) + -- Function: int gsl_sf_erfc_e (double X, gsl_sf_result * RESULT) + These routines compute the complementary error function erfc(x) = + 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2). + + +File: gsl-ref.info, Node: Log Complementary Error Function, Next: Probability functions, Prev: Complementary Error Function, Up: Error Functions + +7.15.3 Log Complementary Error Function +--------------------------------------- + + -- Function: double gsl_sf_log_erfc (double X) + -- Function: int gsl_sf_log_erfc_e (double X, gsl_sf_result * RESULT) + These routines compute the logarithm of the complementary error + function \log(\erfc(x)). + + +File: gsl-ref.info, Node: Probability functions, Prev: Log Complementary Error Function, Up: Error Functions + +7.15.4 Probability functions +---------------------------- + +The probability functions for the Normal or Gaussian distribution are +described in Abramowitz & Stegun, Section 26.2. + + -- Function: double gsl_sf_erf_Z (double X) + -- Function: int gsl_sf_erf_Z_e (double X, gsl_sf_result * RESULT) + These routines compute the Gaussian probability density function + Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2). + + -- Function: double gsl_sf_erf_Q (double X) + -- Function: int gsl_sf_erf_Q_e (double X, gsl_sf_result * RESULT) + These routines compute the upper tail of the Gaussian probability + function Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2). + + The "hazard function" for the normal distribution, also known as the +inverse Mill's ratio, is defined as, + + h(x) = Z(x)/Q(x) = \sqrt{2/\pi} \exp(-x^2 / 2) / \erfc(x/\sqrt 2) + +It decreases rapidly as x approaches -\infty and asymptotes to h(x) +\sim x as x approaches +\infty. + + -- Function: double gsl_sf_hazard (double X) + -- Function: int gsl_sf_hazard_e (double X, gsl_sf_result * RESULT) + These routines compute the hazard function for the normal + distribution. + + +File: gsl-ref.info, Node: Exponential Functions, Next: Exponential Integrals, Prev: Error Functions, Up: Special Functions + +7.16 Exponential Functions +========================== + +The functions described in this section are declared in the header file +`gsl_sf_exp.h'. + +* Menu: + +* Exponential Function:: +* Relative Exponential Functions:: +* Exponentiation With Error Estimate:: + + +File: gsl-ref.info, Node: Exponential Function, Next: Relative Exponential Functions, Up: Exponential Functions + +7.16.1 Exponential Function +--------------------------- + + -- Function: double gsl_sf_exp (double X) + -- Function: int gsl_sf_exp_e (double X, gsl_sf_result * RESULT) + These routines provide an exponential function \exp(x) using GSL + semantics and error checking. + + -- Function: int gsl_sf_exp_e10_e (double X, gsl_sf_result_e10 * + RESULT) + This function computes the exponential \exp(x) using the + `gsl_sf_result_e10' type to return a result with extended range. + This function may be useful if the value of \exp(x) would overflow + the numeric range of `double'. + + -- Function: double gsl_sf_exp_mult (double X, double Y) + -- Function: int gsl_sf_exp_mult_e (double X, double Y, gsl_sf_result + * RESULT) + These routines exponentiate X and multiply by the factor Y to + return the product y \exp(x). + + -- Function: int gsl_sf_exp_mult_e10_e (const double X, const double + Y, gsl_sf_result_e10 * RESULT) + This function computes the product y \exp(x) using the + `gsl_sf_result_e10' type to return a result with extended numeric + range. + + +File: gsl-ref.info, Node: Relative Exponential Functions, Next: Exponentiation With Error Estimate, Prev: Exponential Function, Up: Exponential Functions + +7.16.2 Relative Exponential Functions +------------------------------------- + + -- Function: double gsl_sf_expm1 (double X) + -- Function: int gsl_sf_expm1_e (double X, gsl_sf_result * RESULT) + These routines compute the quantity \exp(x)-1 using an algorithm + that is accurate for small x. + + -- Function: double gsl_sf_exprel (double X) + -- Function: int gsl_sf_exprel_e (double X, gsl_sf_result * RESULT) + These routines compute the quantity (\exp(x)-1)/x using an + algorithm that is accurate for small x. For small x the algorithm + is based on the expansion (\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + + x^3/(2*3*4) + \dots. + + -- Function: double gsl_sf_exprel_2 (double X) + -- Function: int gsl_sf_exprel_2_e (double X, gsl_sf_result * RESULT) + These routines compute the quantity 2(\exp(x)-1-x)/x^2 using an + algorithm that is accurate for small x. For small x the algorithm + is based on the expansion 2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + + x^3/(3*4*5) + \dots. + + -- Function: double gsl_sf_exprel_n (int N, double X) + -- Function: int gsl_sf_exprel_n_e (int N, double X, gsl_sf_result * + RESULT) + These routines compute the N-relative exponential, which is the + N-th generalization of the functions `gsl_sf_exprel' and + `gsl_sf_exprel2'. The N-relative exponential is given by, + + exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!) + = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ... + = 1F1 (1,1+N,x) + + +File: gsl-ref.info, Node: Exponentiation With Error Estimate, Prev: Relative Exponential Functions, Up: Exponential Functions + +7.16.3 Exponentiation With Error Estimate +----------------------------------------- + + -- Function: int gsl_sf_exp_err_e (double X, double DX, gsl_sf_result + * RESULT) + This function exponentiates X with an associated absolute error DX. + + -- Function: int gsl_sf_exp_err_e10_e (double X, double DX, + gsl_sf_result_e10 * RESULT) + This function exponentiates a quantity X with an associated + absolute error DX using the `gsl_sf_result_e10' type to return a + result with extended range. + + -- Function: int gsl_sf_exp_mult_err_e (double X, double DX, double Y, + double DY, gsl_sf_result * RESULT) + This routine computes the product y \exp(x) for the quantities X, + Y with associated absolute errors DX, DY. + + -- Function: int gsl_sf_exp_mult_err_e10_e (double X, double DX, + double Y, double DY, gsl_sf_result_e10 * RESULT) + This routine computes the product y \exp(x) for the quantities X, + Y with associated absolute errors DX, DY using the + `gsl_sf_result_e10' type to return a result with extended range. + + +File: gsl-ref.info, Node: Exponential Integrals, Next: Fermi-Dirac Function, Prev: Exponential Functions, Up: Special Functions + +7.17 Exponential Integrals +========================== + +Information on the exponential integrals can be found in Abramowitz & +Stegun, Chapter 5. These functions are declared in the header file +`gsl_sf_expint.h'. + +* Menu: + +* Exponential Integral:: +* Ei(x):: +* Hyperbolic Integrals:: +* Ei_3(x):: +* Trigonometric Integrals:: +* Arctangent Integral:: + + +File: gsl-ref.info, Node: Exponential Integral, Next: Ei(x), Up: Exponential Integrals + +7.17.1 Exponential Integral +--------------------------- + + -- Function: double gsl_sf_expint_E1 (double X) + -- Function: int gsl_sf_expint_E1_e (double X, gsl_sf_result * RESULT) + These routines compute the exponential integral E_1(x), + + E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t. + + + + -- Function: double gsl_sf_expint_E2 (double X) + -- Function: int gsl_sf_expint_E2_e (double X, gsl_sf_result * RESULT) + These routines compute the second-order exponential integral + E_2(x), + + E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2. + + + + +File: gsl-ref.info, Node: Ei(x), Next: Hyperbolic Integrals, Prev: Exponential Integral, Up: Exponential Integrals + +7.17.2 Ei(x) +------------ + + -- Function: double gsl_sf_expint_Ei (double X) + -- Function: int gsl_sf_expint_Ei_e (double X, gsl_sf_result * RESULT) + These routines compute the exponential integral Ei(x), + + Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t) + + where PV denotes the principal value of the integral. + + +File: gsl-ref.info, Node: Hyperbolic Integrals, Next: Ei_3(x), Prev: Ei(x), Up: Exponential Integrals + +7.17.3 Hyperbolic Integrals +--------------------------- + + -- Function: double gsl_sf_Shi (double X) + -- Function: int gsl_sf_Shi_e (double X, gsl_sf_result * RESULT) + These routines compute the integral Shi(x) = \int_0^x dt + \sinh(t)/t. + + -- Function: double gsl_sf_Chi (double X) + -- Function: int gsl_sf_Chi_e (double X, gsl_sf_result * RESULT) + These routines compute the integral Chi(x) := \Re[ \gamma_E + + \log(x) + \int_0^x dt (\cosh[t]-1)/t] , where \gamma_E is the + Euler constant (available as the macro `M_EULER'). + + +File: gsl-ref.info, Node: Ei_3(x), Next: Trigonometric Integrals, Prev: Hyperbolic Integrals, Up: Exponential Integrals + +7.17.4 Ei_3(x) +-------------- + + -- Function: double gsl_sf_expint_3 (double X) + -- Function: int gsl_sf_expint_3_e (double X, gsl_sf_result * RESULT) + These routines compute the third-order exponential integral + Ei_3(x) = \int_0^xdt \exp(-t^3) for x >= 0. + + +File: gsl-ref.info, Node: Trigonometric Integrals, Next: Arctangent Integral, Prev: Ei_3(x), Up: Exponential Integrals + +7.17.5 Trigonometric Integrals +------------------------------ + + -- Function: double gsl_sf_Si (const double X) + -- Function: int gsl_sf_Si_e (double X, gsl_sf_result * RESULT) + These routines compute the Sine integral Si(x) = \int_0^x dt + \sin(t)/t. + + -- Function: double gsl_sf_Ci (const double X) + -- Function: int gsl_sf_Ci_e (double X, gsl_sf_result * RESULT) + These routines compute the Cosine integral Ci(x) = -\int_x^\infty + dt \cos(t)/t for x > 0. + + +File: gsl-ref.info, Node: Arctangent Integral, Prev: Trigonometric Integrals, Up: Exponential Integrals + +7.17.6 Arctangent Integral +-------------------------- + + -- Function: double gsl_sf_atanint (double X) + -- Function: int gsl_sf_atanint_e (double X, gsl_sf_result * RESULT) + These routines compute the Arctangent integral, which is defined + as AtanInt(x) = \int_0^x dt \arctan(t)/t. + + +File: gsl-ref.info, Node: Fermi-Dirac Function, Next: Gamma and Beta Functions, Prev: Exponential Integrals, Up: Special Functions + +7.18 Fermi-Dirac Function +========================= + +The functions described in this section are declared in the header file +`gsl_sf_fermi_dirac.h'. + +* Menu: + +* Complete Fermi-Dirac Integrals:: +* Incomplete Fermi-Dirac Integrals:: + + +File: gsl-ref.info, Node: Complete Fermi-Dirac Integrals, Next: Incomplete Fermi-Dirac Integrals, Up: Fermi-Dirac Function + +7.18.1 Complete Fermi-Dirac Integrals +------------------------------------- + +The complete Fermi-Dirac integral F_j(x) is given by, + + F_j(x) := (1/r\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1)) + + -- 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). + + +File: gsl-ref.info, Node: Incomplete Fermi-Dirac Integrals, Prev: Complete Fermi-Dirac Integrals, Up: Fermi-Dirac Function + +7.18.2 Incomplete Fermi-Dirac Integrals +--------------------------------------- + +The incomplete Fermi-Dirac integral F_j(x,b) is given by, + + F_j(x,b) := (1/\Gamma(j+1)) \int_b^\infty dt (t^j / (\Exp(t-x) + 1)) + + -- Function: double gsl_sf_fermi_dirac_inc_0 (double X, double B) + -- Function: int gsl_sf_fermi_dirac_inc_0_e (double X, double B, + gsl_sf_result * RESULT) + These routines compute the incomplete Fermi-Dirac integral with an + index of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x). + + +File: gsl-ref.info, Node: Gamma and Beta Functions, Next: Gegenbauer Functions, Prev: Fermi-Dirac Function, Up: Special Functions + +7.19 Gamma and Beta Functions +============================= + +The functions described in this section are declared in the header +file `gsl_sf_gamma.h'. + +* Menu: + +* Gamma Functions:: +* Factorials:: +* Pochhammer Symbol:: +* Incomplete Gamma Functions:: +* Beta Functions:: +* Incomplete Beta Function:: + + +File: gsl-ref.info, Node: Gamma Functions, Next: Factorials, Up: Gamma and Beta Functions + +7.19.1 Gamma Functions +---------------------- + +The Gamma function is defined by the following integral, + + \Gamma(x) = \int_0^\infty dt t^{x-1} \exp(-t) + +It is related to the factorial function by \Gamma(n)=(n-1)! for +positive integer n. Further information on the Gamma function can be +found in Abramowitz & Stegun, Chapter 6. The functions described in +this section are declared in the header file `gsl_sf_gamma.h'. + + -- Function: double gsl_sf_gamma (double X) + -- Function: int gsl_sf_gamma_e (double X, gsl_sf_result * RESULT) + These routines compute the Gamma function \Gamma(x), subject to x + not being a negative integer or zero. The function is computed + using the real Lanczos method. The maximum value of x such that + \Gamma(x) is not considered an overflow is given by the macro + `GSL_SF_GAMMA_XMAX' and is 171.0. + + -- Function: double gsl_sf_lngamma (double X) + -- Function: int gsl_sf_lngamma_e (double X, gsl_sf_result * RESULT) + These routines compute the logarithm of the Gamma function, + \log(\Gamma(x)), subject to x not being a negative integer or + zero. For x<0 the real part of \log(\Gamma(x)) is returned, which + is equivalent to \log(|\Gamma(x)|). The function is computed + using the real Lanczos method. + + -- Function: int gsl_sf_lngamma_sgn_e (double X, gsl_sf_result * + RESULT_LG, double * SGN) + This routine computes the sign of the gamma function and the + logarithm of its magnitude, subject to x not being a negative + integer or zero. The function is computed using the real Lanczos + method. The value of the gamma function can be reconstructed + using the relation \Gamma(x) = sgn * \exp(resultlg). + + -- Function: double gsl_sf_gammastar (double X) + -- Function: int gsl_sf_gammastar_e (double X, gsl_sf_result * RESULT) + These routines compute the regulated Gamma Function \Gamma^*(x) + for x > 0. The regulated gamma function is given by, + + \Gamma^*(x) = \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x)) + = (1 + (1/12x) + ...) for x \to \infty + and is a useful suggestion of Temme. + + -- Function: double gsl_sf_gammainv (double X) + -- Function: int gsl_sf_gammainv_e (double X, gsl_sf_result * RESULT) + These routines compute the reciprocal of the gamma function, + 1/\Gamma(x) using the real Lanczos method. + + -- Function: int gsl_sf_lngamma_complex_e (double ZR, double ZI, + gsl_sf_result * LNR, gsl_sf_result * ARG) + This routine computes \log(\Gamma(z)) for complex z=z_r+i z_i and + z not a negative integer or zero, using the complex Lanczos + method. The returned parameters are lnr = \log|\Gamma(z)| and arg + = \arg(\Gamma(z)) in (-\pi,\pi]. Note that the phase part (ARG) + is not well-determined when |z| is very large, due to inevitable + roundoff in restricting to (-\pi,\pi]. This will result in a + `GSL_ELOSS' error when it occurs. The absolute value part (LNR), + however, never suffers from loss of precision. + + +File: gsl-ref.info, Node: Factorials, Next: Pochhammer Symbol, Prev: Gamma Functions, Up: Gamma and Beta Functions + +7.19.2 Factorials +----------------- + +Although factorials can be computed from the Gamma function, using the +relation n! = \Gamma(n+1) for non-negative integer n, it is usually +more efficient to call the functions in this section, particularly for +small values of n, whose factorial values are maintained in hardcoded +tables. + + -- Function: double gsl_sf_fact (unsigned int N) + -- Function: int gsl_sf_fact_e (unsigned int N, gsl_sf_result * RESULT) + These routines compute the factorial n!. The factorial is related + to the Gamma function by n! = \Gamma(n+1). The maximum value of n + such that n! is not considered an overflow is given by the macro + `GSL_SF_FACT_NMAX' and is 170. + + -- Function: double gsl_sf_doublefact (unsigned int N) + -- Function: int gsl_sf_doublefact_e (unsigned int N, gsl_sf_result * + RESULT) + These routines compute the double factorial n!! = n(n-2)(n-4) + \dots. The maximum value of n such that n!! is not considered an + overflow is given by the macro `GSL_SF_DOUBLEFACT_NMAX' and is 297. + + -- Function: double gsl_sf_lnfact (unsigned int N) + -- Function: int gsl_sf_lnfact_e (unsigned int N, gsl_sf_result * + RESULT) + These routines compute the logarithm of the factorial of N, + \log(n!). The algorithm is faster than computing \ln(\Gamma(n+1)) + via `gsl_sf_lngamma' for n < 170, but defers for larger N. + + -- Function: double gsl_sf_lndoublefact (unsigned int N) + -- Function: int gsl_sf_lndoublefact_e (unsigned int N, gsl_sf_result + * RESULT) + These routines compute the logarithm of the double factorial of N, + \log(n!!). + + -- Function: double gsl_sf_choose (unsigned int N, unsigned int M) + -- Function: int gsl_sf_choose_e (unsigned int N, unsigned int M, + gsl_sf_result * RESULT) + These routines compute the combinatorial factor `n choose m' = + n!/(m!(n-m)!) + + -- Function: double gsl_sf_lnchoose (unsigned int N, unsigned int M) + -- Function: int gsl_sf_lnchoose_e (unsigned int N, unsigned int M, + gsl_sf_result * RESULT) + These routines compute the logarithm of `n choose m'. This is + equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!). + + -- Function: double gsl_sf_taylorcoeff (int N, double X) + -- Function: int gsl_sf_taylorcoeff_e (int N, double X, gsl_sf_result + * RESULT) + These routines compute the Taylor coefficient x^n / n! for x >= 0, + n >= 0. + + +File: gsl-ref.info, Node: Pochhammer Symbol, Next: Incomplete Gamma Functions, Prev: Factorials, Up: Gamma and Beta Functions + +7.19.3 Pochhammer Symbol +------------------------ + + -- Function: double gsl_sf_poch (double A, double X) + -- Function: int gsl_sf_poch_e (double A, double X, gsl_sf_result * + RESULT) + These routines compute the Pochhammer symbol (a)_x = \Gamma(a + + x)/\Gamma(a), subject to a and a+x not being negative integers or + zero. The Pochhammer symbol is also known as the Apell symbol and + sometimes written as (a,x). + + -- Function: double gsl_sf_lnpoch (double A, double X) + -- Function: int gsl_sf_lnpoch_e (double A, double X, gsl_sf_result * + RESULT) + These routines compute the logarithm of the Pochhammer symbol, + \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)) for a > 0, a+x > 0. + + -- Function: int gsl_sf_lnpoch_sgn_e (double A, double X, + gsl_sf_result * RESULT, double * SGN) + These routines compute the sign of the Pochhammer symbol and the + logarithm of its magnitude. The computed parameters are result = + \log(|(a)_x|) and sgn = \sgn((a)_x) where (a)_x = \Gamma(a + + x)/\Gamma(a), subject to a, a+x not being negative integers or + zero. + + -- Function: double gsl_sf_pochrel (double A, double X) + -- Function: int gsl_sf_pochrel_e (double A, double X, gsl_sf_result * + RESULT) + These routines compute the relative Pochhammer symbol ((a)_x - + 1)/x where (a)_x = \Gamma(a + x)/\Gamma(a). + + +File: gsl-ref.info, Node: Incomplete Gamma Functions, Next: Beta Functions, Prev: Pochhammer Symbol, Up: Gamma and Beta Functions + +7.19.4 Incomplete Gamma Functions +--------------------------------- + + -- Function: double gsl_sf_gamma_inc (double A, double X) + -- Function: int gsl_sf_gamma_inc_e (double A, double X, gsl_sf_result + * RESULT) + These functions compute the unnormalized incomplete Gamma Function + \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x + >= 0. + + -- Function: double gsl_sf_gamma_inc_Q (double A, double X) + -- Function: int gsl_sf_gamma_inc_Q_e (double A, double X, + gsl_sf_result * RESULT) + These routines compute the normalized incomplete Gamma Function + Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0, + x >= 0. + + -- Function: double gsl_sf_gamma_inc_P (double A, double X) + -- Function: int gsl_sf_gamma_inc_P_e (double A, double X, + gsl_sf_result * RESULT) + These routines compute the complementary normalized incomplete + Gamma Function P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt + t^{a-1} \exp(-t) for a > 0, x >= 0. + + Note that Abramowitz & Stegun call P(a,x) the incomplete gamma + function (section 6.5). + + +File: gsl-ref.info, Node: Beta Functions, Next: Incomplete Beta Function, Prev: Incomplete Gamma Functions, Up: Gamma and Beta Functions + +7.19.5 Beta Functions +--------------------- + + -- Function: double gsl_sf_beta (double A, double B) + -- Function: int gsl_sf_beta_e (double A, double B, gsl_sf_result * + RESULT) + These routines compute the Beta Function, B(a,b) = + \Gamma(a)\Gamma(b)/\Gamma(a+b) subject to a and b not being + negative integers. + + -- Function: double gsl_sf_lnbeta (double A, double B) + -- Function: int gsl_sf_lnbeta_e (double A, double B, gsl_sf_result * + RESULT) + These routines compute the logarithm of the Beta Function, + \log(B(a,b)) subject to a and b not being negative integers. + + +File: gsl-ref.info, Node: Incomplete Beta Function, Prev: Beta Functions, Up: Gamma and Beta Functions + +7.19.6 Incomplete Beta Function +------------------------------- + + -- Function: double gsl_sf_beta_inc (double A, double B, double X) + -- Function: int gsl_sf_beta_inc_e (double A, double B, double X, + gsl_sf_result * RESULT) + These routines compute the normalized incomplete Beta function + I_x(a,b)=B_x(a,b)/B(a,b) where B_x(a,b) = \int_0^x t^{a-1} + (1-t)^{b-1} dt for a > 0, b > 0, and 0 <= x <= 1. + + +File: gsl-ref.info, Node: Gegenbauer Functions, Next: Hypergeometric Functions, Prev: Gamma and Beta Functions, Up: Special Functions + +7.20 Gegenbauer Functions +========================= + +The Gegenbauer polynomials are defined in Abramowitz & Stegun, Chapter +22, where they are known as Ultraspherical polynomials. The functions +described in this section are declared in the header file +`gsl_sf_gegenbauer.h'. + + -- Function: double gsl_sf_gegenpoly_1 (double LAMBDA, double X) + -- Function: double gsl_sf_gegenpoly_2 (double LAMBDA, double X) + -- Function: double gsl_sf_gegenpoly_3 (double LAMBDA, double X) + -- Function: int gsl_sf_gegenpoly_1_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + -- Function: int gsl_sf_gegenpoly_2_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + -- Function: int gsl_sf_gegenpoly_3_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + These functions evaluate the Gegenbauer polynomials + C^{(\lambda)}_n(x) using explicit representations for n =1, 2, 3. + + -- Function: double gsl_sf_gegenpoly_n (int N, double LAMBDA, double X) + -- Function: int gsl_sf_gegenpoly_n_e (int N, double LAMBDA, double X, + gsl_sf_result * RESULT) + These functions evaluate the Gegenbauer polynomial + C^{(\lambda)}_n(x) for a specific value of N, LAMBDA, X subject to + \lambda > -1/2, n >= 0. + + -- Function: int gsl_sf_gegenpoly_array (int NMAX, double LAMBDA, + double X, double RESULT_ARRAY[]) + This function computes an array of Gegenbauer polynomials + C^{(\lambda)}_n(x) for n = 0, 1, 2, \dots, nmax, subject to + \lambda > -1/2, nmax >= 0. + + +File: gsl-ref.info, Node: Hypergeometric Functions, Next: Laguerre Functions, Prev: Gegenbauer Functions, Up: Special Functions + +7.21 Hypergeometric Functions +============================= + +Hypergeometric functions are described in Abramowitz & Stegun, Chapters +13 and 15. These functions are declared in the header file +`gsl_sf_hyperg.h'. + + -- Function: double gsl_sf_hyperg_0F1 (double C, double X) + -- Function: int gsl_sf_hyperg_0F1_e (double C, double X, + gsl_sf_result * RESULT) + These routines compute the hypergeometric function 0F1(c,x). + + -- Function: double gsl_sf_hyperg_1F1_int (int M, int N, double X) + -- Function: int gsl_sf_hyperg_1F1_int_e (int M, int N, double X, + gsl_sf_result * RESULT) + These routines compute the confluent hypergeometric function + 1F1(m,n,x) = M(m,n,x) for integer parameters M, N. + + -- Function: double gsl_sf_hyperg_1F1 (double A, double B, double X) + -- Function: int gsl_sf_hyperg_1F1_e (double A, double B, double X, + gsl_sf_result * RESULT) + These routines compute the confluent hypergeometric function + 1F1(a,b,x) = M(a,b,x) for general parameters A, B. + + -- Function: double gsl_sf_hyperg_U_int (int M, int N, double X) + -- Function: int gsl_sf_hyperg_U_int_e (int M, int N, double X, + gsl_sf_result * RESULT) + These routines compute the confluent hypergeometric function + U(m,n,x) for integer parameters M, N. + + -- Function: int gsl_sf_hyperg_U_int_e10_e (int M, int N, double X, + gsl_sf_result_e10 * RESULT) + This routine computes the confluent hypergeometric function + U(m,n,x) for integer parameters M, N using the `gsl_sf_result_e10' + type to return a result with extended range. + + -- Function: double gsl_sf_hyperg_U (double A, double B, double X) + -- Function: int gsl_sf_hyperg_U_e (double A, double B, double X, + gsl_sf_result * RESULT) + These routines compute the confluent hypergeometric function + U(a,b,x). + + -- Function: int gsl_sf_hyperg_U_e10_e (double A, double B, double X, + gsl_sf_result_e10 * RESULT) + This routine computes the confluent hypergeometric function + U(a,b,x) using the `gsl_sf_result_e10' type to return a result + with extended range. + + -- Function: double gsl_sf_hyperg_2F1 (double A, double B, double C, + double X) + -- Function: int gsl_sf_hyperg_2F1_e (double A, double B, double C, + double X, gsl_sf_result * RESULT) + These routines compute the Gauss hypergeometric function + 2F1(a,b,c,x) for |x| < 1. + + If the arguments (a,b,c,x) are too close to a singularity then the + function can return the error code `GSL_EMAXITER' when the series + approximation converges too slowly. This occurs in the region of + x=1, c - a - b = m for integer m. + + -- Function: double gsl_sf_hyperg_2F1_conj (double AR, double AI, + double C, double X) + -- Function: int gsl_sf_hyperg_2F1_conj_e (double AR, double AI, + double C, double X, gsl_sf_result * RESULT) + These routines compute the Gauss hypergeometric function 2F1(a_R + + i a_I, a_R - i a_I, c, x) with complex parameters for |x| < 1. + exceptions: + + -- Function: double gsl_sf_hyperg_2F1_renorm (double A, double B, + double C, double X) + -- Function: int gsl_sf_hyperg_2F1_renorm_e (double A, double B, + double C, double X, gsl_sf_result * RESULT) + These routines compute the renormalized Gauss hypergeometric + function 2F1(a,b,c,x) / \Gamma(c) for |x| < 1. + + -- Function: double gsl_sf_hyperg_2F1_conj_renorm (double AR, double + AI, double C, double X) + -- Function: int gsl_sf_hyperg_2F1_conj_renorm_e (double AR, double + AI, double C, double X, gsl_sf_result * RESULT) + These routines compute the renormalized Gauss hypergeometric + function 2F1(a_R + i a_I, a_R - i a_I, c, x) / \Gamma(c) for |x| < + 1. + + -- Function: double gsl_sf_hyperg_2F0 (double A, double B, double X) + -- Function: int gsl_sf_hyperg_2F0_e (double A, double B, double X, + gsl_sf_result * RESULT) + These routines compute the hypergeometric function 2F0(a,b,x). + The series representation is a divergent hypergeometric series. + However, for x < 0 we have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x) + + +File: gsl-ref.info, Node: Laguerre Functions, Next: Lambert W Functions, Prev: Hypergeometric Functions, Up: Special Functions + +7.22 Laguerre Functions +======================= + +The generalized Laguerre polynomials are defined in terms of confluent +hypergeometric functions as L^a_n(x) = ((a+1)_n / n!) 1F1(-n,a+1,x), +and are sometimes referred to as the associated Laguerre polynomials. +They are related to the plain Laguerre polynomials L_n(x) by L^0_n(x) = +L_n(x) and L^k_n(x) = (-1)^k (d^k/dx^k) L_(n+k)(x). For more +information see Abramowitz & Stegun, Chapter 22. + + The functions described in this section are declared in the header +file `gsl_sf_laguerre.h'. + + -- Function: double gsl_sf_laguerre_1 (double A, double X) + -- Function: double gsl_sf_laguerre_2 (double A, double X) + -- Function: double gsl_sf_laguerre_3 (double A, double X) + -- Function: int gsl_sf_laguerre_1_e (double A, double X, + gsl_sf_result * RESULT) + -- Function: int gsl_sf_laguerre_2_e (double A, double X, + gsl_sf_result * RESULT) + -- Function: int gsl_sf_laguerre_3_e (double A, double X, + gsl_sf_result * RESULT) + These routines evaluate the generalized Laguerre polynomials + L^a_1(x), L^a_2(x), L^a_3(x) using explicit representations. + + -- Function: double gsl_sf_laguerre_n (const int N, const double A, + const double X) + -- Function: int gsl_sf_laguerre_n_e (int N, double A, double X, + gsl_sf_result * RESULT) + These routines evaluate the generalized Laguerre polynomials + L^a_n(x) for a > -1, n >= 0. + + + +File: gsl-ref.info, Node: Lambert W Functions, Next: Legendre Functions and Spherical Harmonics, Prev: Laguerre Functions, Up: Special Functions + +7.23 Lambert W Functions +======================== + +Lambert's W functions, W(x), are defined to be solutions of the +equation W(x) \exp(W(x)) = x. This function has multiple branches for x +< 0; however, it has only two real-valued branches. We define W_0(x) to +be the principal branch, where W > -1 for x < 0, and W_{-1}(x) to be +the other real branch, where W < -1 for x < 0. The Lambert functions +are declared in the header file `gsl_sf_lambert.h'. + + -- Function: double gsl_sf_lambert_W0 (double X) + -- Function: int gsl_sf_lambert_W0_e (double X, gsl_sf_result * RESULT) + These compute the principal branch of the Lambert W function, + W_0(x). + + -- Function: double gsl_sf_lambert_Wm1 (double X) + -- Function: int gsl_sf_lambert_Wm1_e (double X, gsl_sf_result * + RESULT) + These compute the secondary real-valued branch of the Lambert W + function, W_{-1}(x). + + +File: gsl-ref.info, Node: Legendre Functions and Spherical Harmonics, Next: Logarithm and Related Functions, Prev: Lambert W Functions, Up: Special Functions + +7.24 Legendre Functions and Spherical Harmonics +=============================================== + +The Legendre Functions and Legendre Polynomials are described in +Abramowitz & Stegun, Chapter 8. These functions are declared in the +header file `gsl_sf_legendre.h'. + +* Menu: + +* Legendre Polynomials:: +* Associated Legendre Polynomials and Spherical Harmonics:: +* Conical Functions:: +* Radial Functions for Hyperbolic Space:: + + +File: gsl-ref.info, Node: Legendre Polynomials, Next: Associated Legendre Polynomials and Spherical Harmonics, Up: Legendre Functions and Spherical Harmonics + +7.24.1 Legendre Polynomials +--------------------------- + + -- Function: double gsl_sf_legendre_P1 (double X) + -- Function: double gsl_sf_legendre_P2 (double X) + -- Function: double gsl_sf_legendre_P3 (double X) + -- Function: int gsl_sf_legendre_P1_e (double X, gsl_sf_result * + RESULT) + -- Function: int gsl_sf_legendre_P2_e (double X, gsl_sf_result * + RESULT) + -- Function: int gsl_sf_legendre_P3_e (double X, gsl_sf_result * + RESULT) + These functions evaluate the Legendre polynomials P_l(x) using + explicit representations for l=1, 2, 3. + + -- Function: double gsl_sf_legendre_Pl (int L, double X) + -- Function: int gsl_sf_legendre_Pl_e (int L, double X, gsl_sf_result + * RESULT) + These functions evaluate the Legendre polynomial P_l(x) for a + specific value of L, X subject to l >= 0, |x| <= 1 + + -- Function: int gsl_sf_legendre_Pl_array (int LMAX, double X, double + RESULT_ARRAY[]) + -- Function: int gsl_sf_legendre_Pl_deriv_array (int LMAX, double X, + double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[]) + These functions compute an array of Legendre polynomials P_l(x), + and optionally their derivatives dP_l(x)/dx, for l = 0, \dots, + lmax, |x| <= 1 + + -- Function: double gsl_sf_legendre_Q0 (double X) + -- Function: int gsl_sf_legendre_Q0_e (double X, gsl_sf_result * + RESULT) + These routines compute the Legendre function Q_0(x) for x > -1, x + != 1. + + -- Function: double gsl_sf_legendre_Q1 (double X) + -- Function: int gsl_sf_legendre_Q1_e (double X, gsl_sf_result * + RESULT) + These routines compute the Legendre function Q_1(x) for x > -1, x + != 1. + + -- Function: double gsl_sf_legendre_Ql (int L, double X) + -- Function: int gsl_sf_legendre_Ql_e (int L, double X, gsl_sf_result + * RESULT) + These routines compute the Legendre function Q_l(x) for x > -1, x + != 1 and l >= 0. + + +File: gsl-ref.info, Node: Associated Legendre Polynomials and Spherical Harmonics, Next: Conical Functions, Prev: Legendre Polynomials, Up: Legendre Functions and Spherical Harmonics + +7.24.2 Associated Legendre Polynomials and Spherical Harmonics +-------------------------------------------------------------- + +The following functions compute the associated Legendre Polynomials +P_l^m(x). Note that this function grows combinatorially with l and can +overflow for l larger than about 150. There is no trouble for small m, +but overflow occurs when m and l are both large. Rather than allow +overflows, these functions refuse to calculate P_l^m(x) and return +`GSL_EOVRFLW' when they can sense that l and m are too big. + + If you want to calculate a spherical harmonic, then _do not_ use +these functions. Instead use `gsl_sf_legendre_sphPlm' below, which +uses a similar recursion, but with the normalized functions. + + -- Function: double gsl_sf_legendre_Plm (int L, int M, double X) + -- Function: int gsl_sf_legendre_Plm_e (int L, int M, double X, + gsl_sf_result * RESULT) + These routines compute the associated Legendre polynomial P_l^m(x) + for m >= 0, l >= m, |x| <= 1. + + -- Function: int gsl_sf_legendre_Plm_array (int LMAX, int M, double X, + double RESULT_ARRAY[]) + -- Function: int gsl_sf_legendre_Plm_deriv_array (int LMAX, int M, + double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[]) + These functions compute an array of Legendre polynomials P_l^m(x), + and optionally their derivatives dP_l^m(x)/dx, for m >= 0, l = + |m|, ..., lmax, |x| <= 1. + + -- Function: double gsl_sf_legendre_sphPlm (int L, int M, double X) + -- Function: int gsl_sf_legendre_sphPlm_e (int L, int M, double X, + gsl_sf_result * RESULT) + These routines compute the normalized associated Legendre + polynomial $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$ + suitable for use in spherical harmonics. The parameters must + satisfy m >= 0, l >= m, |x| <= 1. Theses routines avoid the + overflows that occur for the standard normalization of P_l^m(x). + + -- Function: int gsl_sf_legendre_sphPlm_array (int LMAX, int M, double + X, double RESULT_ARRAY[]) + -- Function: int gsl_sf_legendre_sphPlm_deriv_array (int LMAX, int M, + double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[]) + These functions compute an array of normalized associated Legendre + functions $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$, + and optionally their derivatives, for m >= 0, l = |m|, ..., lmax, + |x| <= 1.0 + + -- Function: int gsl_sf_legendre_array_size (const int LMAX, const int + M) + This function returns the size of RESULT_ARRAY[] needed for the + array versions of P_l^m(x), LMAX - M + 1. + + +File: gsl-ref.info, Node: Conical Functions, Next: Radial Functions for Hyperbolic Space, Prev: Associated Legendre Polynomials and Spherical Harmonics, Up: Legendre Functions and Spherical Harmonics + +7.24.3 Conical Functions +------------------------ + +The Conical Functions P^\mu_{-(1/2)+i\lambda}(x) and +Q^\mu_{-(1/2)+i\lambda} are described in Abramowitz & Stegun, Section +8.12. + + -- Function: double gsl_sf_conicalP_half (double LAMBDA, double X) + -- Function: int gsl_sf_conicalP_half_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + These routines compute the irregular Spherical Conical Function + P^{1/2}_{-1/2 + i \lambda}(x) for x > -1. + + -- Function: double gsl_sf_conicalP_mhalf (double LAMBDA, double X) + -- Function: int gsl_sf_conicalP_mhalf_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + These routines compute the regular Spherical Conical Function + P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1. + + -- Function: double gsl_sf_conicalP_0 (double LAMBDA, double X) + -- Function: int gsl_sf_conicalP_0_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + These routines compute the conical function P^0_{-1/2 + i + \lambda}(x) for x > -1. + + -- Function: double gsl_sf_conicalP_1 (double LAMBDA, double X) + -- Function: int gsl_sf_conicalP_1_e (double LAMBDA, double X, + gsl_sf_result * RESULT) + These routines compute the conical function P^1_{-1/2 + i + \lambda}(x) for x > -1. + + -- Function: double gsl_sf_conicalP_sph_reg (int L, double LAMBDA, + double X) + -- Function: int gsl_sf_conicalP_sph_reg_e (int L, double LAMBDA, + double X, gsl_sf_result * RESULT) + These routines compute the Regular Spherical Conical Function + P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1. + + -- Function: double gsl_sf_conicalP_cyl_reg (int M, double LAMBDA, + double X) + -- Function: int gsl_sf_conicalP_cyl_reg_e (int M, double LAMBDA, + double X, gsl_sf_result * RESULT) + These routines compute the Regular Cylindrical Conical Function + P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1. + + +File: gsl-ref.info, Node: Radial Functions for Hyperbolic Space, Prev: Conical Functions, Up: Legendre Functions and Spherical Harmonics + +7.24.4 Radial Functions for Hyperbolic Space +-------------------------------------------- + +The following spherical functions are specializations of Legendre +functions which give the regular eigenfunctions of the Laplacian on a +3-dimensional hyperbolic space H3d. Of particular interest is the flat +limit, \lambda \to \infty, \eta \to 0, \lambda\eta fixed. + + -- Function: double gsl_sf_legendre_H3d_0 (double LAMBDA, double ETA) + -- Function: int gsl_sf_legendre_H3d_0_e (double LAMBDA, double ETA, + gsl_sf_result * RESULT) + These routines compute the zeroth radial eigenfunction of the + Laplacian on the 3-dimensional hyperbolic space, + L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta)) + for \eta >= 0. In the flat limit this takes the form + L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta). + + -- Function: double gsl_sf_legendre_H3d_1 (double LAMBDA, double ETA) + -- Function: int gsl_sf_legendre_H3d_1_e (double LAMBDA, double ETA, + gsl_sf_result * RESULT) + These routines compute the first radial eigenfunction of the + Laplacian on the 3-dimensional hyperbolic space, + L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1} \sin(\lambda + \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda + \cot(\lambda\eta)) for \eta >= 0. In the flat limit this takes + the form L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta). + + -- Function: double gsl_sf_legendre_H3d (int L, double LAMBDA, double + ETA) + -- Function: int gsl_sf_legendre_H3d_e (int L, double LAMBDA, double + ETA, gsl_sf_result * RESULT) + These routines compute the L-th radial eigenfunction of the + Laplacian on the 3-dimensional hyperbolic space \eta >= 0, l >= 0. + In the flat limit this takes the form L^{H3d}_l(\lambda,\eta) = + j_l(\lambda\eta). + + -- Function: int gsl_sf_legendre_H3d_array (int LMAX, double LAMBDA, + double ETA, double RESULT_ARRAY[]) + This function computes an array of radial eigenfunctions + L^{H3d}_l(\lambda, \eta) for 0 <= l <= lmax. + + +File: gsl-ref.info, Node: Logarithm and Related Functions, Next: Mathieu Functions, Prev: Legendre Functions and Spherical Harmonics, Up: Special Functions + +7.25 Logarithm and Related Functions +==================================== + +Information on the properties of the Logarithm function can be found in +Abramowitz & Stegun, Chapter 4. The functions described in this section +are declared in the header file `gsl_sf_log.h'. + + -- Function: double gsl_sf_log (double X) + -- Function: int gsl_sf_log_e (double X, gsl_sf_result * RESULT) + These routines compute the logarithm of X, \log(x), for x > 0. + + -- Function: double gsl_sf_log_abs (double X) + -- Function: int gsl_sf_log_abs_e (double X, gsl_sf_result * RESULT) + These routines compute the logarithm of the magnitude of X, + \log(|x|), for x \ne 0. + + -- Function: int gsl_sf_complex_log_e (double ZR, double ZI, + gsl_sf_result * LNR, gsl_sf_result * THETA) + This routine computes the complex logarithm of z = z_r + i z_i. + The results are returned as LNR, THETA such that \exp(lnr + i + \theta) = z_r + i z_i, where \theta lies in the range [-\pi,\pi]. + + -- Function: double gsl_sf_log_1plusx (double X) + -- Function: int gsl_sf_log_1plusx_e (double X, gsl_sf_result * RESULT) + These routines compute \log(1 + x) for x > -1 using an algorithm + that is accurate for small x. + + -- Function: double gsl_sf_log_1plusx_mx (double X) + -- Function: int gsl_sf_log_1plusx_mx_e (double X, gsl_sf_result * + RESULT) + These routines compute \log(1 + x) - x for x > -1 using an + algorithm that is accurate for small x. + + +File: gsl-ref.info, Node: Mathieu Functions, Next: Power Function, Prev: Logarithm and Related Functions, Up: Special Functions + +7.26 Mathieu Functions +====================== + +The routines described in this section compute the angular and radial +Mathieu functions, and their characteristic values. Mathieu functions +are the solutions of the following two differential equations: + + d^2y/dv^2 + (a - 2q\cos 2v)y = 0 + d^2f/du^2 - (a - 2q\cosh 2u)f = 0 + +The angular Mathieu functions ce_r(x,q), se_r(x,q) are the even and odd +periodic solutions of the first equation, which is known as Mathieu's +equation. These exist only for the discrete sequence of characteristic +values a=a_r(q) (even-periodic) and a=b_r(q) (odd-periodic). + + The radial Mathieu functions Mc^{(j)}_{r}(z,q), Ms^{(j)}_{r}(z,q) +are the solutions of the second equation, which is referred to as +Mathieu's modified equation. The radial Mathieu functions of the +first, second, third and fourth kind are denoted by the parameter j, +which takes the value 1, 2, 3 or 4. + + For more information on the Mathieu functions, see Abramowitz and +Stegun, Chapter 20. These functions are defined in the header file +`gsl_sf_mathieu.h'. + +* Menu: + +* Mathieu Function Workspace:: +* Mathieu Function Characteristic Values:: +* Angular Mathieu Functions:: +* Radial Mathieu Functions:: + + +File: gsl-ref.info, Node: Mathieu Function Workspace, Next: Mathieu Function Characteristic Values, Up: Mathieu Functions + +7.26.1 Mathieu Function Workspace +--------------------------------- + +The Mathieu functions can be computed for a single order or for +multiple orders, using array-based routines. The array-based routines +require a preallocated workspace. + + -- Function: gsl_sf_mathieu_workspace * gsl_sf_mathieu_alloc (size_t + N, double QMAX) + This function returns a workspace for the array versions of the + Mathieu routines. The arguments N and QMAX specify the maximum + order and q-value of Mathieu functions which can be computed with + this workspace. + + + -- Function: void gsl_sf_mathieu_free (gsl_sf_mathieu_workspace *WORK) + This function frees the workspace WORK. + + +File: gsl-ref.info, Node: Mathieu Function Characteristic Values, Next: Angular Mathieu Functions, Prev: Mathieu Function Workspace, Up: Mathieu Functions + +7.26.2 Mathieu Function Characteristic Values +--------------------------------------------- + + -- Function: int gsl_sf_mathieu_a (int N, double Q, gsl_sf_result + *RESULT) + -- Function: int gsl_sf_mathieu_b (int N, double Q, gsl_sf_result + *RESULT) + These routines compute the characteristic values a_n(q), b_n(q) of + the Mathieu functions ce_n(q,x) and se_n(q,x), respectively. + + -- Function: int gsl_sf_mathieu_a_array (int ORDER_MIN, int ORDER_MAX, + double Q, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + -- Function: int gsl_sf_mathieu_b_array (int ORDER_MIN, int ORDER_MAX, + double Q, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + These routines compute a series of Mathieu characteristic values + a_n(q), b_n(q) for n from ORDER_MIN to ORDER_MAX inclusive, + storing the results in the array RESULT_ARRAY. + + +File: gsl-ref.info, Node: Angular Mathieu Functions, Next: Radial Mathieu Functions, Prev: Mathieu Function Characteristic Values, Up: Mathieu Functions + +7.26.3 Angular Mathieu Functions +-------------------------------- + + -- Function: int gsl_sf_mathieu_ce (int N, double Q, double X, + gsl_sf_result *RESULT) + -- Function: int gsl_sf_mathieu_se (int N, double Q, double X, + gsl_sf_result *RESULT) + These routines compute the angular Mathieu functions ce_n(q,x) and + se_n(q,x), respectively. + + -- Function: int gsl_sf_mathieu_ce_array (int NMIN, int NMAX, double + Q, double X, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + -- Function: int gsl_sf_mathieu_se_array (int NMIN, int NMAX, double + Q, double X, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + These routines compute a series of the angular Mathieu functions + ce_n(q,x) and se_n(q,x) of order n from NMIN to NMAX inclusive, + storing the results in the array RESULT_ARRAY. + + +File: gsl-ref.info, Node: Radial Mathieu Functions, Prev: Angular Mathieu Functions, Up: Mathieu Functions + +7.26.4 Radial Mathieu Functions +------------------------------- + + -- Function: int gsl_sf_mathieu_Mc (int J, int N, double Q, double X, + gsl_sf_result *RESULT) + -- Function: int gsl_sf_mathieu_Ms (int J, int N, double Q, double X, + gsl_sf_result *RESULT) + These routines compute the radial J-th kind Mathieu functions + Mc_n^{(j)}(q,x) and Ms_n^{(j)}(q,x) of order N. + + The allowed values of J are 1 and 2. The functions for j = 3,4 + can be computed as M_n^{(3)} = M_n^{(1)} + iM_n^{(2)} and + M_n^{(4)} = M_n^{(1)} - iM_n^{(2)}, where M_n^{(j)} = Mc_n^{(j)} or + Ms_n^{(j)}. + + -- Function: int gsl_sf_mathieu_Mc_array (int J, int NMIN, int NMAX, + double Q, double X, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + -- Function: int gsl_sf_mathieu_Ms_array (int J, int NMIN, int NMAX, + double Q, double X, gsl_sf_mathieu_workspace *WORK, double + RESULT_ARRAY[]) + These routines compute a series of the radial Mathieu functions of + kind J, with order from NMIN to NMAX inclusive, storing the + results in the array RESULT_ARRAY. + + +File: gsl-ref.info, Node: Power Function, Next: Psi (Digamma) Function, Prev: Mathieu Functions, Up: Special Functions + +7.27 Power Function +=================== + +The following functions are equivalent to the function `gsl_pow_int' +(*note Small integer powers::) with an error estimate. These functions +are declared in the header file `gsl_sf_pow_int.h'. + + -- Function: double gsl_sf_pow_int (double X, int N) + -- Function: int gsl_sf_pow_int_e (double X, int N, gsl_sf_result * + RESULT) + These routines compute the power x^n for integer N. The power is + computed using the minimum number of multiplications. For example, + x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications. + For reasons of efficiency, these functions do not check for + overflow or underflow conditions. + + #include <gsl/gsl_sf_pow_int.h> + /* compute 3.0**12 */ + double y = gsl_sf_pow_int(3.0, 12); + + +File: gsl-ref.info, Node: Psi (Digamma) Function, Next: Synchrotron Functions, Prev: Power Function, Up: Special Functions + +7.28 Psi (Digamma) Function +=========================== + +The polygamma functions of order n are defined by + + \psi^{(n)}(x) = (d/dx)^n \psi(x) = (d/dx)^{n+1} \log(\Gamma(x)) + +where \psi(x) = \Gamma'(x)/\Gamma(x) is known as the digamma function. +These functions are declared in the header file `gsl_sf_psi.h'. + +* Menu: + +* Digamma Function:: +* Trigamma Function:: +* Polygamma Function:: + + +File: gsl-ref.info, Node: Digamma Function, Next: Trigamma Function, Up: Psi (Digamma) Function + +7.28.1 Digamma Function +----------------------- + + -- Function: double gsl_sf_psi_int (int N) + -- Function: int gsl_sf_psi_int_e (int N, gsl_sf_result * RESULT) + These routines compute the digamma function \psi(n) for positive + integer N. The digamma function is also called the Psi function. + + -- Function: double gsl_sf_psi (double X) + -- Function: int gsl_sf_psi_e (double X, gsl_sf_result * RESULT) + These routines compute the digamma function \psi(x) for general x, + x \ne 0. + + -- Function: double gsl_sf_psi_1piy (double Y) + -- Function: int gsl_sf_psi_1piy_e (double Y, gsl_sf_result * RESULT) + These routines compute the real part of the digamma function on + the line 1+i y, \Re[\psi(1 + i y)]. + + +File: gsl-ref.info, Node: Trigamma Function, Next: Polygamma Function, Prev: Digamma Function, Up: Psi (Digamma) Function + +7.28.2 Trigamma Function +------------------------ + + -- Function: double gsl_sf_psi_1_int (int N) + -- Function: int gsl_sf_psi_1_int_e (int N, gsl_sf_result * RESULT) + These routines compute the Trigamma function \psi'(n) for positive + integer n. + + -- Function: double gsl_sf_psi_1 (double X) + -- Function: int gsl_sf_psi_1_e (double X, gsl_sf_result * RESULT) + These routines compute the Trigamma function \psi'(x) for general + x. + + +File: gsl-ref.info, Node: Polygamma Function, Prev: Trigamma Function, Up: Psi (Digamma) Function + +7.28.3 Polygamma Function +------------------------- + + -- Function: double gsl_sf_psi_n (int N, double X) + -- Function: int gsl_sf_psi_n_e (int N, double X, gsl_sf_result * + RESULT) + These routines compute the polygamma function \psi^{(n)}(x) for n + >= 0, x > 0. + + +File: gsl-ref.info, Node: Synchrotron Functions, Next: Transport Functions, Prev: Psi (Digamma) Function, Up: Special Functions + +7.29 Synchrotron Functions +========================== + +The functions described in this section are declared in the header file +`gsl_sf_synchrotron.h'. + + -- Function: double gsl_sf_synchrotron_1 (double X) + -- Function: int gsl_sf_synchrotron_1_e (double X, gsl_sf_result * + RESULT) + These routines compute the first synchrotron function x + \int_x^\infty dt K_{5/3}(t) for x >= 0. + + -- Function: double gsl_sf_synchrotron_2 (double X) + -- Function: int gsl_sf_synchrotron_2_e (double X, gsl_sf_result * + RESULT) + These routines compute the second synchrotron function x + K_{2/3}(x) for x >= 0. + + +File: gsl-ref.info, Node: Transport Functions, Next: Trigonometric Functions, Prev: Synchrotron Functions, Up: Special Functions + +7.30 Transport Functions +======================== + +The transport functions J(n,x) are defined by the integral +representations J(n,x) := \int_0^x dt t^n e^t /(e^t - 1)^2. They are +declared in the header file `gsl_sf_transport.h'. + + -- Function: double gsl_sf_transport_2 (double X) + -- Function: int gsl_sf_transport_2_e (double X, gsl_sf_result * + RESULT) + These routines compute the transport function J(2,x). + + -- Function: double gsl_sf_transport_3 (double X) + -- Function: int gsl_sf_transport_3_e (double X, gsl_sf_result * + RESULT) + These routines compute the transport function J(3,x). + + -- Function: double gsl_sf_transport_4 (double X) + -- Function: int gsl_sf_transport_4_e (double X, gsl_sf_result * + RESULT) + These routines compute the transport function J(4,x). + + -- Function: double gsl_sf_transport_5 (double X) + -- Function: int gsl_sf_transport_5_e (double X, gsl_sf_result * + RESULT) + These routines compute the transport function J(5,x). + + +File: gsl-ref.info, Node: Trigonometric Functions, Next: Zeta Functions, Prev: Transport Functions, Up: Special Functions + +7.31 Trigonometric Functions +============================ + +The library includes its own trigonometric functions in order to provide +consistency across platforms and reliable error estimates. These +functions are declared in the header file `gsl_sf_trig.h'. + +* Menu: + +* Circular Trigonometric Functions:: +* Trigonometric Functions for Complex Arguments:: +* Hyperbolic Trigonometric Functions:: +* Conversion Functions:: +* Restriction Functions:: +* Trigonometric Functions With Error Estimates:: + + +File: gsl-ref.info, Node: Circular Trigonometric Functions, Next: Trigonometric Functions for Complex Arguments, Up: Trigonometric Functions + +7.31.1 Circular Trigonometric Functions +--------------------------------------- + + -- Function: double gsl_sf_sin (double X) + -- Function: int gsl_sf_sin_e (double X, gsl_sf_result * RESULT) + These routines compute the sine function \sin(x). + + -- Function: double gsl_sf_cos (double X) + -- Function: int gsl_sf_cos_e (double X, gsl_sf_result * RESULT) + These routines compute the cosine function \cos(x). + + -- Function: double gsl_sf_hypot (double X, double Y) + -- Function: int gsl_sf_hypot_e (double X, double Y, gsl_sf_result * + RESULT) + These routines compute the hypotenuse function \sqrt{x^2 + y^2} + avoiding overflow and underflow. + + -- Function: double gsl_sf_sinc (double X) + -- Function: int gsl_sf_sinc_e (double X, gsl_sf_result * RESULT) + These routines compute \sinc(x) = \sin(\pi x) / (\pi x) for any + value of X. + + +File: gsl-ref.info, Node: Trigonometric Functions for Complex Arguments, Next: Hyperbolic Trigonometric Functions, Prev: Circular Trigonometric Functions, Up: Trigonometric Functions + +7.31.2 Trigonometric Functions for Complex Arguments +---------------------------------------------------- + + -- Function: int gsl_sf_complex_sin_e (double ZR, double ZI, + gsl_sf_result * SZR, gsl_sf_result * SZI) + This function computes the complex sine, \sin(z_r + i z_i) storing + the real and imaginary parts in SZR, SZI. + + -- Function: int gsl_sf_complex_cos_e (double ZR, double ZI, + gsl_sf_result * CZR, gsl_sf_result * CZI) + This function computes the complex cosine, \cos(z_r + i z_i) + storing the real and imaginary parts in SZR, SZI. + + -- Function: int gsl_sf_complex_logsin_e (double ZR, double ZI, + gsl_sf_result * LSZR, gsl_sf_result * LSZI) + This function computes the logarithm of the complex sine, + \log(\sin(z_r + i z_i)) storing the real and imaginary parts in + SZR, SZI. + + +File: gsl-ref.info, Node: Hyperbolic Trigonometric Functions, Next: Conversion Functions, Prev: Trigonometric Functions for Complex Arguments, Up: Trigonometric Functions + +7.31.3 Hyperbolic Trigonometric Functions +----------------------------------------- + + -- Function: double gsl_sf_lnsinh (double X) + -- Function: int gsl_sf_lnsinh_e (double X, gsl_sf_result * RESULT) + These routines compute \log(\sinh(x)) for x > 0. + + -- Function: double gsl_sf_lncosh (double X) + -- Function: int gsl_sf_lncosh_e (double X, gsl_sf_result * RESULT) + These routines compute \log(\cosh(x)) for any X. + + +File: gsl-ref.info, Node: Conversion Functions, Next: Restriction Functions, Prev: Hyperbolic Trigonometric Functions, Up: Trigonometric Functions + +7.31.4 Conversion Functions +--------------------------- + + -- Function: int gsl_sf_polar_to_rect (double R, double THETA, + gsl_sf_result * X, gsl_sf_result * Y); + This function converts the polar coordinates (R,THETA) to + rectilinear coordinates (X,Y), x = r\cos(\theta), y = + r\sin(\theta). + + -- Function: int gsl_sf_rect_to_polar (double X, double Y, + gsl_sf_result * R, gsl_sf_result * THETA) + This function converts the rectilinear coordinates (X,Y) to polar + coordinates (R,THETA), such that x = r\cos(\theta), y = + r\sin(\theta). The argument THETA lies in the range [-\pi, \pi]. + + +File: gsl-ref.info, Node: Restriction Functions, Next: Trigonometric Functions With Error Estimates, Prev: Conversion Functions, Up: Trigonometric Functions + +7.31.5 Restriction Functions +---------------------------- + + -- Function: double gsl_sf_angle_restrict_symm (double THETA) + -- Function: int gsl_sf_angle_restrict_symm_e (double * THETA) + These routines force the angle THETA to lie in the range + (-\pi,\pi]. + + Note that the mathematical value of \pi is slightly greater than + `M_PI', so the machine numbers `M_PI' and `-M_PI' are included in + the range. + + -- Function: double gsl_sf_angle_restrict_pos (double THETA) + -- Function: int gsl_sf_angle_restrict_pos_e (double * THETA) + These routines force the angle THETA to lie in the range [0, 2\pi). + + Note that the mathematical value of 2\pi is slightly greater than + `2*M_PI', so the machine number `2*M_PI' is included in the range. + + + +File: gsl-ref.info, Node: Trigonometric Functions With Error Estimates, Prev: Restriction Functions, Up: Trigonometric Functions + +7.31.6 Trigonometric Functions With Error Estimates +--------------------------------------------------- + + -- Function: int gsl_sf_sin_err_e (double X, double DX, gsl_sf_result + * RESULT) + This routine computes the sine of an angle X with an associated + absolute error DX, \sin(x \pm dx). Note that this function is + provided in the error-handling form only since its purpose is to + compute the propagated error. + + -- Function: int gsl_sf_cos_err_e (double X, double DX, gsl_sf_result + * RESULT) + This routine computes the cosine of an angle X with an associated + absolute error DX, \cos(x \pm dx). Note that this function is + provided in the error-handling form only since its purpose is to + compute the propagated error. + + +File: gsl-ref.info, Node: Zeta Functions, Next: Special Functions Examples, Prev: Trigonometric Functions, Up: Special Functions + +7.32 Zeta Functions +=================== + +The Riemann zeta function is defined in Abramowitz & Stegun, Section +23.2. The functions described in this section are declared in the +header file `gsl_sf_zeta.h'. + +* Menu: + +* Riemann Zeta Function:: +* Riemann Zeta Function Minus One:: +* Hurwitz Zeta Function:: +* Eta Function:: + + +File: gsl-ref.info, Node: Riemann Zeta Function, Next: Riemann Zeta Function Minus One, Up: Zeta Functions + +7.32.1 Riemann Zeta Function +---------------------------- + +The Riemann zeta function is defined by the infinite sum \zeta(s) = +\sum_{k=1}^\infty k^{-s}. + + -- Function: double gsl_sf_zeta_int (int N) + -- Function: int gsl_sf_zeta_int_e (int N, gsl_sf_result * RESULT) + These routines compute the Riemann zeta function \zeta(n) for + integer N, n \ne 1. + + -- Function: double gsl_sf_zeta (double S) + -- Function: int gsl_sf_zeta_e (double S, gsl_sf_result * RESULT) + These routines compute the Riemann zeta function \zeta(s) for + arbitrary S, s \ne 1. + + +File: gsl-ref.info, Node: Riemann Zeta Function Minus One, Next: Hurwitz Zeta Function, Prev: Riemann Zeta Function, Up: Zeta Functions + +7.32.2 Riemann Zeta Function Minus One +-------------------------------------- + +For large positive argument, the Riemann zeta function approaches one. +In this region the fractional part is interesting, and therefore we +need a function to evaluate it explicitly. + + -- Function: double gsl_sf_zetam1_int (int N) + -- Function: int gsl_sf_zetam1_int_e (int N, gsl_sf_result * RESULT) + These routines compute \zeta(n) - 1 for integer N, n \ne 1. + + -- Function: double gsl_sf_zetam1 (double S) + -- Function: int gsl_sf_zetam1_e (double S, gsl_sf_result * RESULT) + These routines compute \zeta(s) - 1 for arbitrary S, s \ne 1. + + +File: gsl-ref.info, Node: Hurwitz Zeta Function, Next: Eta Function, Prev: Riemann Zeta Function Minus One, Up: Zeta Functions + +7.32.3 Hurwitz Zeta Function +---------------------------- + +The Hurwitz zeta function is defined by \zeta(s,q) = \sum_0^\infty +(k+q)^{-s}. + + -- Function: double gsl_sf_hzeta (double S, double Q) + -- Function: int gsl_sf_hzeta_e (double S, double Q, gsl_sf_result * + RESULT) + These routines compute the Hurwitz zeta function \zeta(s,q) for s + > 1, q > 0. + + +File: gsl-ref.info, Node: Eta Function, Prev: Hurwitz Zeta Function, Up: Zeta Functions + +7.32.4 Eta Function +------------------- + +The eta function is defined by \eta(s) = (1-2^{1-s}) \zeta(s). + + -- Function: double gsl_sf_eta_int (int N) + -- Function: int gsl_sf_eta_int_e (int N, gsl_sf_result * RESULT) + These routines compute the eta function \eta(n) for integer N. + + -- Function: double gsl_sf_eta (double S) + -- Function: int gsl_sf_eta_e (double S, gsl_sf_result * RESULT) + These routines compute the eta function \eta(s) for arbitrary S. + + +File: gsl-ref.info, Node: Special Functions Examples, Next: Special Functions References and Further Reading, Prev: Zeta Functions, Up: Special Functions + +7.33 Examples +============= + +The following example demonstrates the use of the error handling form of +the special functions, in this case to compute the Bessel function +J_0(5.0), + + #include <stdio.h> + #include <gsl/gsl_errno.h> + #include <gsl/gsl_sf_bessel.h> + + int + main (void) + { + double x = 5.0; + gsl_sf_result result; + + double expected = -0.17759677131433830434739701; + + int status = gsl_sf_bessel_J0_e (x, &result); + + printf ("status = %s\n", gsl_strerror(status)); + printf ("J0(5.0) = %.18f\n" + " +/- % .18f\n", + result.val, result.err); + printf ("exact = %.18f\n", expected); + return status; + } + +Here are the results of running the program, + + $ ./a.out + status = success + J0(5.0) = -0.177596771314338292 + +/- 0.000000000000000193 + exact = -0.177596771314338292 + +The next program computes the same quantity using the natural form of +the function. In this case the error term RESULT.ERR and return status +are not accessible. + + #include <stdio.h> + #include <gsl/gsl_sf_bessel.h> + + int + main (void) + { + double x = 5.0; + double expected = -0.17759677131433830434739701; + + double y = gsl_sf_bessel_J0 (x); + + printf ("J0(5.0) = %.18f\n", y); + printf ("exact = %.18f\n", expected); + return 0; + } + +The results of the function are the same, + + $ ./a.out + J0(5.0) = -0.177596771314338292 + exact = -0.177596771314338292 + + +File: gsl-ref.info, Node: Special Functions References and Further Reading, Prev: Special Functions Examples, Up: Special Functions + +7.34 References and Further Reading +=================================== + +The library follows the conventions of `Abramowitz & Stegun' where +possible, + Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions' + +The following papers contain information on the algorithms used to +compute the special functions, + MISCFUN: A software package to compute uncommon special functions. + `ACM Trans. Math. Soft.', vol. 22, 1996, 288-301 + + G.N. Watson, A Treatise on the Theory of Bessel Functions, 2nd + Edition (Cambridge University Press, 1944). + + G. Nemeth, Mathematical Approximations of Special Functions, Nova + Science Publishers, ISBN 1-56072-052-2 + + B.C. Carlson, Special Functions of Applied Mathematics (1977) + + W.J. Thompson, Atlas for Computing Mathematical Functions, John + Wiley & Sons, New York (1997). + + Y.Y. Luke, Algorithms for the Computation of Mathematical + Functions, Academic Press, New York (1977). + + + +File: gsl-ref.info, Node: Vectors and Matrices, Next: Permutations, Prev: Special Functions, Up: Top + +8 Vectors and Matrices +********************** + +The functions described in this chapter provide a simple vector and +matrix interface to ordinary C arrays. The memory management of these +arrays is implemented using a single underlying type, known as a block. +By writing your functions in terms of vectors and matrices you can pass +a single structure containing both data and dimensions as an argument +without needing additional function parameters. The structures are +compatible with the vector and matrix formats used by BLAS routines. + +* Menu: + +* Data types:: +* Blocks:: +* Vectors:: +* Matrices:: +* Vector and Matrix References and Further Reading:: + + +File: gsl-ref.info, Node: Data types, Next: Blocks, Up: Vectors and Matrices + +8.1 Data types +============== + +All the functions are available for each of the standard data-types. +The versions for `double' have the prefix `gsl_block', `gsl_vector' and +`gsl_matrix'. Similarly the versions for single-precision `float' +arrays have the prefix `gsl_block_float', `gsl_vector_float' and +`gsl_matrix_float'. The full list of available types is given below, + + gsl_block double + gsl_block_float float + gsl_block_long_double long double + gsl_block_int int + gsl_block_uint unsigned int + gsl_block_long long + gsl_block_ulong unsigned long + gsl_block_short short + gsl_block_ushort unsigned short + gsl_block_char char + gsl_block_uchar unsigned char + gsl_block_complex complex double + gsl_block_complex_float complex float + gsl_block_complex_long_double complex long double + +Corresponding types exist for the `gsl_vector' and `gsl_matrix' +functions. + + +File: gsl-ref.info, Node: Blocks, Next: Vectors, Prev: Data types, Up: Vectors and Matrices + +8.2 Blocks +========== + +For consistency all memory is allocated through a `gsl_block' +structure. The structure contains two components, the size of an area +of memory and a pointer to the memory. The `gsl_block' structure looks +like this, + + typedef struct + { + size_t size; + double * data; + } gsl_block; + +Vectors and matrices are made by "slicing" an underlying block. A slice +is a set of elements formed from an initial offset and a combination of +indices and step-sizes. In the case of a matrix the step-size for the +column index represents the row-length. The step-size for a vector is +known as the "stride". + + The functions for allocating and deallocating blocks are defined in +`gsl_block.h' + +* Menu: + +* Block allocation:: +* Reading and writing blocks:: +* Example programs for blocks:: + + +File: gsl-ref.info, Node: Block allocation, Next: Reading and writing blocks, Up: Blocks + +8.2.1 Block allocation +---------------------- + +The functions for allocating memory to a block follow the style of +`malloc' and `free'. In addition they also perform their own error +checking. If there is insufficient memory available to allocate a +block then the functions call the GSL error handler (with an error +number of `GSL_ENOMEM') in addition to returning a null pointer. Thus +if you use the library error handler to abort your program then it +isn't necessary to check every `alloc'. + + -- Function: gsl_block * gsl_block_alloc (size_t N) + This function allocates memory for a block of N double-precision + elements, returning a pointer to the block struct. The block is + not initialized and so the values of its elements are undefined. + Use the function `gsl_block_calloc' if you want to ensure that all + the elements are initialized to zero. + + A null pointer is returned if insufficient memory is available to + create the block. + + -- Function: gsl_block * gsl_block_calloc (size_t N) + This function allocates memory for a block and initializes all the + elements of the block to zero. + + -- Function: void gsl_block_free (gsl_block * B) + This function frees the memory used by a block B previously + allocated with `gsl_block_alloc' or `gsl_block_calloc'. The block + B must be a valid block object (a null pointer is not allowed). + + +File: gsl-ref.info, Node: Reading and writing blocks, Next: Example programs for blocks, Prev: Block allocation, Up: Blocks + +8.2.2 Reading and writing blocks +-------------------------------- + +The library provides functions for reading and writing blocks to a file +as binary data or formatted text. + + -- Function: int gsl_block_fwrite (FILE * STREAM, const gsl_block * B) + This function writes the elements of the block B to the stream + STREAM in binary format. The return value is 0 for success and + `GSL_EFAILED' if there was a problem writing to the file. Since + the data is written in the native binary format it may not be + portable between different architectures. + + -- Function: int gsl_block_fread (FILE * STREAM, gsl_block * B) + This function reads into the block B from the open stream STREAM + in binary format. The block B must be preallocated with the + correct length since the function uses the size of B to determine + how many bytes to read. The return value is 0 for success and + `GSL_EFAILED' if there was a problem reading from the file. The + data is assumed to have been written in the native binary format + on the same architecture. + + -- Function: int gsl_block_fprintf (FILE * STREAM, const gsl_block * + B, const char * FORMAT) + This function writes the elements of the block B line-by-line to + the stream STREAM using the format specifier FORMAT, which should + be one of the `%g', `%e' or `%f' formats for floating point + numbers and `%d' for integers. The function returns 0 for success + and `GSL_EFAILED' if there was a problem writing to the file. + + -- Function: int gsl_block_fscanf (FILE * STREAM, gsl_block * B) + This function reads formatted data from the stream STREAM into the + block B. The block B must be preallocated with the correct length + since the function uses the size of B to determine how many + numbers to read. The function returns 0 for success and + `GSL_EFAILED' if there was a problem reading from the file. + + +File: gsl-ref.info, Node: Example programs for blocks, Prev: Reading and writing blocks, Up: Blocks + +8.2.3 Example programs for blocks +--------------------------------- + +The following program shows how to allocate a block, + + #include <stdio.h> + #include <gsl/gsl_block.h> + + int + main (void) + { + gsl_block * b = gsl_block_alloc (100); + + printf ("length of block = %u\n", b->size); + printf ("block data address = %#x\n", b->data); + + gsl_block_free (b); + return 0; + } + +Here is the output from the program, + + length of block = 100 + block data address = 0x804b0d8 + + +File: gsl-ref.info, Node: Vectors, Next: Matrices, Prev: Blocks, Up: Vectors and Matrices + +8.3 Vectors +=========== + +Vectors are defined by a `gsl_vector' structure which describes a slice +of a block. Different vectors can be created which point to the same +block. A vector slice is a set of equally-spaced elements of an area +of memory. + + The `gsl_vector' structure contains five components, the "size", the +"stride", a pointer to the memory where the elements are stored, DATA, +a pointer to the block owned by the vector, BLOCK, if any, and an +ownership flag, OWNER. The structure is very simple and looks like +this, + + typedef struct + { + size_t size; + size_t stride; + double * data; + gsl_block * block; + int owner; + } gsl_vector; + +The SIZE is simply the number of vector elements. The range of valid +indices runs from 0 to `size-1'. The STRIDE is the step-size from one +element to the next in physical memory, measured in units of the +appropriate datatype. The pointer DATA gives the location of the first +element of the vector in memory. The pointer BLOCK stores the location +of the memory block in which the vector elements are located (if any). +If the vector owns this block then the OWNER field is set to one and +the block will be deallocated when the vector is freed. If the vector +points to a block owned by another object then the OWNER field is zero +and any underlying block will not be deallocated with the vector. + + The functions for allocating and accessing vectors are defined in +`gsl_vector.h' + +* Menu: + +* Vector allocation:: +* Accessing vector elements:: +* Initializing vector elements:: +* Reading and writing vectors:: +* Vector views:: +* Copying vectors:: +* Exchanging elements:: +* Vector operations:: +* Finding maximum and minimum elements of vectors:: +* Vector properties:: +* Example programs for vectors:: + + +File: gsl-ref.info, Node: Vector allocation, Next: Accessing vector elements, Up: Vectors + +8.3.1 Vector allocation +----------------------- + +The functions for allocating memory to a vector follow the style of +`malloc' and `free'. In addition they also perform their own error +checking. If there is insufficient memory available to allocate a +vector then the functions call the GSL error handler (with an error +number of `GSL_ENOMEM') in addition to returning a null pointer. Thus +if you use the library error handler to abort your program then it +isn't necessary to check every `alloc'. + + -- Function: gsl_vector * gsl_vector_alloc (size_t N) + This function creates a vector of length N, returning a pointer to + a newly initialized vector struct. A new block is allocated for the + elements of the vector, and stored in the BLOCK component of the + vector struct. The block is "owned" by the vector, and will be + deallocated when the vector is deallocated. + + -- Function: gsl_vector * gsl_vector_calloc (size_t N) + This function allocates memory for a vector of length N and + initializes all the elements of the vector to zero. + + -- Function: void gsl_vector_free (gsl_vector * V) + This function frees a previously allocated vector V. If the + vector was created using `gsl_vector_alloc' then the block + underlying the vector will also be deallocated. If the vector has + been created from another object then the memory is still owned by + that object and will not be deallocated. The vector V must be a + valid vector object (a null pointer is not allowed). + + +File: gsl-ref.info, Node: Accessing vector elements, Next: Initializing vector elements, Prev: Vector allocation, Up: Vectors + +8.3.2 Accessing vector elements +------------------------------- + +Unlike FORTRAN compilers, C compilers do not usually provide support +for range checking of vectors and matrices. Range checking is +available in the GNU C Compiler bounds-checking extension, but it is not +part of the default installation of GCC. The functions +`gsl_vector_get' and `gsl_vector_set' can perform portable range +checking for you and report an error if you attempt to access elements +outside the allowed range. + + The functions for accessing the elements of a vector or matrix are +defined in `gsl_vector.h' and declared `extern inline' to eliminate +function-call overhead. You must compile your program with the macro +`HAVE_INLINE' defined to use these functions. + + If necessary you can turn off range checking completely without +modifying any source files by recompiling your program with the +preprocessor definition `GSL_RANGE_CHECK_OFF'. Provided your compiler +supports inline functions the effect of turning off range checking is +to replace calls to `gsl_vector_get(v,i)' by `v->data[i*v->stride]' and +calls to `gsl_vector_set(v,i,x)' by `v->data[i*v->stride]=x'. Thus +there should be no performance penalty for using the range checking +functions when range checking is turned off. + + -- Function: double gsl_vector_get (const gsl_vector * V, size_t I) + This function returns the I-th element of a vector V. If I lies + outside the allowed range of 0 to N-1 then the error handler is + invoked and 0 is returned. + + -- Function: void gsl_vector_set (gsl_vector * V, size_t I, double X) + This function sets the value of the I-th element of a vector V to + X. If I lies outside the allowed range of 0 to N-1 then the error + handler is invoked. + + -- Function: double * gsl_vector_ptr (gsl_vector * V, size_t I) + -- Function: const double * gsl_vector_const_ptr (const gsl_vector * + V, size_t I) + These functions return a pointer to the I-th element of a vector + V. If I lies outside the allowed range of 0 to N-1 then the error + handler is invoked and a null pointer is returned. + + +File: gsl-ref.info, Node: Initializing vector elements, Next: Reading and writing vectors, Prev: Accessing vector elements, Up: Vectors + +8.3.3 Initializing vector elements +---------------------------------- + + -- Function: void gsl_vector_set_all (gsl_vector * V, double X) + This function sets all the elements of the vector V to the value X. + + -- Function: void gsl_vector_set_zero (gsl_vector * V) + This function sets all the elements of the vector V to zero. + + -- Function: int gsl_vector_set_basis (gsl_vector * V, size_t I) + This function makes a basis vector by setting all the elements of + the vector V to zero except for the I-th element which is set to + one. + + +File: gsl-ref.info, Node: Reading and writing vectors, Next: Vector views, Prev: Initializing vector elements, Up: Vectors + +8.3.4 Reading and writing vectors +--------------------------------- + +The library provides functions for reading and writing vectors to a file +as binary data or formatted text. + + -- Function: int gsl_vector_fwrite (FILE * STREAM, const gsl_vector * + V) + This function writes the elements of the vector V to the stream + STREAM in binary format. The return value is 0 for success and + `GSL_EFAILED' if there was a problem writing to the file. Since + the data is written in the native binary format it may not be + portable between different architectures. + + -- Function: int gsl_vector_fread (FILE * STREAM, gsl_vector * V) + This function reads into the vector V from the open stream STREAM + in binary format. The vector V must be preallocated with the + correct length since the function uses the size of V to determine + how many bytes to read. The return value is 0 for success and + `GSL_EFAILED' if there was a problem reading from the file. The + data is assumed to have been written in the native binary format + on the same architecture. + + -- Function: int gsl_vector_fprintf (FILE * STREAM, const gsl_vector * + V, const char * FORMAT) + This function writes the elements of the vector V line-by-line to + the stream STREAM using the format specifier FORMAT, which should + be one of the `%g', `%e' or `%f' formats for floating point + numbers and `%d' for integers. The function returns 0 for success + and `GSL_EFAILED' if there was a problem writing to the file. + + -- Function: int gsl_vector_fscanf (FILE * STREAM, gsl_vector * V) + This function reads formatted data from the stream STREAM into the + vector V. The vector V must be preallocated with the correct + length since the function uses the size of V to determine how many + numbers to read. The function returns 0 for success and + `GSL_EFAILED' if there was a problem reading from the file. + + +File: gsl-ref.info, Node: Vector views, Next: Copying vectors, Prev: Reading and writing vectors, Up: Vectors + +8.3.5 Vector views +------------------ + +In addition to creating vectors from slices of blocks it is also +possible to slice vectors and create vector views. For example, a +subvector of another vector can be described with a view, or two views +can be made which provide access to the even and odd elements of a +vector. + + A vector view is a temporary object, stored on the stack, which can +be used to operate on a subset of vector elements. Vector views can be +defined for both constant and non-constant vectors, using separate types +that preserve constness. A vector view has the type `gsl_vector_view' +and a constant vector view has the type `gsl_vector_const_view'. In +both cases the elements of the view can be accessed as a `gsl_vector' +using the `vector' component of the view object. A pointer to a vector +of type `gsl_vector *' or `const gsl_vector *' can be obtained by +taking the address of this component with the `&' operator. + + When using this pointer it is important to ensure that the view +itself remains in scope--the simplest way to do so is by always writing +the pointer as `&'VIEW`.vector', and never storing this value in +another variable. + + -- Function: gsl_vector_view gsl_vector_subvector (gsl_vector * V, + size_t OFFSET, size_t N) + -- Function: gsl_vector_const_view gsl_vector_const_subvector (const + gsl_vector * V, size_t OFFSET, size_t N) + These functions return a vector view of a subvector of another + vector V. The start of the new vector is offset by OFFSET elements + from the start of the original vector. The new vector has N + elements. Mathematically, the I-th element of the new vector V' + is given by, + + v'(i) = v->data[(offset + i)*v->stride] + + where the index I runs from 0 to `n-1'. + + The `data' pointer of the returned vector struct is set to null if + the combined parameters (OFFSET,N) overrun the end of the original + vector. + + The new vector is only a view of the block underlying the original + vector, V. The block containing the elements of V is not owned by + the new vector. When the view goes out of scope the original + vector V and its block will continue to exist. The original + memory can only be deallocated by freeing the original vector. Of + course, the original vector should not be deallocated while the + view is still in use. + + The function `gsl_vector_const_subvector' is equivalent to + `gsl_vector_subvector' but can be used for vectors which are + declared `const'. + + -- Function: gsl_vector_view gsl_vector_subvector_with_stride + (gsl_vector * V, size_t OFFSET, size_t STRIDE, size_t N) + -- Function: gsl_vector_const_view +gsl_vector_const_subvector_with_stride (const gsl_vector * V, size_t + OFFSET, size_t STRIDE, size_t N) + These functions return a vector view of a subvector of another + vector V with an additional stride argument. The subvector is + formed in the same way as for `gsl_vector_subvector' but the new + vector has N elements with a step-size of STRIDE from one element + to the next in the original vector. Mathematically, the I-th + element of the new vector V' is given by, + + v'(i) = v->data[(offset + i*stride)*v->stride] + + where the index I runs from 0 to `n-1'. + + Note that subvector views give direct access to the underlying + elements of the original vector. For example, the following code + will zero the even elements of the vector `v' of length `n', while + leaving the odd elements untouched, + + gsl_vector_view v_even + = gsl_vector_subvector_with_stride (v, 0, 2, n/2); + gsl_vector_set_zero (&v_even.vector); + + A vector view can be passed to any subroutine which takes a vector + argument just as a directly allocated vector would be, using + `&'VIEW`.vector'. For example, the following code computes the + norm of the odd elements of `v' using the BLAS routine DNRM2, + + gsl_vector_view v_odd + = gsl_vector_subvector_with_stride (v, 1, 2, n/2); + double r = gsl_blas_dnrm2 (&v_odd.vector); + + The function `gsl_vector_const_subvector_with_stride' is equivalent + to `gsl_vector_subvector_with_stride' but can be used for vectors + which are declared `const'. + + -- Function: gsl_vector_view gsl_vector_complex_real + (gsl_vector_complex * V) + -- Function: gsl_vector_const_view gsl_vector_complex_const_real + (const gsl_vector_complex * V) + These functions return a vector view of the real parts of the + complex vector V. + + The function `gsl_vector_complex_const_real' is equivalent to + `gsl_vector_complex_real' but can be used for vectors which are + declared `const'. + + -- Function: gsl_vector_view gsl_vector_complex_imag + (gsl_vector_complex * V) + -- Function: gsl_vector_const_view gsl_vector_complex_const_imag + (const gsl_vector_complex * V) + These functions return a vector view of the imaginary parts of the + complex vector V. + + The function `gsl_vector_complex_const_imag' is equivalent to + `gsl_vector_complex_imag' but can be used for vectors which are + declared `const'. + + -- Function: gsl_vector_view gsl_vector_view_array (double * BASE, + size_t N) + -- Function: gsl_vector_const_view gsl_vector_const_view_array (const + double * BASE, size_t N) + These functions return a vector view of an array. The start of + the new vector is given by BASE and has N elements. + Mathematically, the I-th element of the new vector V' is given by, + + v'(i) = base[i] + + where the index I runs from 0 to `n-1'. + + The array containing the elements of V is not owned by the new + vector view. When the view goes out of scope the original array + will continue to exist. The original memory can only be + deallocated by freeing the original pointer BASE. Of course, the + original array should not be deallocated while the view is still + in use. + + The function `gsl_vector_const_view_array' is equivalent to + `gsl_vector_view_array' but can be used for arrays which are + declared `const'. + + -- Function: gsl_vector_view gsl_vector_view_array_with_stride (double + * BASE, size_t STRIDE, size_t N) + -- Function: gsl_vector_const_view +gsl_vector_const_view_array_with_stride (const double * BASE, size_t + STRIDE, size_t N) + These functions return a vector view of an array BASE with an + additional stride argument. The subvector is formed in the same + way as for `gsl_vector_view_array' but the new vector has N + elements with a step-size of STRIDE from one element to the next + in the original array. Mathematically, the I-th element of the new + vector V' is given by, + + v'(i) = base[i*stride] + + where the index I runs from 0 to `n-1'. + + Note that the view gives direct access to the underlying elements + of the original array. A vector view can be passed to any + subroutine which takes a vector argument just as a directly + allocated vector would be, using `&'VIEW`.vector'. + + The function `gsl_vector_const_view_array_with_stride' is + equivalent to `gsl_vector_view_array_with_stride' but can be used + for arrays which are declared `const'. + + +File: gsl-ref.info, Node: Copying vectors, Next: Exchanging elements, Prev: Vector views, Up: Vectors + +8.3.6 Copying vectors +--------------------- + +Common operations on vectors such as addition and multiplication are +available in the BLAS part of the library (*note BLAS Support::). +However, it is useful to have a small number of utility functions which +do not require the full BLAS code. The following functions fall into +this category. + + -- Function: int gsl_vector_memcpy (gsl_vector * DEST, const + gsl_vector * SRC) + This function copies the elements of the vector SRC into the + vector DEST. The two vectors must have the same length. + + -- Function: int gsl_vector_swap (gsl_vector * V, gsl_vector * W) + This function exchanges the elements of the vectors V and W by + copying. The two vectors must have the same length. + + +File: gsl-ref.info, Node: Exchanging elements, Next: Vector operations, Prev: Copying vectors, Up: Vectors + +8.3.7 Exchanging elements +------------------------- + +The following function can be used to exchange, or permute, the elements +of a vector. + + -- Function: int gsl_vector_swap_elements (gsl_vector * V, size_t I, + size_t J) + This function exchanges the I-th and J-th elements of the vector V + in-place. + + -- Function: int gsl_vector_reverse (gsl_vector * V) + This function reverses the order of the elements of the vector V. + + +File: gsl-ref.info, Node: Vector operations, Next: Finding maximum and minimum elements of vectors, Prev: Exchanging elements, Up: Vectors + +8.3.8 Vector operations +----------------------- + +The following operations are only defined for real vectors. + + -- Function: int gsl_vector_add (gsl_vector * A, const gsl_vector * B) + This function adds the elements of vector B to the elements of + vector A, a'_i = a_i + b_i. The two vectors must have the same + length. + + -- Function: int gsl_vector_sub (gsl_vector * A, const gsl_vector * B) + This function subtracts the elements of vector B from the elements + of vector A, a'_i = a_i - b_i. The two vectors must have the same + length. + + -- Function: int gsl_vector_mul (gsl_vector * A, const gsl_vector * B) + This function multiplies the elements of vector A by the elements + of vector B, a'_i = a_i * b_i. The two vectors must have the same + length. + + -- Function: int gsl_vector_div (gsl_vector * A, const gsl_vector * B) + This function divides the elements of vector A by the elements of + vector B, a'_i = a_i / b_i. The two vectors must have the same + length. + + -- Function: int gsl_vector_scale (gsl_vector * A, const double X) + This function multiplies the elements of vector A by the constant + factor X, a'_i = x a_i. + + -- Function: int gsl_vector_add_constant (gsl_vector * A, const double + X) + This function adds the constant value X to the elements of the + vector A, a'_i = a_i + x. + + +File: gsl-ref.info, Node: Finding maximum and minimum elements of vectors, Next: Vector properties, Prev: Vector operations, Up: Vectors + +8.3.9 Finding maximum and minimum elements of vectors +----------------------------------------------------- + + -- Function: double gsl_vector_max (const gsl_vector * V) + This function returns the maximum value in the vector V. + + -- Function: double gsl_vector_min (const gsl_vector * V) + This function returns the minimum value in the vector V. + + -- Function: void gsl_vector_minmax (const gsl_vector * V, double * + MIN_OUT, double * MAX_OUT) + This function returns the minimum and maximum values in the vector + V, storing them in MIN_OUT and MAX_OUT. + + -- Function: size_t gsl_vector_max_index (const gsl_vector * V) + This function returns the index of the maximum value in the vector + V. When there are several equal maximum elements then the lowest + index is returned. + + -- Function: size_t gsl_vector_min_index (const gsl_vector * V) + This function returns the index of the minimum value in the vector + V. When there are several equal minimum elements then the lowest + index is returned. + + -- Function: void gsl_vector_minmax_index (const gsl_vector * V, + size_t * IMIN, size_t * IMAX) + This function returns the indices of the minimum and maximum + values in the vector V, storing them in IMIN and IMAX. When there + are several equal minimum or maximum elements then the lowest + indices are returned. + + +File: gsl-ref.info, Node: Vector properties, Next: Example programs for vectors, Prev: Finding maximum and minimum elements of vectors, Up: Vectors + +8.3.10 Vector properties +------------------------ + + -- Function: int gsl_vector_isnull (const gsl_vector * V) + -- Function: int gsl_vector_ispos (const gsl_vector * V) + -- Function: int gsl_vector_isneg (const gsl_vector * V) + These functions return 1 if all the elements of the vector V are + zero, strictly positive, or strictly negative respectively, and 0 + otherwise. To test for a non-negative vector, use the expression + `!gsl_vector_isneg(v)'. + + +File: gsl-ref.info, Node: Example programs for vectors, Prev: Vector properties, Up: Vectors + +8.3.11 Example programs for vectors +----------------------------------- + +This program shows how to allocate, initialize and read from a vector +using the functions `gsl_vector_alloc', `gsl_vector_set' and +`gsl_vector_get'. + + #include <stdio.h> + #include <gsl/gsl_vector.h> + + int + main (void) + { + int i; + gsl_vector * v = gsl_vector_alloc (3); + + for (i = 0; i < 3; i++) + { + gsl_vector_set (v, i, 1.23 + i); + } + + for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */ + { + printf ("v_%d = %g\n", i, gsl_vector_get (v, i)); + } + + gsl_vector_free (v); + return 0; + } + +Here is the output from the program. The final loop attempts to read +outside the range of the vector `v', and the error is trapped by the +range-checking code in `gsl_vector_get'. + + $ ./a.out + v_0 = 1.23 + v_1 = 2.23 + v_2 = 3.23 + gsl: vector_source.c:12: ERROR: index out of range + Default GSL error handler invoked. + Aborted (core dumped) + +The next program shows how to write a vector to a file. + + #include <stdio.h> + #include <gsl/gsl_vector.h> + + int + main (void) + { + int i; + gsl_vector * v = gsl_vector_alloc (100); + + for (i = 0; i < 100; i++) + { + gsl_vector_set (v, i, 1.23 + i); + } + + { + FILE * f = fopen ("test.dat", "w"); + gsl_vector_fprintf (f, v, "%.5g"); + fclose (f); + } + + gsl_vector_free (v); + return 0; + } + +After running this program the file `test.dat' should contain the +elements of `v', written using the format specifier `%.5g'. The vector +could then be read back in using the function `gsl_vector_fscanf (f, +v)' as follows: + + #include <stdio.h> + #include <gsl/gsl_vector.h> + + int + main (void) + { + int i; + gsl_vector * v = gsl_vector_alloc (10); + + { + FILE * f = fopen ("test.dat", "r"); + gsl_vector_fscanf (f, v); + fclose (f); + } + + for (i = 0; i < 10; i++) + { + printf ("%g\n", gsl_vector_get(v, i)); + } + + gsl_vector_free (v); + return 0; + } + + +File: gsl-ref.info, Node: Matrices, Next: Vector and Matrix References and Further Reading, Prev: Vectors, Up: Vectors and Matrices + +8.4 Matrices +============ + +Matrices are defined by a `gsl_matrix' structure which describes a +generalized slice of a block. Like a vector it represents a set of +elements in an area of memory, but uses two indices instead of one. + + The `gsl_matrix' structure contains six components, the two +dimensions of the matrix, a physical dimension, a pointer to the memory +where the elements of the matrix are stored, DATA, a pointer to the +block owned by the matrix BLOCK, if any, and an ownership flag, OWNER. +The physical dimension determines the memory layout and can differ from +the matrix dimension to allow the use of submatrices. The `gsl_matrix' +structure is very simple and looks like this, + + typedef struct + { + size_t size1; + size_t size2; + size_t tda; + double * data; + gsl_block * block; + int owner; + } gsl_matrix; + +Matrices are stored in row-major order, meaning that each row of +elements forms a contiguous block in memory. This is the standard +"C-language ordering" of two-dimensional arrays. Note that FORTRAN +stores arrays in column-major order. The number of rows is SIZE1. The +range of valid row indices runs from 0 to `size1-1'. Similarly SIZE2 +is the number of columns. The range of valid column indices runs from +0 to `size2-1'. The physical row dimension TDA, or "trailing +dimension", specifies the size of a row of the matrix as laid out in +memory. + + For example, in the following matrix SIZE1 is 3, SIZE2 is 4, and TDA +is 8. The physical memory layout of the matrix begins in the top left +hand-corner and proceeds from left to right along each row in turn. + + 00 01 02 03 XX XX XX XX + 10 11 12 13 XX XX XX XX + 20 21 22 23 XX XX XX XX + +Each unused memory location is represented by "`XX'". The pointer DATA +gives the location of the first element of the matrix in memory. The +pointer BLOCK stores the location of the memory block in which the +elements of the matrix are located (if any). If the matrix owns this +block then the OWNER field is set to one and the block will be +deallocated when the matrix is freed. If the matrix is only a slice of +a block owned by another object then the OWNER field is zero and any +underlying block will not be freed. + + The functions for allocating and accessing matrices are defined in +`gsl_matrix.h' + +* Menu: + +* Matrix allocation:: +* Accessing matrix elements:: +* Initializing matrix elements:: +* Reading and writing matrices:: +* Matrix views:: +* Creating row and column views:: +* Copying matrices:: +* Copying rows and columns:: +* Exchanging rows and columns:: +* Matrix operations:: +* Finding maximum and minimum elements of matrices:: +* Matrix properties:: +* Example programs for matrices:: + + +File: gsl-ref.info, Node: Matrix allocation, Next: Accessing matrix elements, Up: Matrices + +8.4.1 Matrix allocation +----------------------- + +The functions for allocating memory to a matrix follow the style of +`malloc' and `free'. They also perform their own error checking. If +there is insufficient memory available to allocate a vector then the +functions call the GSL error handler (with an error number of +`GSL_ENOMEM') in addition to returning a null pointer. Thus if you use +the library error handler to abort your program then it isn't necessary +to check every `alloc'. + + -- Function: gsl_matrix * gsl_matrix_alloc (size_t N1, size_t N2) + This function creates a matrix of size N1 rows by N2 columns, + returning a pointer to a newly initialized matrix struct. A new + block is allocated for the elements of the matrix, and stored in + the BLOCK component of the matrix struct. The block is "owned" by + the matrix, and will be deallocated when the matrix is deallocated. + + -- Function: gsl_matrix * gsl_matrix_calloc (size_t N1, size_t N2) + This function allocates memory for a matrix of size N1 rows by N2 + columns and initializes all the elements of the matrix to zero. + + -- Function: void gsl_matrix_free (gsl_matrix * M) + This function frees a previously allocated matrix M. If the + matrix was created using `gsl_matrix_alloc' then the block + underlying the matrix will also be deallocated. If the matrix has + been created from another object then the memory is still owned by + that object and will not be deallocated. The matrix M must be a + valid matrix object (a null pointer is not allowed). + + +File: gsl-ref.info, Node: Accessing matrix elements, Next: Initializing matrix elements, Prev: Matrix allocation, Up: Matrices + +8.4.2 Accessing matrix elements +------------------------------- + +The functions for accessing the elements of a matrix use the same range +checking system as vectors. You can turn off range checking by +recompiling your program with the preprocessor definition +`GSL_RANGE_CHECK_OFF'. + + The elements of the matrix are stored in "C-order", where the second +index moves continuously through memory. More precisely, the element +accessed by the function `gsl_matrix_get(m,i,j)' and +`gsl_matrix_set(m,i,j,x)' is + + m->data[i * m->tda + j] + +where TDA is the physical row-length of the matrix. + + -- Function: double gsl_matrix_get (const gsl_matrix * M, size_t I, + size_t J) + This function returns the (i,j)-th element of a matrix M. If I or + J lie outside the allowed range of 0 to N1-1 and 0 to N2-1 then + the error handler is invoked and 0 is returned. + + -- Function: void gsl_matrix_set (gsl_matrix * M, size_t I, size_t J, + double X) + This function sets the value of the (i,j)-th element of a matrix M + to X. If I or J lies outside the allowed range of 0 to N1-1 and 0 + to N2-1 then the error handler is invoked. + + -- Function: double * gsl_matrix_ptr (gsl_matrix * M, size_t I, size_t + J) + -- Function: const double * gsl_matrix_const_ptr (const gsl_matrix * + M, size_t I, size_t J) + These functions return a pointer to the (i,j)-th element of a + matrix M. If I or J lie outside the allowed range of 0 to N1-1 + and 0 to N2-1 then the error handler is invoked and a null pointer + is returned. + + +File: gsl-ref.info, Node: Initializing matrix elements, Next: Reading and writing matrices, Prev: Accessing matrix elements, Up: Matrices + +8.4.3 Initializing matrix elements +---------------------------------- + + -- Function: void gsl_matrix_set_all (gsl_matrix * M, double X) + This function sets all the elements of the matrix M to the value X. + + -- Function: void gsl_matrix_set_zero (gsl_matrix * M) + This function sets all the elements of the matrix M to zero. + + -- Function: void gsl_matrix_set_identity (gsl_matrix * M) + This function sets the elements of the matrix M to the + corresponding elements of the identity matrix, m(i,j) = + \delta(i,j), i.e. a unit diagonal with all off-diagonal elements + zero. This applies to both square and rectangular matrices. + + +File: gsl-ref.info, Node: Reading and writing matrices, Next: Matrix views, Prev: Initializing matrix elements, Up: Matrices + +8.4.4 Reading and writing matrices +---------------------------------- + +The library provides functions for reading and writing matrices to a +file as binary data or formatted text. + + -- Function: int gsl_matrix_fwrite (FILE * STREAM, const gsl_matrix * + M) + This function writes the elements of the matrix M to the stream + STREAM in binary format. The return value is 0 for success and + `GSL_EFAILED' if there was a problem writing to the file. Since + the data is written in the native binary format it may not be + portable between different architectures. + + -- Function: int gsl_matrix_fread (FILE * STREAM, gsl_matrix * M) + This function reads into the matrix M from the open stream STREAM + in binary format. The matrix M must be preallocated with the + correct dimensions since the function uses the size of M to + determine how many bytes to read. The return value is 0 for + success and `GSL_EFAILED' if there was a problem reading from the + file. The data is assumed to have been written in the native + binary format on the same architecture. + + -- Function: int gsl_matrix_fprintf (FILE * STREAM, const gsl_matrix * + M, const char * FORMAT) + This function writes the elements of the matrix M line-by-line to + the stream STREAM using the format specifier FORMAT, which should + be one of the `%g', `%e' or `%f' formats for floating point + numbers and `%d' for integers. The function returns 0 for success + and `GSL_EFAILED' if there was a problem writing to the file. + + -- Function: int gsl_matrix_fscanf (FILE * STREAM, gsl_matrix * M) + This function reads formatted data from the stream STREAM into the + matrix M. The matrix M must be preallocated with the correct + dimensions since the function uses the size of M to determine how + many numbers to read. The function returns 0 for success and + `GSL_EFAILED' if there was a problem reading from the file. + + +File: gsl-ref.info, Node: Matrix views, Next: Creating row and column views, Prev: Reading and writing matrices, Up: Matrices + +8.4.5 Matrix views +------------------ + +A matrix view is a temporary object, stored on the stack, which can be +used to operate on a subset of matrix elements. Matrix views can be +defined for both constant and non-constant matrices using separate types +that preserve constness. A matrix view has the type `gsl_matrix_view' +and a constant matrix view has the type `gsl_matrix_const_view'. In +both cases the elements of the view can by accessed using the `matrix' +component of the view object. A pointer `gsl_matrix *' or `const +gsl_matrix *' can be obtained by taking the address of the `matrix' +component with the `&' operator. In addition to matrix views it is +also possible to create vector views of a matrix, such as row or column +views. + + -- Function: gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * M, + size_t K1, size_t K2, size_t N1, size_t N2) + -- Function: gsl_matrix_const_view gsl_matrix_const_submatrix (const + gsl_matrix * M, size_t K1, size_t K2, size_t N1, size_t N2) + These functions return a matrix view of a submatrix of the matrix + M. The upper-left element of the submatrix is the element (K1,K2) + of the original matrix. The submatrix has N1 rows and N2 columns. + The physical number of columns in memory given by TDA is + unchanged. Mathematically, the (i,j)-th element of the new matrix + is given by, + + m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j] + + where the index I runs from 0 to `n1-1' and the index J runs from + 0 to `n2-1'. + + The `data' pointer of the returned matrix struct is set to null if + the combined parameters (I,J,N1,N2,TDA) overrun the ends of the + original matrix. + + The new matrix view is only a view of the block underlying the + existing matrix, M. The block containing the elements of M is not + owned by the new matrix view. When the view goes out of scope the + original matrix M and its block will continue to exist. The + original memory can only be deallocated by freeing the original + matrix. Of course, the original matrix should not be deallocated + while the view is still in use. + + The function `gsl_matrix_const_submatrix' is equivalent to + `gsl_matrix_submatrix' but can be used for matrices which are + declared `const'. + + -- Function: gsl_matrix_view gsl_matrix_view_array (double * BASE, + size_t N1, size_t N2) + -- Function: gsl_matrix_const_view gsl_matrix_const_view_array (const + double * BASE, size_t N1, size_t N2) + These functions return a matrix view of the array BASE. The + matrix has N1 rows and N2 columns. The physical number of columns + in memory is also given by N2. Mathematically, the (i,j)-th + element of the new matrix is given by, + + m'(i,j) = base[i*n2 + j] + + where the index I runs from 0 to `n1-1' and the index J runs from + 0 to `n2-1'. + + The new matrix is only a view of the array BASE. When the view + goes out of scope the original array BASE will continue to exist. + The original memory can only be deallocated by freeing the original + array. Of course, the original array should not be deallocated + while the view is still in use. + + The function `gsl_matrix_const_view_array' is equivalent to + `gsl_matrix_view_array' but can be used for matrices which are + declared `const'. + + -- Function: gsl_matrix_view gsl_matrix_view_array_with_tda (double * + BASE, size_t N1, size_t N2, size_t TDA) + -- Function: gsl_matrix_const_view + gsl_matrix_const_view_array_with_tda (const double * BASE, + size_t N1, size_t N2, size_t TDA) + These functions return a matrix view of the array BASE with a + physical number of columns TDA which may differ from the + corresponding dimension of the matrix. The matrix has N1 rows and + N2 columns, and the physical number of columns in memory is given + by TDA. Mathematically, the (i,j)-th element of the new matrix is + given by, + + m'(i,j) = base[i*tda + j] + + where the index I runs from 0 to `n1-1' and the index J runs from + 0 to `n2-1'. + + The new matrix is only a view of the array BASE. When the view + goes out of scope the original array BASE will continue to exist. + The original memory can only be deallocated by freeing the original + array. Of course, the original array should not be deallocated + while the view is still in use. + + The function `gsl_matrix_const_view_array_with_tda' is equivalent + to `gsl_matrix_view_array_with_tda' but can be used for matrices + which are declared `const'. + + -- Function: gsl_matrix_view gsl_matrix_view_vector (gsl_vector * V, + size_t N1, size_t N2) + -- Function: gsl_matrix_const_view gsl_matrix_const_view_vector (const + gsl_vector * V, size_t N1, size_t N2) + These functions return a matrix view of the vector V. The matrix + has N1 rows and N2 columns. The vector must have unit stride. The + physical number of columns in memory is also given by N2. + Mathematically, the (i,j)-th element of the new matrix is given by, + + m'(i,j) = v->data[i*n2 + j] + + where the index I runs from 0 to `n1-1' and the index J runs from + 0 to `n2-1'. + + The new matrix is only a view of the vector V. When the view goes + out of scope the original vector V will continue to exist. The + original memory can only be deallocated by freeing the original + vector. Of course, the original vector should not be deallocated + while the view is still in use. + + The function `gsl_matrix_const_view_vector' is equivalent to + `gsl_matrix_view_vector' but can be used for matrices which are + declared `const'. + + -- Function: gsl_matrix_view gsl_matrix_view_vector_with_tda + (gsl_vector * V, size_t N1, size_t N2, size_t TDA) + -- Function: gsl_matrix_const_view +gsl_matrix_const_view_vector_with_tda (const gsl_vector * V, size_t N1, + size_t N2, size_t TDA) + These functions return a matrix view of the vector V with a + physical number of columns TDA which may differ from the + corresponding matrix dimension. The vector must have unit stride. + The matrix has N1 rows and N2 columns, and the physical number of + columns in memory is given by TDA. Mathematically, the (i,j)-th + element of the new matrix is given by, + + m'(i,j) = v->data[i*tda + j] + + where the index I runs from 0 to `n1-1' and the index J runs from + 0 to `n2-1'. + + The new matrix is only a view of the vector V. When the view goes + out of scope the original vector V will continue to exist. The + original memory can only be deallocated by freeing the original + vector. Of course, the original vector should not be deallocated + while the view is still in use. + + The function `gsl_matrix_const_view_vector_with_tda' is equivalent + to `gsl_matrix_view_vector_with_tda' but can be used for matrices + which are declared `const'. + + +File: gsl-ref.info, Node: Creating row and column views, Next: Copying matrices, Prev: Matrix views, Up: Matrices + +8.4.6 Creating row and column views +----------------------------------- + +In general there are two ways to access an object, by reference or by +copying. The functions described in this section create vector views +which allow access to a row or column of a matrix by reference. +Modifying elements of the view is equivalent to modifying the matrix, +since both the vector view and the matrix point to the same memory +block. + + -- Function: gsl_vector_view gsl_matrix_row (gsl_matrix * M, size_t I) + -- Function: gsl_vector_const_view gsl_matrix_const_row (const + gsl_matrix * M, size_t I) + These functions return a vector view of the I-th row of the matrix + M. The `data' pointer of the new vector is set to null if I is + out of range. + + The function `gsl_vector_const_row' is equivalent to + `gsl_matrix_row' but can be used for matrices which are declared + `const'. + + -- Function: gsl_vector_view gsl_matrix_column (gsl_matrix * M, size_t + J) + -- Function: gsl_vector_const_view gsl_matrix_const_column (const + gsl_matrix * M, size_t J) + These functions return a vector view of the J-th column of the + matrix M. The `data' pointer of the new vector is set to null if + J is out of range. + + The function `gsl_vector_const_column' is equivalent to + `gsl_matrix_column' but can be used for matrices which are declared + `const'. + + -- Function: gsl_vector_view gsl_matrix_diagonal (gsl_matrix * M) + -- Function: gsl_vector_const_view gsl_matrix_const_diagonal (const + gsl_matrix * M) + These functions returns a vector view of the diagonal of the matrix + M. The matrix M is not required to be square. For a rectangular + matrix the length of the diagonal is the same as the smaller + dimension of the matrix. + + The function `gsl_matrix_const_diagonal' is equivalent to + `gsl_matrix_diagonal' but can be used for matrices which are + declared `const'. + + -- Function: gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * M, + size_t K) + -- Function: gsl_vector_const_view gsl_matrix_const_subdiagonal (const + gsl_matrix * M, size_t K) + These functions return a vector view of the K-th subdiagonal of + the matrix M. The matrix M is not required to be square. The + diagonal of the matrix corresponds to k = 0. + + The function `gsl_matrix_const_subdiagonal' is equivalent to + `gsl_matrix_subdiagonal' but can be used for matrices which are + declared `const'. + + -- Function: gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * M, + size_t K) + -- Function: gsl_vector_const_view gsl_matrix_const_superdiagonal + (const gsl_matrix * M, size_t K) + These functions return a vector view of the K-th superdiagonal of + the matrix M. The matrix M is not required to be square. The + diagonal of the matrix corresponds to k = 0. + + The function `gsl_matrix_const_superdiagonal' is equivalent to + `gsl_matrix_superdiagonal' but can be used for matrices which are + declared `const'. + + +File: gsl-ref.info, Node: Copying matrices, Next: Copying rows and columns, Prev: Creating row and column views, Up: Matrices + +8.4.7 Copying matrices +---------------------- + + -- Function: int gsl_matrix_memcpy (gsl_matrix * DEST, const + gsl_matrix * SRC) + This function copies the elements of the matrix SRC into the + matrix DEST. The two matrices must have the same size. + + -- Function: int gsl_matrix_swap (gsl_matrix * M1, gsl_matrix * M2) + This function exchanges the elements of the matrices M1 and M2 by + copying. The two matrices must have the same size. + + +File: gsl-ref.info, Node: Copying rows and columns, Next: Exchanging rows and columns, Prev: Copying matrices, Up: Matrices + +8.4.8 Copying rows and columns +------------------------------ + +The functions described in this section copy a row or column of a matrix +into a vector. This allows the elements of the vector and the matrix to +be modified independently. Note that if the matrix and the vector point +to overlapping regions of memory then the result will be undefined. The +same effect can be achieved with more generality using +`gsl_vector_memcpy' with vector views of rows and columns. + + -- Function: int gsl_matrix_get_row (gsl_vector * V, const gsl_matrix + * M, size_t I) + This function copies the elements of the I-th row of the matrix M + into the vector V. The length of the vector must be the same as + the length of the row. + + -- Function: int gsl_matrix_get_col (gsl_vector * V, const gsl_matrix + * M, size_t J) + This function copies the elements of the J-th column of the matrix + M into the vector V. The length of the vector must be the same as + the length of the column. + + -- Function: int gsl_matrix_set_row (gsl_matrix * M, size_t I, const + gsl_vector * V) + This function copies the elements of the vector V into the I-th + row of the matrix M. The length of the vector must be the same as + the length of the row. + + -- Function: int gsl_matrix_set_col (gsl_matrix * M, size_t J, const + gsl_vector * V) + This function copies the elements of the vector V into the J-th + column of the matrix M. The length of the vector must be the same + as the length of the column. + + +File: gsl-ref.info, Node: Exchanging rows and columns, Next: Matrix operations, Prev: Copying rows and columns, Up: Matrices + +8.4.9 Exchanging rows and columns +--------------------------------- + +The following functions can be used to exchange the rows and columns of +a matrix. + + -- Function: int gsl_matrix_swap_rows (gsl_matrix * M, size_t I, + size_t J) + This function exchanges the I-th and J-th rows of the matrix M + in-place. + + -- Function: int gsl_matrix_swap_columns (gsl_matrix * M, size_t I, + size_t J) + This function exchanges the I-th and J-th columns of the matrix M + in-place. + + -- Function: int gsl_matrix_swap_rowcol (gsl_matrix * M, size_t I, + size_t J) + This function exchanges the I-th row and J-th column of the matrix + M in-place. The matrix must be square for this operation to be + possible. + + -- Function: int gsl_matrix_transpose_memcpy (gsl_matrix * DEST, const + gsl_matrix * SRC) + This function makes the matrix DEST the transpose of the matrix + SRC by copying the elements of SRC into DEST. This function works + for all matrices provided that the dimensions of the matrix DEST + match the transposed dimensions of the matrix SRC. + + -- Function: int gsl_matrix_transpose (gsl_matrix * M) + This function replaces the matrix M by its transpose by copying + the elements of the matrix in-place. The matrix must be square + for this operation to be possible. + + +File: gsl-ref.info, Node: Matrix operations, Next: Finding maximum and minimum elements of matrices, Prev: Exchanging rows and columns, Up: Matrices + +8.4.10 Matrix operations +------------------------ + +The following operations are defined for real and complex matrices. + + -- Function: int gsl_matrix_add (gsl_matrix * A, const gsl_matrix * B) + This function adds the elements of matrix B to the elements of + matrix A, a'(i,j) = a(i,j) + b(i,j). The two matrices must have the + same dimensions. + + -- Function: int gsl_matrix_sub (gsl_matrix * A, const gsl_matrix * B) + This function subtracts the elements of matrix B from the elements + of matrix A, a'(i,j) = a(i,j) - b(i,j). The two matrices must have + the same dimensions. + + -- Function: int gsl_matrix_mul_elements (gsl_matrix * A, const + gsl_matrix * B) + This function multiplies the elements of matrix A by the elements + of matrix B, a'(i,j) = a(i,j) * b(i,j). The two matrices must have + the same dimensions. + + -- Function: int gsl_matrix_div_elements (gsl_matrix * A, const + gsl_matrix * B) + This function divides the elements of matrix A by the elements of + matrix B, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the + same dimensions. + + -- Function: int gsl_matrix_scale (gsl_matrix * A, const double X) + This function multiplies the elements of matrix A by the constant + factor X, a'(i,j) = x a(i,j). + + -- Function: int gsl_matrix_add_constant (gsl_matrix * A, const double + X) + This function adds the constant value X to the elements of the + matrix A, a'(i,j) = a(i,j) + x. + + +File: gsl-ref.info, Node: Finding maximum and minimum elements of matrices, Next: Matrix properties, Prev: Matrix operations, Up: Matrices + +8.4.11 Finding maximum and minimum elements of matrices +------------------------------------------------------- + +The following operations are only defined for real matrices. + + -- Function: double gsl_matrix_max (const gsl_matrix * M) + This function returns the maximum value in the matrix M. + + -- Function: double gsl_matrix_min (const gsl_matrix * M) + This function returns the minimum value in the matrix M. + + -- Function: void gsl_matrix_minmax (const gsl_matrix * M, double * + MIN_OUT, double * MAX_OUT) + This function returns the minimum and maximum values in the matrix + M, storing them in MIN_OUT and MAX_OUT. + + -- Function: void gsl_matrix_max_index (const gsl_matrix * M, size_t * + IMAX, size_t * JMAX) + This function returns the indices of the maximum value in the + matrix M, storing them in IMAX and JMAX. When there are several + equal maximum elements then the first element found is returned, + searching in row-major order. + + -- Function: void gsl_matrix_min_index (const gsl_matrix * M, size_t * + IMIN, size_t * JMIN) + This function returns the indices of the minimum value in the + matrix M, storing them in IMIN and JMIN. When there are several + equal minimum elements then the first element found is returned, + searching in row-major order. + + -- Function: void gsl_matrix_minmax_index (const gsl_matrix * M, + size_t * IMIN, size_t * JMIN, size_t * IMAX, size_t * JMAX) + This function returns the indices of the minimum and maximum + values in the matrix M, storing them in (IMIN,JMIN) and + (IMAX,JMAX). When there are several equal minimum or maximum + elements then the first elements found are returned, searching in + row-major order. + + +File: gsl-ref.info, Node: Matrix properties, Next: Example programs for matrices, Prev: Finding maximum and minimum elements of matrices, Up: Matrices + +8.4.12 Matrix properties +------------------------ + + -- Function: int gsl_matrix_isnull (const gsl_matrix * M) + -- Function: int gsl_matrix_ispos (const gsl_matrix * M) + -- Function: int gsl_matrix_isneg (const gsl_matrix * M) + These functions return 1 if all the elements of the matrix M are + zero, strictly positive, or strictly negative respectively, and 0 + otherwise. To test for a non-negative matrix, use the expression + `!gsl_matrix_isneg(m)'. To test whether a matrix is + positive-definite, use the Cholesky decomposition (*note Cholesky + Decomposition::). + + +File: gsl-ref.info, Node: Example programs for matrices, Prev: Matrix properties, Up: Matrices + +8.4.13 Example programs for matrices +------------------------------------ + +The program below shows how to allocate, initialize and read from a +matrix using the functions `gsl_matrix_alloc', `gsl_matrix_set' and +`gsl_matrix_get'. + + #include <stdio.h> + #include <gsl/gsl_matrix.h> + + int + main (void) + { + int i, j; + gsl_matrix * m = gsl_matrix_alloc (10, 3); + + for (i = 0; i < 10; i++) + for (j = 0; j < 3; j++) + gsl_matrix_set (m, i, j, 0.23 + 100*i + j); + + for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */ + for (j = 0; j < 3; j++) + printf ("m(%d,%d) = %g\n", i, j, + gsl_matrix_get (m, i, j)); + + gsl_matrix_free (m); + + return 0; + } + +Here is the output from the program. The final loop attempts to read +outside the range of the matrix `m', and the error is trapped by the +range-checking code in `gsl_matrix_get'. + + $ ./a.out + m(0,0) = 0.23 + m(0,1) = 1.23 + m(0,2) = 2.23 + m(1,0) = 100.23 + m(1,1) = 101.23 + m(1,2) = 102.23 + ... + m(9,2) = 902.23 + gsl: matrix_source.c:13: ERROR: first index out of range + Default GSL error handler invoked. + Aborted (core dumped) + +The next program shows how to write a matrix to a file. + + #include <stdio.h> + #include <gsl/gsl_matrix.h> + + int + main (void) + { + int i, j, k = 0; + gsl_matrix * m = gsl_matrix_alloc (100, 100); + gsl_matrix * a = gsl_matrix_alloc (100, 100); + + for (i = 0; i < 100; i++) + for (j = 0; j < 100; j++) + gsl_matrix_set (m, i, j, 0.23 + i + j); + + { + FILE * f = fopen ("test.dat", "wb"); + gsl_matrix_fwrite (f, m); + fclose (f); + } + + { + FILE * f = fopen ("test.dat", "rb"); + gsl_matrix_fread (f, a); + fclose (f); + } + + for (i = 0; i < 100; i++) + for (j = 0; j < 100; j++) + { + double mij = gsl_matrix_get (m, i, j); + double aij = gsl_matrix_get (a, i, j); + if (mij != aij) k++; + } + + gsl_matrix_free (m); + gsl_matrix_free (a); + + printf ("differences = %d (should be zero)\n", k); + return (k > 0); + } + +After running this program the file `test.dat' should contain the +elements of `m', written in binary format. The matrix which is read +back in using the function `gsl_matrix_fread' should be exactly equal +to the original matrix. + + The following program demonstrates the use of vector views. The +program computes the column norms of a matrix. + + #include <math.h> + #include <stdio.h> + #include <gsl/gsl_matrix.h> + #include <gsl/gsl_blas.h> + + int + main (void) + { + size_t i,j; + + gsl_matrix *m = gsl_matrix_alloc (10, 10); + + for (i = 0; i < 10; i++) + for (j = 0; j < 10; j++) + gsl_matrix_set (m, i, j, sin (i) + cos (j)); + + for (j = 0; j < 10; j++) + { + gsl_vector_view column = gsl_matrix_column (m, j); + double d; + + d = gsl_blas_dnrm2 (&column.vector); + + printf ("matrix column %d, norm = %g\n", j, d); + } + + gsl_matrix_free (m); + + return 0; + } + +Here is the output of the program, + + $ ./a.out + matrix column 0, norm = 4.31461 + matrix column 1, norm = 3.1205 + matrix column 2, norm = 2.19316 + matrix column 3, norm = 3.26114 + matrix column 4, norm = 2.53416 + matrix column 5, norm = 2.57281 + matrix column 6, norm = 4.20469 + matrix column 7, norm = 3.65202 + matrix column 8, norm = 2.08524 + matrix column 9, norm = 3.07313 + +The results can be confirmed using GNU OCTAVE, + + $ octave + GNU Octave, version 2.0.16.92 + octave> m = sin(0:9)' * ones(1,10) + + ones(10,1) * cos(0:9); + octave> sqrt(sum(m.^2)) + ans = + 4.3146 3.1205 2.1932 3.2611 2.5342 2.5728 + 4.2047 3.6520 2.0852 3.0731 + + +File: gsl-ref.info, Node: Vector and Matrix References and Further Reading, Prev: Matrices, Up: Vectors and Matrices + +8.5 References and Further Reading +================================== + +The block, vector and matrix objects in GSL follow the `valarray' model +of C++. A description of this model can be found in the following +reference, + + B. Stroustrup, `The C++ Programming Language' (3rd Ed), Section + 22.4 Vector Arithmetic. Addison-Wesley 1997, ISBN 0-201-88954-4. + + +File: gsl-ref.info, Node: Permutations, Next: Combinations, Prev: Vectors and Matrices, Up: Top + +9 Permutations +************** + +This chapter describes functions for creating and manipulating +permutations. A permutation p is represented by an array of n integers +in the range 0 to n-1, where each value p_i occurs once and only once. +The application of a permutation p to a vector v yields a new vector v' +where v'_i = v_{p_i}. For example, the array (0,1,3,2) represents a +permutation which exchanges the last two elements of a four element +vector. The corresponding identity permutation is (0,1,2,3). + + Note that the permutations produced by the linear algebra routines +correspond to the exchange of matrix columns, and so should be +considered as applying to row-vectors in the form v' = v P rather than +column-vectors, when permuting the elements of a vector. + + The functions described in this chapter are defined in the header +file `gsl_permutation.h'. + +* Menu: + +* The Permutation struct:: +* Permutation allocation:: +* Accessing permutation elements:: +* Permutation properties:: +* Permutation functions:: +* Applying Permutations:: +* Reading and writing permutations:: +* Permutations in cyclic form:: +* Permutation Examples:: +* Permutation References and Further Reading:: + + +File: gsl-ref.info, Node: The Permutation struct, Next: Permutation allocation, Up: Permutations + +9.1 The Permutation struct +========================== + +A permutation is defined by a structure containing two components, the +size of the permutation and a pointer to the permutation array. The +elements of the permutation array are all of type `size_t'. The +`gsl_permutation' structure looks like this, + + typedef struct + { + size_t size; + size_t * data; + } gsl_permutation; + + + +File: gsl-ref.info, Node: Permutation allocation, Next: Accessing permutation elements, Prev: The Permutation struct, Up: Permutations + +9.2 Permutation allocation +========================== + + -- Function: gsl_permutation * gsl_permutation_alloc (size_t N) + This function allocates memory for a new permutation of size N. + The permutation is not initialized and its elements are undefined. + Use the function `gsl_permutation_calloc' if you want to create a + permutation which is initialized to the identity. A null pointer is + returned if insufficient memory is available to create the + permutation. + + -- Function: gsl_permutation * gsl_permutation_calloc (size_t N) + This function allocates memory for a new permutation of size N and + initializes it to the identity. A null pointer is returned if + insufficient memory is available to create the permutation. + + -- Function: void gsl_permutation_init (gsl_permutation * P) + This function initializes the permutation P to the identity, i.e. + (0,1,2,...,n-1). + + -- Function: void gsl_permutation_free (gsl_permutation * P) + This function frees all the memory used by the permutation P. + + -- Function: int gsl_permutation_memcpy (gsl_permutation * DEST, const + gsl_permutation * SRC) + This function copies the elements of the permutation SRC into the + permutation DEST. The two permutations must have the same size. + + +File: gsl-ref.info, Node: Accessing permutation elements, Next: Permutation properties, Prev: Permutation allocation, Up: Permutations + +9.3 Accessing permutation elements +================================== + +The following functions can be used to access and manipulate +permutations. + + -- Function: size_t gsl_permutation_get (const gsl_permutation * P, + const size_t I) + This function returns the value of the I-th element of the + permutation P. If I lies outside the allowed range of 0 to N-1 + then the error handler is invoked and 0 is returned. + + -- Function: int gsl_permutation_swap (gsl_permutation * P, const + size_t I, const size_t J) + This function exchanges the I-th and J-th elements of the + permutation P. + + +File: gsl-ref.info, Node: Permutation properties, Next: Permutation functions, Prev: Accessing permutation elements, Up: Permutations + +9.4 Permutation properties +========================== + + -- Function: size_t gsl_permutation_size (const gsl_permutation * P) + This function returns the size of the permutation P. + + -- Function: size_t * gsl_permutation_data (const gsl_permutation * P) + This function returns a pointer to the array of elements in the + permutation P. + + -- Function: int gsl_permutation_valid (gsl_permutation * P) + This function checks that the permutation P is valid. The N + elements should contain each of the numbers 0 to N-1 once and only + once. + + +File: gsl-ref.info, Node: Permutation functions, Next: Applying Permutations, Prev: Permutation properties, Up: Permutations + +9.5 Permutation functions +========================= + + -- Function: void gsl_permutation_reverse (gsl_permutation * P) + This function reverses the elements of the permutation P. + + -- Function: int gsl_permutation_inverse (gsl_permutation * INV, const + gsl_permutation * P) + This function computes the inverse of the permutation P, storing + the result in INV. + + -- Function: int gsl_permutation_next (gsl_permutation * P) + This function advances the permutation P to the next permutation + in lexicographic order and returns `GSL_SUCCESS'. If no further + permutations are available it returns `GSL_FAILURE' and leaves P + unmodified. Starting with the identity permutation and repeatedly + applying this function will iterate through all possible + permutations of a given order. + + -- Function: int gsl_permutation_prev (gsl_permutation * P) + This function steps backwards from the permutation P to the + previous permutation in lexicographic order, returning + `GSL_SUCCESS'. If no previous permutation is available it returns + `GSL_FAILURE' and leaves P unmodified. + + +File: gsl-ref.info, Node: Applying Permutations, Next: Reading and writing permutations, Prev: Permutation functions, Up: Permutations + +9.6 Applying Permutations +========================= + + -- Function: int gsl_permute (const size_t * P, double * DATA, size_t + STRIDE, size_t N) + This function applies the permutation P to the array DATA of size + N with stride STRIDE. + + -- Function: int gsl_permute_inverse (const size_t * P, double * DATA, + size_t STRIDE, size_t N) + This function applies the inverse of the permutation P to the + array DATA of size N with stride STRIDE. + + -- Function: int gsl_permute_vector (const gsl_permutation * P, + gsl_vector * V) + This function applies the permutation P to the elements of the + vector V, considered as a row-vector acted on by a permutation + matrix from the right, v' = v P. The j-th column of the + permutation matrix P is given by the p_j-th column of the identity + matrix. The permutation P and the vector V must have the same + length. + + -- Function: int gsl_permute_vector_inverse (const gsl_permutation * + P, gsl_vector * V) + This function applies the inverse of the permutation P to the + elements of the vector V, considered as a row-vector acted on by + an inverse permutation matrix from the right, v' = v P^T. Note + that for permutation matrices the inverse is the same as the + transpose. The j-th column of the permutation matrix P is given by + the p_j-th column of the identity matrix. The permutation P and + the vector V must have the same length. + + -- Function: int gsl_permutation_mul (gsl_permutation * P, const + gsl_permutation * PA, const gsl_permutation * PB) + This function combines the two permutations PA and PB into a + single permutation P, where p = pa . pb. The permutation P is + equivalent to applying pb first and then PA. + + +File: gsl-ref.info, Node: Reading and writing permutations, Next: Permutations in cyclic form, Prev: Applying Permutations, Up: Permutations + +9.7 Reading and writing permutations +==================================== + +The library provides functions for reading and writing permutations to a +file as binary data or formatted text. + + -- Function: int gsl_permutation_fwrite (FILE * STREAM, const + gsl_permutation * P) + This function writes the elements of the permutation P to the + stream STREAM in binary format. The function returns + `GSL_EFAILED' if there was a problem writing to the file. Since + the data is written in the native binary format it may not be + portable between different architectures. + + -- Function: int gsl_permutation_fread (FILE * STREAM, gsl_permutation + * P) + This function reads into the permutation P from the open stream + STREAM in binary format. The permutation P must be preallocated + with the correct length since the function uses the size of P to + determine how many bytes to read. The function returns + `GSL_EFAILED' if there was a problem reading from the file. The + data is assumed to have been written in the native binary format + on the same architecture. + + -- Function: int gsl_permutation_fprintf (FILE * STREAM, const + gsl_permutation * P, const char * FORMAT) + This function writes the elements of the permutation P + line-by-line to the stream STREAM using the format specifier + FORMAT, which should be suitable for a type of SIZE_T. On a GNU + system the type modifier `Z' represents `size_t', so `"%Zu\n"' is + a suitable format. The function returns `GSL_EFAILED' if there + was a problem writing to the file. + + -- Function: int gsl_permutation_fscanf (FILE * STREAM, + gsl_permutation * P) + This function reads formatted data from the stream STREAM into the + permutation P. The permutation P must be preallocated with the + correct length since the function uses the size of P to determine + how many numbers to read. The function returns `GSL_EFAILED' if + there was a problem reading from the file. + + +File: gsl-ref.info, Node: Permutations in cyclic form, Next: Permutation Examples, Prev: Reading and writing permutations, Up: Permutations + +9.8 Permutations in cyclic form +=============================== + +A permutation can be represented in both "linear" and "cyclic" +notations. The functions described in this section convert between the +two forms. The linear notation is an index mapping, and has already +been described above. The cyclic notation expresses a permutation as a +series of circular rearrangements of groups of elements, or "cycles". + + For example, under the cycle (1 2 3), 1 is replaced by 2, 2 is +replaced by 3 and 3 is replaced by 1 in a circular fashion. Cycles of +different sets of elements can be combined independently, for example +(1 2 3) (4 5) combines the cycle (1 2 3) with the cycle (4 5), which is +an exchange of elements 4 and 5. A cycle of length one represents an +element which is unchanged by the permutation and is referred to as a +"singleton". + + It can be shown that every permutation can be decomposed into +combinations of cycles. The decomposition is not unique, but can always +be rearranged into a standard "canonical form" by a reordering of +elements. The library uses the canonical form defined in Knuth's `Art +of Computer Programming' (Vol 1, 3rd Ed, 1997) Section 1.3.3, p.178. + + The procedure for obtaining the canonical form given by Knuth is, + + 1. Write all singleton cycles explicitly + + 2. Within each cycle, put the smallest number first + + 3. Order the cycles in decreasing order of the first number in the + cycle. + +For example, the linear representation (2 4 3 0 1) is represented as (1 +4) (0 2 3) in canonical form. The permutation corresponds to an +exchange of elements 1 and 4, and rotation of elements 0, 2 and 3. + + The important property of the canonical form is that it can be +reconstructed from the contents of each cycle without the brackets. In +addition, by removing the brackets it can be considered as a linear +representation of a different permutation. In the example given above +the permutation (2 4 3 0 1) would become (1 4 0 2 3). This mapping has +many applications in the theory of permutations. + + -- Function: int gsl_permutation_linear_to_canonical (gsl_permutation + * Q, const gsl_permutation * P) + This function computes the canonical form of the permutation P and + stores it in the output argument Q. + + -- Function: int gsl_permutation_canonical_to_linear (gsl_permutation + * P, const gsl_permutation * Q) + This function converts a permutation Q in canonical form back into + linear form storing it in the output argument P. + + -- Function: size_t gsl_permutation_inversions (const gsl_permutation + * P) + This function counts the number of inversions in the permutation + P. An inversion is any pair of elements that are not in order. + For example, the permutation 2031 has three inversions, + corresponding to the pairs (2,0) (2,1) and (3,1). The identity + permutation has no inversions. + + -- Function: size_t gsl_permutation_linear_cycles (const + gsl_permutation * P) + This function counts the number of cycles in the permutation P, + given in linear form. + + -- Function: size_t gsl_permutation_canonical_cycles (const + gsl_permutation * Q) + This function counts the number of cycles in the permutation Q, + given in canonical form. + + +File: gsl-ref.info, Node: Permutation Examples, Next: Permutation References and Further Reading, Prev: Permutations in cyclic form, Up: Permutations + +9.9 Examples +============ + +The example program below creates a random permutation (by shuffling the +elements of the identity) and finds its inverse. + + #include <stdio.h> + #include <gsl/gsl_rng.h> + #include <gsl/gsl_randist.h> + #include <gsl/gsl_permutation.h> + + int + main (void) + { + const size_t N = 10; + const gsl_rng_type * T; + gsl_rng * r; + + gsl_permutation * p = gsl_permutation_alloc (N); + gsl_permutation * q = gsl_permutation_alloc (N); + + gsl_rng_env_setup(); + T = gsl_rng_default; + r = gsl_rng_alloc (T); + + printf ("initial permutation:"); + gsl_permutation_init (p); + gsl_permutation_fprintf (stdout, p, " %u"); + printf ("\n"); + + printf (" random permutation:"); + gsl_ran_shuffle (r, p->data, N, sizeof(size_t)); + gsl_permutation_fprintf (stdout, p, " %u"); + printf ("\n"); + + printf ("inverse permutation:"); + gsl_permutation_inverse (q, p); + gsl_permutation_fprintf (stdout, q, " %u"); + printf ("\n"); + + gsl_permutation_free (p); + gsl_permutation_free (q); + gsl_rng_free (r); + + return 0; + } + +Here is the output from the program, + + $ ./a.out + initial permutation: 0 1 2 3 4 5 6 7 8 9 + random permutation: 1 3 5 2 7 6 0 4 9 8 + inverse permutation: 6 0 3 1 7 2 5 4 9 8 + +The random permutation `p[i]' and its inverse `q[i]' are related +through the identity `p[q[i]] = i', which can be verified from the +output. + + The next example program steps forwards through all possible third +order permutations, starting from the identity, + + #include <stdio.h> + #include <gsl/gsl_permutation.h> + + int + main (void) + { + gsl_permutation * p = gsl_permutation_alloc (3); + + gsl_permutation_init (p); + + do + { + gsl_permutation_fprintf (stdout, p, " %u"); + printf ("\n"); + } + while (gsl_permutation_next(p) == GSL_SUCCESS); + + gsl_permutation_free (p); + + return 0; + } + +Here is the output from the program, + + $ ./a.out + 0 1 2 + 0 2 1 + 1 0 2 + 1 2 0 + 2 0 1 + 2 1 0 + +The permutations are generated in lexicographic order. To reverse the +sequence, begin with the final permutation (which is the reverse of the +identity) and replace `gsl_permutation_next' with +`gsl_permutation_prev'. + + +File: gsl-ref.info, Node: Permutation References and Further Reading, Prev: Permutation Examples, Up: Permutations + +9.10 References and Further Reading +=================================== + +The subject of permutations is covered extensively in Knuth's `Sorting +and Searching', + + Donald E. Knuth, `The Art of Computer Programming: Sorting and + Searching' (Vol 3, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850. + +For the definition of the "canonical form" see, + + Donald E. Knuth, `The Art of Computer Programming: Fundamental + Algorithms' (Vol 1, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850. + Section 1.3.3, `An Unusual Correspondence', p.178-179. + + +File: gsl-ref.info, Node: Combinations, Next: Sorting, Prev: Permutations, Up: Top + +10 Combinations +*************** + +This chapter describes functions for creating and manipulating +combinations. A combination c is represented by an array of k integers +in the range 0 to n-1, where each value c_i occurs at most once. The +combination c corresponds to indices of k elements chosen from an n +element vector. Combinations are useful for iterating over all +k-element subsets of a set. + + The functions described in this chapter are defined in the header +file `gsl_combination.h'. + +* Menu: + +* The Combination struct:: +* Combination allocation:: +* Accessing combination elements:: +* Combination properties:: +* Combination functions:: +* Reading and writing combinations:: +* Combination Examples:: +* Combination References and Further Reading:: + + +File: gsl-ref.info, Node: The Combination struct, Next: Combination allocation, Up: Combinations + +10.1 The Combination struct +=========================== + +A combination is defined by a structure containing three components, the +values of n and k, and a pointer to the combination array. The +elements of the combination array are all of type `size_t', and are +stored in increasing order. The `gsl_combination' structure looks like +this, + + typedef struct + { + size_t n; + size_t k; + size_t *data; + } gsl_combination; + + + +File: gsl-ref.info, Node: Combination allocation, Next: Accessing combination elements, Prev: The Combination struct, Up: Combinations + +10.2 Combination allocation +=========================== + + -- Function: gsl_combination * gsl_combination_alloc (size_t N, size_t + K) + This function allocates memory for a new combination with + parameters N, K. The combination is not initialized and its + elements are undefined. Use the function `gsl_combination_calloc' + if you want to create a combination which is initialized to the + lexicographically first combination. A null pointer is returned if + insufficient memory is available to create the combination. + + -- Function: gsl_combination * gsl_combination_calloc (size_t N, + size_t K) + This function allocates memory for a new combination with + parameters N, K and initializes it to the lexicographically first + combination. A null pointer is returned if insufficient memory is + available to create the combination. + + -- Function: void gsl_combination_init_first (gsl_combination * C) + This function initializes the combination C to the + lexicographically first combination, i.e. (0,1,2,...,k-1). + + -- Function: void gsl_combination_init_last (gsl_combination * C) + This function initializes the combination C to the + lexicographically last combination, i.e. (n-k,n-k+1,...,n-1). + + -- Function: void gsl_combination_free (gsl_combination * C) + This function frees all the memory used by the combination C. + + -- Function: int gsl_combination_memcpy (gsl_combination * DEST, const + gsl_combination * SRC) + This function copies the elements of the combination SRC into the + combination DEST. The two combinations must have the same size. + + +File: gsl-ref.info, Node: Accessing combination elements, Next: Combination properties, Prev: Combination allocation, Up: Combinations + +10.3 Accessing combination elements +=================================== + +The following function can be used to access the elements of a +combination. + + -- Function: size_t gsl_combination_get (const gsl_combination * C, + const size_t I) + This function returns the value of the I-th element of the + combination C. If I lies outside the allowed range of 0 to K-1 + then the error handler is invoked and 0 is returned. + + +File: gsl-ref.info, Node: Combination properties, Next: Combination functions, Prev: Accessing combination elements, Up: Combinations + +10.4 Combination properties +=========================== + + -- Function: size_t gsl_combination_n (const gsl_combination * C) + This function returns the range (n) of the combination C. + + -- Function: size_t gsl_combination_k (const gsl_combination * C) + This function returns the number of elements (k) in the + combination C. + + -- Function: size_t * gsl_combination_data (const gsl_combination * C) + This function returns a pointer to the array of elements in the + combination C. + + -- Function: int gsl_combination_valid (gsl_combination * C) + This function checks that the combination C is valid. The K + elements should lie in the range 0 to N-1, with each value + occurring once at most and in increasing order. + + +File: gsl-ref.info, Node: Combination functions, Next: Reading and writing combinations, Prev: Combination properties, Up: Combinations + +10.5 Combination functions +========================== + + -- Function: int gsl_combination_next (gsl_combination * C) + This function advances the combination C to the next combination + in lexicographic order and returns `GSL_SUCCESS'. If no further + combinations are available it returns `GSL_FAILURE' and leaves C + unmodified. Starting with the first combination and repeatedly + applying this function will iterate through all possible + combinations of a given order. + + -- Function: int gsl_combination_prev (gsl_combination * C) + This function steps backwards from the combination C to the + previous combination in lexicographic order, returning + `GSL_SUCCESS'. If no previous combination is available it returns + `GSL_FAILURE' and leaves C unmodified. + + +File: gsl-ref.info, Node: Reading and writing combinations, Next: Combination Examples, Prev: Combination functions, Up: Combinations + +10.6 Reading and writing combinations +===================================== + +The library provides functions for reading and writing combinations to a +file as binary data or formatted text. + + -- Function: int gsl_combination_fwrite (FILE * STREAM, const + gsl_combination * C) + This function writes the elements of the combination C to the + stream STREAM in binary format. The function returns + `GSL_EFAILED' if there was a problem writing to the file. Since + the data is written in the native binary format it may not be + portable between different architectures. + + -- Function: int gsl_combination_fread (FILE * STREAM, gsl_combination + * C) + This function reads elements from the open stream STREAM into the + combination C in binary format. The combination C must be + preallocated with correct values of n and k since the function + uses the size of C to determine how many bytes to read. The + function returns `GSL_EFAILED' if there was a problem reading from + the file. The data is assumed to have been written in the native + binary format on the same architecture. + + -- Function: int gsl_combination_fprintf (FILE * STREAM, const + gsl_combination * C, const char * FORMAT) + This function writes the elements of the combination C + line-by-line to the stream STREAM using the format specifier + FORMAT, which should be suitable for a type of SIZE_T. On a GNU + system the type modifier `Z' represents `size_t', so `"%Zu\n"' is + a suitable format. The function returns `GSL_EFAILED' if there + was a problem writing to the file. + + -- Function: int gsl_combination_fscanf (FILE * STREAM, + gsl_combination * C) + This function reads formatted data from the stream STREAM into the + combination C. The combination C must be preallocated with + correct values of n and k since the function uses the size of C to + determine how many numbers to read. The function returns + `GSL_EFAILED' if there was a problem reading from the file. + + +File: gsl-ref.info, Node: Combination Examples, Next: Combination References and Further Reading, Prev: Reading and writing combinations, Up: Combinations + +10.7 Examples +============= + +The example program below prints all subsets of the set {0,1,2,3} +ordered by size. Subsets of the same size are ordered +lexicographically. + + #include <stdio.h> + #include <gsl/gsl_combination.h> + + int + main (void) + { + gsl_combination * c; + size_t i; + + printf ("All subsets of {0,1,2,3} by size:\n") ; + for (i = 0; i <= 4; i++) + { + c = gsl_combination_calloc (4, i); + do + { + printf ("{"); + gsl_combination_fprintf (stdout, c, " %u"); + printf (" }\n"); + } + while (gsl_combination_next (c) == GSL_SUCCESS); + gsl_combination_free (c); + } + + return 0; + } + +Here is the output from the program, + + $ ./a.out + All subsets of {0,1,2,3} by size: + { } + { 0 } + { 1 } + { 2 } + { 3 } + { 0 1 } + { 0 2 } + { 0 3 } + { 1 2 } + { 1 3 } + { 2 3 } + { 0 1 2 } + { 0 1 3 } + { 0 2 3 } + { 1 2 3 } + { 0 1 2 3 } + +All 16 subsets are generated, and the subsets of each size are sorted +lexicographically. + + +File: gsl-ref.info, Node: Combination References and Further Reading, Prev: Combination Examples, Up: Combinations + +10.8 References and Further Reading +=================================== + +Further information on combinations can be found in, + + Donald L. Kreher, Douglas R. Stinson, `Combinatorial Algorithms: + Generation, Enumeration and Search', 1998, CRC Press LLC, ISBN + 084933988X + + + +File: gsl-ref.info, Node: Sorting, Next: BLAS Support, Prev: Combinations, Up: Top + +11 Sorting +********** + +This chapter describes functions for sorting data, both directly and +indirectly (using an index). All the functions use the "heapsort" +algorithm. Heapsort is an O(N \log N) algorithm which operates +in-place and does not require any additional storage. It also provides +consistent performance, the running time for its worst-case (ordered +data) being not significantly longer than the average and best cases. +Note that the heapsort algorithm does not preserve the relative ordering +of equal elements--it is an "unstable" sort. However the resulting +order of equal elements will be consistent across different platforms +when using these functions. + +* Menu: + +* Sorting objects:: +* Sorting vectors:: +* Selecting the k smallest or largest elements:: +* Computing the rank:: +* Sorting Examples:: +* Sorting References and Further Reading:: + + +File: gsl-ref.info, Node: Sorting objects, Next: Sorting vectors, Up: Sorting + +11.1 Sorting objects +==================== + +The following function provides a simple alternative to the standard +library function `qsort'. It is intended for systems lacking `qsort', +not as a replacement for it. The function `qsort' should be used +whenever possible, as it will be faster and can provide stable ordering +of equal elements. Documentation for `qsort' is available in the `GNU +C Library Reference Manual'. + + The functions described in this section are defined in the header +file `gsl_heapsort.h'. + + -- Function: void gsl_heapsort (void * ARRAY, size_t COUNT, size_t + SIZE, gsl_comparison_fn_t COMPARE) + This function sorts the COUNT elements of the array ARRAY, each of + size SIZE, into ascending order using the comparison function + COMPARE. The type of the comparison function is defined by, + + int (*gsl_comparison_fn_t) (const void * a, + const void * b) + + A comparison function should return a negative integer if the first + argument is less than the second argument, `0' if the two arguments + are equal and a positive integer if the first argument is greater + than the second argument. + + For example, the following function can be used to sort doubles + into ascending numerical order. + + int + compare_doubles (const double * a, + const double * b) + { + if (*a > *b) + return 1; + else if (*a < *b) + return -1; + else + return 0; + } + + The appropriate function call to perform the sort is, + + gsl_heapsort (array, count, sizeof(double), + compare_doubles); + + Note that unlike `qsort' the heapsort algorithm cannot be made into + a stable sort by pointer arithmetic. The trick of comparing + pointers for equal elements in the comparison function does not + work for the heapsort algorithm. The heapsort algorithm performs + an internal rearrangement of the data which destroys its initial + ordering. + + -- Function: int gsl_heapsort_index (size_t * P, const void * ARRAY, + size_t COUNT, size_t SIZE, gsl_comparison_fn_t COMPARE) + This function indirectly sorts the COUNT elements of the array + ARRAY, each of size SIZE, into ascending order using the + comparison function COMPARE. The resulting permutation is stored + in P, an array of length N. The elements of P give the index of + the array element which would have been stored in that position if + the array had been sorted in place. The first element of P gives + the index of the least element in ARRAY, and the last element of P + gives the index of the greatest element in ARRAY. The array + itself is not changed. + + +File: gsl-ref.info, Node: Sorting vectors, Next: Selecting the k smallest or largest elements, Prev: Sorting objects, Up: Sorting + +11.2 Sorting vectors +==================== + +The following functions will sort the elements of an array or vector, +either directly or indirectly. They are defined for all real and +integer types using the normal suffix rules. For example, the `float' +versions of the array functions are `gsl_sort_float' and +`gsl_sort_float_index'. The corresponding vector functions are +`gsl_sort_vector_float' and `gsl_sort_vector_float_index'. The +prototypes are available in the header files `gsl_sort_float.h' +`gsl_sort_vector_float.h'. The complete set of prototypes can be +included using the header files `gsl_sort.h' and `gsl_sort_vector.h'. + + There are no functions for sorting complex arrays or vectors, since +the ordering of complex numbers is not uniquely defined. To sort a +complex vector by magnitude compute a real vector containing the +magnitudes of the complex elements, and sort this vector indirectly. +The resulting index gives the appropriate ordering of the original +complex vector. + + -- Function: void gsl_sort (double * DATA, size_t STRIDE, size_t N) + This function sorts the N elements of the array DATA with stride + STRIDE into ascending numerical order. + + -- Function: void gsl_sort_vector (gsl_vector * V) + This function sorts the elements of the vector V into ascending + numerical order. + + -- Function: void gsl_sort_index (size_t * P, const double * DATA, + size_t STRIDE, size_t N) + This function indirectly sorts the N elements of the array DATA + with stride STRIDE into ascending order, storing the resulting + permutation in P. The array P must be allocated with a sufficient + length to store the N elements of the permutation. The elements + of P give the index of the array element which would have been + stored in that position if the array had been sorted in place. + The array DATA is not changed. + + -- Function: int gsl_sort_vector_index (gsl_permutation * P, const + gsl_vector * V) + This function indirectly sorts the elements of the vector V into + ascending order, storing the resulting permutation in P. The + elements of P give the index of the vector element which would + have been stored in that position if the vector had been sorted in + place. The first element of P gives the index of the least element + in V, and the last element of P gives the index of the greatest + element in V. The vector V is not changed. + |