1 This is gsl-ref.info, produced by makeinfo version 4.8 from
4 INFO-DIR-SECTION Scientific software
6 * gsl-ref: (gsl-ref). GNU Scientific Library - Reference
9 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
10 2005, 2006, 2007 The GSL Team.
12 Permission is granted to copy, distribute and/or modify this document
13 under the terms of the GNU Free Documentation License, Version 1.2 or
14 any later version published by the Free Software Foundation; with the
15 Invariant Sections being "GNU General Public License" and "Free Software
16 Needs Free Documentation", the Front-Cover text being "A GNU Manual",
17 and with the Back-Cover Text being (a) (see below). A copy of the
18 license is included in the section entitled "GNU Free Documentation
21 (a) The Back-Cover Text is: "You have the freedom to copy and modify
25 File: gsl-ref.info, Node: Top, Next: Introduction, Prev: (dir), Up: (dir)
30 This file documents the GNU Scientific Library (GSL), a collection of
31 numerical routines for scientific computing. It corresponds to release
32 1.11 of the library. Please report any errors in this manual to
35 More information about GSL can be found at the project homepage,
36 `http://www.gnu.org/software/gsl/'.
38 Printed copies of this manual can be purchased from Network Theory
39 Ltd at `http://www.network-theory.co.uk/gsl/manual/'. The money raised
40 from sales of the manual helps support the development of GSL.
42 A Japanese translation of this manual is available from the GSL
43 project homepage thanks to Daisuke Tominaga.
45 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
46 2005, 2006, 2007 The GSL Team.
48 Permission is granted to copy, distribute and/or modify this document
49 under the terms of the GNU Free Documentation License, Version 1.2 or
50 any later version published by the Free Software Foundation; with the
51 Invariant Sections being "GNU General Public License" and "Free Software
52 Needs Free Documentation", the Front-Cover text being "A GNU Manual",
53 and with the Back-Cover Text being (a) (see below). A copy of the
54 license is included in the section entitled "GNU Free Documentation
57 (a) The Back-Cover Text is: "You have the freedom to copy and modify
66 * Mathematical Functions::
70 * Vectors and Matrices::
77 * Fast Fourier Transforms::
78 * Numerical Integration::
79 * Random Number Generation::
80 * Quasi-Random Sequences::
81 * Random Number Distributions::
85 * Monte Carlo Integration::
86 * Simulated Annealing::
87 * Ordinary Differential Equations::
89 * Numerical Differentiation::
90 * Chebyshev Approximations::
91 * Series Acceleration::
92 * Wavelet Transforms::
93 * Discrete Hankel Transforms::
94 * One dimensional Root-Finding::
95 * One dimensional Minimization::
96 * Multidimensional Root-Finding::
97 * Multidimensional Minimization::
98 * Least-Squares Fitting::
99 * Nonlinear Least-Squares Fitting::
101 * Physical Constants::
102 * IEEE floating-point arithmetic::
103 * Debugging Numerical Programs::
104 * Contributors to GSL::
106 * GSL CBLAS Library::
107 * Free Software Needs Free Documentation::
108 * GNU General Public License::
109 * GNU Free Documentation License::
116 File: gsl-ref.info, Node: Introduction, Next: Using the library, Prev: Top, Up: Top
121 The GNU Scientific Library (GSL) is a collection of routines for
122 numerical computing. The routines have been written from scratch in C,
123 and present a modern Applications Programming Interface (API) for C
124 programmers, allowing wrappers to be written for very high level
125 languages. The source code is distributed under the GNU General Public
130 * Routines available in GSL::
131 * GSL is Free Software::
135 * Further Information::
136 * Conventions used in this manual::
139 File: gsl-ref.info, Node: Routines available in GSL, Next: GSL is Free Software, Up: Introduction
141 1.1 Routines available in GSL
142 =============================
144 The library covers a wide range of topics in numerical computing.
145 Routines are available for the following areas,
147 Complex Numbers Roots of Polynomials
148 Special Functions Vectors and Matrices
149 Permutations Combinations
151 Linear Algebra CBLAS Library
152 Fast Fourier Transforms Eigensystems
153 Random Numbers Quadrature
154 Random Distributions Quasi-Random Sequences
155 Histograms Statistics
156 Monte Carlo Integration N-Tuples
157 Differential Equations Simulated Annealing
158 Numerical Differentiation Interpolation
159 Series Acceleration Chebyshev Approximations
160 Root-Finding Discrete Hankel Transforms
161 Least-Squares Fitting Minimization
162 IEEE Floating-Point Physical Constants
163 Basis Splines Wavelets
165 The use of these routines is described in this manual. Each chapter
166 provides detailed definitions of the functions, followed by example
167 programs and references to the articles on which the algorithms are
170 Where possible the routines have been based on reliable public-domain
171 packages such as FFTPACK and QUADPACK, which the developers of GSL have
172 reimplemented in C with modern coding conventions.
175 File: gsl-ref.info, Node: GSL is Free Software, Next: Obtaining GSL, Prev: Routines available in GSL, Up: Introduction
177 1.2 GSL is Free Software
178 ========================
180 The subroutines in the GNU Scientific Library are "free software"; this
181 means that everyone is free to use them, and to redistribute them in
182 other free programs. The library is not in the public domain; it is
183 copyrighted and there are conditions on its distribution. These
184 conditions are designed to permit everything that a good cooperating
185 citizen would want to do. What is not allowed is to try to prevent
186 others from further sharing any version of the software that they might
189 Specifically, we want to make sure that you have the right to share
190 copies of programs that you are given which use the GNU Scientific
191 Library, that you receive their source code or else can get it if you
192 want it, that you can change these programs or use pieces of them in new
193 free programs, and that you know you can do these things.
195 To make sure that everyone has such rights, we have to forbid you to
196 deprive anyone else of these rights. For example, if you distribute
197 copies of any code which uses the GNU Scientific Library, you must give
198 the recipients all the rights that you have received. You must make
199 sure that they, too, receive or can get the source code, both to the
200 library and the code which uses it. And you must tell them their
201 rights. This means that the library should not be redistributed in
202 proprietary programs.
204 Also, for our own protection, we must make certain that everyone
205 finds out that there is no warranty for the GNU Scientific Library. If
206 these programs are modified by someone else and passed on, we want their
207 recipients to know that what they have is not what we distributed, so
208 that any problems introduced by others will not reflect on our
211 The precise conditions for the distribution of software related to
212 the GNU Scientific Library are found in the GNU General Public License
213 (*note GNU General Public License::). Further information about this
214 license is available from the GNU Project webpage `Frequently Asked
215 Questions about the GNU GPL',
217 `http://www.gnu.org/copyleft/gpl-faq.html'
219 The Free Software Foundation also operates a license consulting service
220 for commercial users (contact details available from
221 `http://www.fsf.org/').
224 File: gsl-ref.info, Node: Obtaining GSL, Next: No Warranty, Prev: GSL is Free Software, Up: Introduction
229 The source code for the library can be obtained in different ways, by
230 copying it from a friend, purchasing it on CDROM or downloading it from
231 the internet. A list of public ftp servers which carry the source code
232 can be found on the GNU website,
234 `http://www.gnu.org/software/gsl/'
236 The preferred platform for the library is a GNU system, which allows it
237 to take advantage of additional features in the GNU C compiler and GNU C
238 library. However, the library is fully portable and should compile on
239 most systems with a C compiler. Precompiled versions of the library
240 can be purchased from commercial redistributors listed on the website
243 Announcements of new releases, updates and other relevant events are
244 made on the `info-gsl@gnu.org' mailing list. To subscribe to this
245 low-volume list, send an email of the following form:
247 To: info-gsl-request@gnu.org
250 You will receive a response asking you to reply in order to confirm
254 File: gsl-ref.info, Node: No Warranty, Next: Reporting Bugs, Prev: Obtaining GSL, Up: Introduction
259 The software described in this manual has no warranty, it is provided
260 "as is". It is your responsibility to validate the behavior of the
261 routines and their accuracy using the source code provided, or to
262 purchase support and warranties from commercial redistributors. Consult
263 the GNU General Public license for further details (*note GNU General
267 File: gsl-ref.info, Node: Reporting Bugs, Next: Further Information, Prev: No Warranty, Up: Introduction
272 A list of known bugs can be found in the `BUGS' file included in the
273 GSL distribution. Details of compilation problems can be found in the
276 If you find a bug which is not listed in these files, please report
277 it to <bug-gsl@gnu.org>.
279 All bug reports should include:
281 * The version number of GSL
283 * The hardware and operating system
285 * The compiler used, including version number and compilation options
287 * A description of the bug behavior
289 * A short program which exercises the bug
291 It is useful if you can check whether the same problem occurs when the
292 library is compiled without optimization. Thank you.
294 Any errors or omissions in this manual can also be reported to the
298 File: gsl-ref.info, Node: Further Information, Next: Conventions used in this manual, Prev: Reporting Bugs, Up: Introduction
300 1.6 Further Information
301 =======================
303 Additional information, including online copies of this manual, links to
304 related projects, and mailing list archives are available from the
305 website mentioned above.
307 Any questions about the use and installation of the library can be
308 asked on the mailing list `help-gsl@gnu.org'. To subscribe to this
309 list, send an email of the following form:
311 To: help-gsl-request@gnu.org
314 This mailing list can be used to ask questions not covered by this
315 manual, and to contact the developers of the library.
317 If you would like to refer to the GNU Scientific Library in a journal
318 article, the recommended way is to cite this reference manual, e.g. `M.
319 Galassi et al, GNU Scientific Library Reference Manual (2nd Ed.), ISBN
322 If you want to give a url, use "`http://www.gnu.org/software/gsl/'".
325 File: gsl-ref.info, Node: Conventions used in this manual, Prev: Further Information, Up: Introduction
327 1.7 Conventions used in this manual
328 ===================================
330 This manual contains many examples which can be typed at the keyboard.
331 A command entered at the terminal is shown like this,
335 The first character on the line is the terminal prompt, and should not
336 be typed. The dollar sign `$' is used as the standard prompt in this
337 manual, although some systems may use a different character.
339 The examples assume the use of the GNU operating system. There may
340 be minor differences in the output on other systems. The commands for
341 setting environment variables use the Bourne shell syntax of the
342 standard GNU shell (`bash').
345 File: gsl-ref.info, Node: Using the library, Next: Error Handling, Prev: Introduction, Up: Top
350 This chapter describes how to compile programs that use GSL, and
351 introduces its conventions.
355 * An Example Program::
356 * Compiling and Linking::
358 * ANSI C Compliance::
361 * Portability functions::
362 * Alternative optimized functions::
363 * Support for different numeric types::
364 * Compatibility with C++::
365 * Aliasing of arrays::
367 * Deprecated Functions::
371 File: gsl-ref.info, Node: An Example Program, Next: Compiling and Linking, Up: Using the library
373 2.1 An Example Program
374 ======================
376 The following short program demonstrates the use of the library by
377 computing the value of the Bessel function J_0(x) for x=5,
380 #include <gsl/gsl_sf_bessel.h>
386 double y = gsl_sf_bessel_J0 (x);
387 printf ("J0(%g) = %.18e\n", x, y);
391 The output is shown below, and should be correct to double-precision
394 J0(5) = -1.775967713143382920e-01
396 The steps needed to compile this program are described in the following
399 ---------- Footnotes ----------
401 (1) The last few digits may vary slightly depending on the compiler
402 and platform used--this is normal.
405 File: gsl-ref.info, Node: Compiling and Linking, Next: Shared Libraries, Prev: An Example Program, Up: Using the library
407 2.2 Compiling and Linking
408 =========================
410 The library header files are installed in their own `gsl' directory.
411 You should write any preprocessor include statements with a `gsl/'
412 directory prefix thus,
414 #include <gsl/gsl_math.h>
416 If the directory is not installed on the standard search path of your
417 compiler you will also need to provide its location to the preprocessor
418 as a command line flag. The default location of the `gsl' directory is
419 `/usr/local/include/gsl'. A typical compilation command for a source
420 file `example.c' with the GNU C compiler `gcc' is,
422 $ gcc -Wall -I/usr/local/include -c example.c
424 This results in an object file `example.o'. The default include path
425 for `gcc' searches `/usr/local/include' automatically so the `-I'
426 option can actually be omitted when GSL is installed in its default
431 * Linking programs with the library::
432 * Linking with an alternative BLAS library::
435 File: gsl-ref.info, Node: Linking programs with the library, Next: Linking with an alternative BLAS library, Up: Compiling and Linking
437 2.2.1 Linking programs with the library
438 ---------------------------------------
440 The library is installed as a single file, `libgsl.a'. A shared
441 version of the library `libgsl.so' is also installed on systems that
442 support shared libraries. The default location of these files is
443 `/usr/local/lib'. If this directory is not on the standard search path
444 of your linker you will also need to provide its location as a command
447 To link against the library you need to specify both the main
448 library and a supporting CBLAS library, which provides standard basic
449 linear algebra subroutines. A suitable CBLAS implementation is
450 provided in the library `libgslcblas.a' if your system does not provide
451 one. The following example shows how to link an application with the
454 $ gcc -L/usr/local/lib example.o -lgsl -lgslcblas -lm
456 The default library path for `gcc' searches `/usr/local/lib'
457 automatically so the `-L' option can be omitted when GSL is installed
458 in its default location.
461 File: gsl-ref.info, Node: Linking with an alternative BLAS library, Prev: Linking programs with the library, Up: Compiling and Linking
463 2.2.2 Linking with an alternative BLAS library
464 ----------------------------------------------
466 The following command line shows how you would link the same application
467 with an alternative CBLAS library called `libcblas',
469 $ gcc example.o -lgsl -lcblas -lm
471 For the best performance an optimized platform-specific CBLAS library
472 should be used for `-lcblas'. The library must conform to the CBLAS
473 standard. The ATLAS package provides a portable high-performance BLAS
474 library with a CBLAS interface. It is free software and should be
475 installed for any work requiring fast vector and matrix operations.
476 The following command line will link with the ATLAS library and its
479 $ gcc example.o -lgsl -lcblas -latlas -lm
481 For more information see *Note BLAS Support::.
484 File: gsl-ref.info, Node: Shared Libraries, Next: ANSI C Compliance, Prev: Compiling and Linking, Up: Using the library
489 To run a program linked with the shared version of the library the
490 operating system must be able to locate the corresponding `.so' file at
491 runtime. If the library cannot be found, the following error will
495 ./a.out: error while loading shared libraries:
496 libgsl.so.0: cannot open shared object file: No such
499 To avoid this error, define the shell variable `LD_LIBRARY_PATH' to
500 include the directory where the library is installed.
502 For example, in the Bourne shell (`/bin/sh' or `/bin/bash'), the
503 library search path can be set with the following commands:
505 $ LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
506 $ export LD_LIBRARY_PATH
509 In the C-shell (`/bin/csh' or `/bin/tcsh') the equivalent command is,
511 % setenv LD_LIBRARY_PATH /usr/local/lib:$LD_LIBRARY_PATH
513 The standard prompt for the C-shell in the example above is the percent
514 character `%', and should not be typed as part of the command.
516 To save retyping these commands each session they should be placed
517 in an individual or system-wide login file.
519 To compile a statically linked version of the program, use the
520 `-static' flag in `gcc',
522 $ gcc -static example.o -lgsl -lgslcblas -lm
525 File: gsl-ref.info, Node: ANSI C Compliance, Next: Inline functions, Prev: Shared Libraries, Up: Using the library
527 2.4 ANSI C Compliance
528 =====================
530 The library is written in ANSI C and is intended to conform to the ANSI
531 C standard (C89). It should be portable to any system with a working
534 The library does not rely on any non-ANSI extensions in the
535 interface it exports to the user. Programs you write using GSL can be
536 ANSI compliant. Extensions which can be used in a way compatible with
537 pure ANSI C are supported, however, via conditional compilation. This
538 allows the library to take advantage of compiler extensions on those
539 platforms which support them.
541 When an ANSI C feature is known to be broken on a particular system
542 the library will exclude any related functions at compile-time. This
543 should make it impossible to link a program that would use these
544 functions and give incorrect results.
546 To avoid namespace conflicts all exported function names and
547 variables have the prefix `gsl_', while exported macros have the prefix
551 File: gsl-ref.info, Node: Inline functions, Next: Long double, Prev: ANSI C Compliance, Up: Using the library
556 The `inline' keyword is not part of the original ANSI C standard (C89)
557 and the library does not export any inline function definitions by
558 default. However, the library provides optional inline versions of
559 performance-critical functions by conditional compilation. The inline
560 versions of these functions can be included by defining the macro
561 `HAVE_INLINE' when compiling an application,
563 $ gcc -Wall -c -DHAVE_INLINE example.c
565 If you use `autoconf' this macro can be defined automatically. If you
566 do not define the macro `HAVE_INLINE' then the slower non-inlined
567 versions of the functions will be used instead.
569 Note that the actual usage of the inline keyword is `extern inline',
570 which eliminates unnecessary function definitions in GCC. If the form
571 `extern inline' causes problems with other compilers a stricter
572 autoconf test can be used, see *Note Autoconf Macros::.
575 File: gsl-ref.info, Node: Long double, Next: Portability functions, Prev: Inline functions, Up: Using the library
580 In general, the algorithms in the library are written for double
581 precision only. The `long double' type is not supported for actual
584 One reason for this choice is that the precision of `long double' is
585 platform dependent. The IEEE standard only specifies the minimum
586 precision of extended precision numbers, while the precision of
587 `double' is the same on all platforms.
589 However, it is sometimes necessary to interact with external data in
590 long-double format, so the vector and matrix datatypes include
591 long-double versions.
593 It should be noted that in some system libraries the `stdio.h'
594 formatted input/output functions `printf' and `scanf' are not
595 implemented correctly for `long double'. Undefined or incorrect
596 results are avoided by testing these functions during the `configure'
597 stage of library compilation and eliminating certain GSL functions
598 which depend on them if necessary. The corresponding line in the
599 `configure' output looks like this,
601 checking whether printf works with long double... no
603 Consequently when `long double' formatted input/output does not work on
604 a given system it should be impossible to link a program which uses GSL
605 functions dependent on this.
607 If it is necessary to work on a system which does not support
608 formatted `long double' input/output then the options are to use binary
609 formats or to convert `long double' results into `double' for reading
613 File: gsl-ref.info, Node: Portability functions, Next: Alternative optimized functions, Prev: Long double, Up: Using the library
615 2.7 Portability functions
616 =========================
618 To help in writing portable applications GSL provides some
619 implementations of functions that are found in other libraries, such as
620 the BSD math library. You can write your application to use the native
621 versions of these functions, and substitute the GSL versions via a
622 preprocessor macro if they are unavailable on another platform.
624 For example, after determining whether the BSD function `hypot' is
625 available you can include the following macro definitions in a file
626 `config.h' with your application,
628 /* Substitute gsl_hypot for missing system hypot */
631 #define hypot gsl_hypot
634 The application source files can then use the include command `#include
635 <config.h>' to replace each occurrence of `hypot' by `gsl_hypot' when
636 `hypot' is not available. This substitution can be made automatically
637 if you use `autoconf', see *Note Autoconf Macros::.
639 In most circumstances the best strategy is to use the native
640 versions of these functions when available, and fall back to GSL
641 versions otherwise, since this allows your application to take
642 advantage of any platform-specific optimizations in the system library.
643 This is the strategy used within GSL itself.
646 File: gsl-ref.info, Node: Alternative optimized functions, Next: Support for different numeric types, Prev: Portability functions, Up: Using the library
648 2.8 Alternative optimized functions
649 ===================================
651 The main implementation of some functions in the library will not be
652 optimal on all architectures. For example, there are several ways to
653 compute a Gaussian random variate and their relative speeds are
654 platform-dependent. In cases like this the library provides alternative
655 implementations of these functions with the same interface. If you
656 write your application using calls to the standard implementation you
657 can select an alternative version later via a preprocessor definition.
658 It is also possible to introduce your own optimized functions this way
659 while retaining portability. The following lines demonstrate the use of
660 a platform-dependent choice of methods for sampling from the Gaussian
664 #define gsl_ran_gaussian gsl_ran_gaussian_ratio_method
667 #define gsl_ran_gaussian my_gaussian
670 These lines would be placed in the configuration header file `config.h'
671 of the application, which should then be included by all the source
672 files. Note that the alternative implementations will not produce
673 bit-for-bit identical results, and in the case of random number
674 distributions will produce an entirely different stream of random
678 File: gsl-ref.info, Node: Support for different numeric types, Next: Compatibility with C++, Prev: Alternative optimized functions, Up: Using the library
680 2.9 Support for different numeric types
681 =======================================
683 Many functions in the library are defined for different numeric types.
684 This feature is implemented by varying the name of the function with a
685 type-related modifier--a primitive form of C++ templates. The modifier
686 is inserted into the function name after the initial module prefix.
687 The following table shows the function names defined for all the
688 numeric types of an imaginary module `gsl_foo' with function `fn',
691 gsl_foo_long_double_fn long double
692 gsl_foo_float_fn float
694 gsl_foo_ulong_fn unsigned long
696 gsl_foo_uint_fn unsigned int
697 gsl_foo_short_fn short
698 gsl_foo_ushort_fn unsigned short
700 gsl_foo_uchar_fn unsigned char
702 The normal numeric precision `double' is considered the default and
703 does not require a suffix. For example, the function `gsl_stats_mean'
704 computes the mean of double precision numbers, while the function
705 `gsl_stats_int_mean' computes the mean of integers.
707 A corresponding scheme is used for library defined types, such as
708 `gsl_vector' and `gsl_matrix'. In this case the modifier is appended
709 to the type name. For example, if a module defines a new
710 type-dependent struct or typedef `gsl_foo' it is modified for other
711 types in the following way,
714 gsl_foo_long_double long double
717 gsl_foo_ulong unsigned long
719 gsl_foo_uint unsigned int
721 gsl_foo_ushort unsigned short
723 gsl_foo_uchar unsigned char
725 When a module contains type-dependent definitions the library provides
726 individual header files for each type. The filenames are modified as
727 shown in the below. For convenience the default header includes the
728 definitions for all the types. To include only the double precision
729 header file, or any other specific type, use its individual filename.
731 #include <gsl/gsl_foo.h> All types
732 #include <gsl/gsl_foo_double.h> double
733 #include <gsl/gsl_foo_long_double.h> long double
734 #include <gsl/gsl_foo_float.h> float
735 #include <gsl/gsl_foo_long.h> long
736 #include <gsl/gsl_foo_ulong.h> unsigned long
737 #include <gsl/gsl_foo_int.h> int
738 #include <gsl/gsl_foo_uint.h> unsigned int
739 #include <gsl/gsl_foo_short.h> short
740 #include <gsl/gsl_foo_ushort.h> unsigned short
741 #include <gsl/gsl_foo_char.h> char
742 #include <gsl/gsl_foo_uchar.h> unsigned char
745 File: gsl-ref.info, Node: Compatibility with C++, Next: Aliasing of arrays, Prev: Support for different numeric types, Up: Using the library
747 2.10 Compatibility with C++
748 ===========================
750 The library header files automatically define functions to have `extern
751 "C"' linkage when included in C++ programs. This allows the functions
752 to be called directly from C++.
754 To use C++ exception handling within user-defined functions passed to
755 the library as parameters, the library must be built with the
756 additional `CFLAGS' compilation option `-fexceptions'.
759 File: gsl-ref.info, Node: Aliasing of arrays, Next: Thread-safety, Prev: Compatibility with C++, Up: Using the library
761 2.11 Aliasing of arrays
762 =======================
764 The library assumes that arrays, vectors and matrices passed as
765 modifiable arguments are not aliased and do not overlap with each other.
766 This removes the need for the library to handle overlapping memory
767 regions as a special case, and allows additional optimizations to be
768 used. If overlapping memory regions are passed as modifiable arguments
769 then the results of such functions will be undefined. If the arguments
770 will not be modified (for example, if a function prototype declares them
771 as `const' arguments) then overlapping or aliased memory regions can be
775 File: gsl-ref.info, Node: Thread-safety, Next: Deprecated Functions, Prev: Aliasing of arrays, Up: Using the library
780 The library can be used in multi-threaded programs. All the functions
781 are thread-safe, in the sense that they do not use static variables.
782 Memory is always associated with objects and not with functions. For
783 functions which use "workspace" objects as temporary storage the
784 workspaces should be allocated on a per-thread basis. For functions
785 which use "table" objects as read-only memory the tables can be used by
786 multiple threads simultaneously. Table arguments are always declared
787 `const' in function prototypes, to indicate that they may be safely
788 accessed by different threads.
790 There are a small number of static global variables which are used to
791 control the overall behavior of the library (e.g. whether to use
792 range-checking, the function to call on fatal error, etc). These
793 variables are set directly by the user, so they should be initialized
794 once at program startup and not modified by different threads.
797 File: gsl-ref.info, Node: Deprecated Functions, Next: Code Reuse, Prev: Thread-safety, Up: Using the library
799 2.13 Deprecated Functions
800 =========================
802 From time to time, it may be necessary for the definitions of some
803 functions to be altered or removed from the library. In these
804 circumstances the functions will first be declared "deprecated" and
805 then removed from subsequent versions of the library. Functions that
806 are deprecated can be disabled in the current release by setting the
807 preprocessor definition `GSL_DISABLE_DEPRECATED'. This allows existing
808 code to be tested for forwards compatibility.
811 File: gsl-ref.info, Node: Code Reuse, Prev: Deprecated Functions, Up: Using the library
816 Where possible the routines in the library have been written to avoid
817 dependencies between modules and files. This should make it possible to
818 extract individual functions for use in your own applications, without
819 needing to have the whole library installed. You may need to define
820 certain macros such as `GSL_ERROR' and remove some `#include'
821 statements in order to compile the files as standalone units. Reuse of
822 the library code in this way is encouraged, subject to the terms of the
823 GNU General Public License.
826 File: gsl-ref.info, Node: Error Handling, Next: Mathematical Functions, Prev: Using the library, Up: Top
831 This chapter describes the way that GSL functions report and handle
832 errors. By examining the status information returned by every function
833 you can determine whether it succeeded or failed, and if it failed you
834 can find out what the precise cause of failure was. You can also define
835 your own error handling functions to modify the default behavior of the
838 The functions described in this section are declared in the header
846 * Using GSL error reporting in your own functions::
847 * Error Reporting Examples::
850 File: gsl-ref.info, Node: Error Reporting, Next: Error Codes, Up: Error Handling
855 The library follows the thread-safe error reporting conventions of the
856 POSIX Threads library. Functions return a non-zero error code to
857 indicate an error and `0' to indicate success.
859 int status = gsl_function (...)
861 if (status) { /* an error occurred */
863 /* status value specifies the type of error */
866 The routines report an error whenever they cannot perform the task
867 requested of them. For example, a root-finding function would return a
868 non-zero error code if could not converge to the requested accuracy, or
869 exceeded a limit on the number of iterations. Situations like this are
870 a normal occurrence when using any mathematical library and you should
871 check the return status of the functions that you call.
873 Whenever a routine reports an error the return value specifies the
874 type of error. The return value is analogous to the value of the
875 variable `errno' in the C library. The caller can examine the return
876 code and decide what action to take, including ignoring the error if it
877 is not considered serious.
879 In addition to reporting errors by return codes the library also has
880 an error handler function `gsl_error'. This function is called by
881 other library functions when they report an error, just before they
882 return to the caller. The default behavior of the error handler is to
883 print a message and abort the program,
885 gsl: file.c:67: ERROR: invalid argument supplied by user
886 Default GSL error handler invoked.
889 The purpose of the `gsl_error' handler is to provide a function
890 where a breakpoint can be set that will catch library errors when
891 running under the debugger. It is not intended for use in production
892 programs, which should handle any errors using the return codes.
895 File: gsl-ref.info, Node: Error Codes, Next: Error Handlers, Prev: Error Reporting, Up: Error Handling
900 The error code numbers returned by library functions are defined in the
901 file `gsl_errno.h'. They all have the prefix `GSL_' and expand to
902 non-zero constant integer values. Error codes above 1024 are reserved
903 for applications, and are not used by the library. Many of the error
904 codes use the same base name as the corresponding error code in the C
905 library. Here are some of the most common error codes,
907 -- Macro: int GSL_EDOM
908 Domain error; used by mathematical functions when an argument
909 value does not fall into the domain over which the function is
910 defined (like EDOM in the C library)
912 -- Macro: int GSL_ERANGE
913 Range error; used by mathematical functions when the result value
914 is not representable because of overflow or underflow (like ERANGE
917 -- Macro: int GSL_ENOMEM
918 No memory available. The system cannot allocate more virtual
919 memory because its capacity is full (like ENOMEM in the C
920 library). This error is reported when a GSL routine encounters
921 problems when trying to allocate memory with `malloc'.
923 -- Macro: int GSL_EINVAL
924 Invalid argument. This is used to indicate various kinds of
925 problems with passing the wrong argument to a library function
926 (like EINVAL in the C library).
928 The error codes can be converted into an error message using the
929 function `gsl_strerror'.
931 -- Function: const char * gsl_strerror (const int GSL_ERRNO)
932 This function returns a pointer to a string describing the error
933 code GSL_ERRNO. For example,
935 printf ("error: %s\n", gsl_strerror (status));
937 would print an error message like `error: output range error' for a
938 status value of `GSL_ERANGE'.
941 File: gsl-ref.info, Node: Error Handlers, Next: Using GSL error reporting in your own functions, Prev: Error Codes, Up: Error Handling
946 The default behavior of the GSL error handler is to print a short
947 message and call `abort'. When this default is in use programs will
948 stop with a core-dump whenever a library routine reports an error.
949 This is intended as a fail-safe default for programs which do not check
950 the return status of library routines (we don't encourage you to write
953 If you turn off the default error handler it is your responsibility
954 to check the return values of routines and handle them yourself. You
955 can also customize the error behavior by providing a new error handler.
956 For example, an alternative error handler could log all errors to a
957 file, ignore certain error conditions (such as underflows), or start the
958 debugger and attach it to the current process when an error occurs.
960 All GSL error handlers have the type `gsl_error_handler_t', which is
961 defined in `gsl_errno.h',
963 -- Data Type: gsl_error_handler_t
964 This is the type of GSL error handler functions. An error handler
965 will be passed four arguments which specify the reason for the
966 error (a string), the name of the source file in which it occurred
967 (also a string), the line number in that file (an integer) and the
968 error number (an integer). The source file and line number are
969 set at compile time using the `__FILE__' and `__LINE__' directives
970 in the preprocessor. An error handler function returns type
971 `void'. Error handler functions should be defined like this,
973 void handler (const char * reason,
978 To request the use of your own error handler you need to call the
979 function `gsl_set_error_handler' which is also declared in
982 -- Function: gsl_error_handler_t * gsl_set_error_handler
983 (gsl_error_handler_t * NEW_HANDLER)
984 This function sets a new error handler, NEW_HANDLER, for the GSL
985 library routines. The previous handler is returned (so that you
986 can restore it later). Note that the pointer to a user defined
987 error handler function is stored in a static variable, so there
988 can be only one error handler per program. This function should
989 be not be used in multi-threaded programs except to set up a
990 program-wide error handler from a master thread. The following
991 example shows how to set and restore a new error handler,
993 /* save original handler, install new handler */
994 old_handler = gsl_set_error_handler (&my_handler);
996 /* code uses new handler */
999 /* restore original handler */
1000 gsl_set_error_handler (old_handler);
1002 To use the default behavior (`abort' on error) set the error
1005 old_handler = gsl_set_error_handler (NULL);
1007 -- Function: gsl_error_handler_t * gsl_set_error_handler_off ()
1008 This function turns off the error handler by defining an error
1009 handler which does nothing. This will cause the program to
1010 continue after any error, so the return values from any library
1011 routines must be checked. This is the recommended behavior for
1012 production programs. The previous handler is returned (so that
1013 you can restore it later).
1015 The error behavior can be changed for specific applications by
1016 recompiling the library with a customized definition of the `GSL_ERROR'
1017 macro in the file `gsl_errno.h'.
1020 File: gsl-ref.info, Node: Using GSL error reporting in your own functions, Next: Error Reporting Examples, Prev: Error Handlers, Up: Error Handling
1022 3.4 Using GSL error reporting in your own functions
1023 ===================================================
1025 If you are writing numerical functions in a program which also uses GSL
1026 code you may find it convenient to adopt the same error reporting
1027 conventions as in the library.
1029 To report an error you need to call the function `gsl_error' with a
1030 string describing the error and then return an appropriate error code
1031 from `gsl_errno.h', or a special value, such as `NaN'. For convenience
1032 the file `gsl_errno.h' defines two macros which carry out these steps:
1034 -- Macro: GSL_ERROR (REASON, GSL_ERRNO)
1035 This macro reports an error using the GSL conventions and returns a
1036 status value of `gsl_errno'. It expands to the following code
1039 gsl_error (reason, __FILE__, __LINE__, gsl_errno);
1042 The macro definition in `gsl_errno.h' actually wraps the code in a
1043 `do { ... } while (0)' block to prevent possible parsing problems.
1045 Here is an example of how the macro could be used to report that a
1046 routine did not achieve a requested tolerance. To report the error the
1047 routine needs to return the error code `GSL_ETOL'.
1049 if (residual > tolerance)
1051 GSL_ERROR("residual exceeds tolerance", GSL_ETOL);
1054 -- Macro: GSL_ERROR_VAL (REASON, GSL_ERRNO, VALUE)
1055 This macro is the same as `GSL_ERROR' but returns a user-defined
1056 value of VALUE instead of an error code. It can be used for
1057 mathematical functions that return a floating point value.
1059 The following example shows how to return a `NaN' at a mathematical
1060 singularity using the `GSL_ERROR_VAL' macro,
1064 GSL_ERROR_VAL("argument lies on singularity",
1065 GSL_ERANGE, GSL_NAN);
1069 File: gsl-ref.info, Node: Error Reporting Examples, Prev: Using GSL error reporting in your own functions, Up: Error Handling
1074 Here is an example of some code which checks the return value of a
1075 function where an error might be reported,
1078 #include <gsl/gsl_errno.h>
1079 #include <gsl/gsl_fft_complex.h>
1085 gsl_set_error_handler_off();
1087 status = gsl_fft_complex_radix2_forward (data, stride, n);
1090 if (status == GSL_EINVAL) {
1091 fprintf (stderr, "invalid argument, n=%d\n", n);
1093 fprintf (stderr, "failed, gsl_errno=%d\n",
1100 The function `gsl_fft_complex_radix2' only accepts integer lengths
1101 which are a power of two. If the variable `n' is not a power of two
1102 then the call to the library function will return `GSL_EINVAL',
1103 indicating that the length argument is invalid. The function call to
1104 `gsl_set_error_handler_off' stops the default error handler from
1105 aborting the program. The `else' clause catches any other possible
1109 File: gsl-ref.info, Node: Mathematical Functions, Next: Complex Numbers, Prev: Error Handling, Up: Top
1111 4 Mathematical Functions
1112 ************************
1114 This chapter describes basic mathematical functions. Some of these
1115 functions are present in system libraries, but the alternative versions
1116 given here can be used as a substitute when the system functions are not
1119 The functions and macros described in this chapter are defined in the
1120 header file `gsl_math.h'.
1124 * Mathematical Constants::
1125 * Infinities and Not-a-number::
1126 * Elementary Functions::
1127 * Small integer powers::
1128 * Testing the Sign of Numbers::
1129 * Testing for Odd and Even Numbers::
1130 * Maximum and Minimum functions::
1131 * Approximate Comparison of Floating Point Numbers::
1134 File: gsl-ref.info, Node: Mathematical Constants, Next: Infinities and Not-a-number, Up: Mathematical Functions
1136 4.1 Mathematical Constants
1137 ==========================
1139 The library ensures that the standard BSD mathematical constants are
1140 defined. For reference, here is a list of the constants:
1143 The base of exponentials, e
1146 The base-2 logarithm of e, \log_2 (e)
1149 The base-10 logarithm of e, \log_10 (e)
1152 The square root of two, \sqrt 2
1155 The square root of one-half, \sqrt{1/2}
1158 The square root of three, \sqrt 3
1161 The constant pi, \pi
1164 Pi divided by two, \pi/2
1167 Pi divided by four, \pi/4
1170 The square root of pi, \sqrt\pi
1173 Two divided by the square root of pi, 2/\sqrt\pi
1176 The reciprocal of pi, 1/\pi
1179 Twice the reciprocal of pi, 2/\pi
1182 The natural logarithm of ten, \ln(10)
1185 The natural logarithm of two, \ln(2)
1188 The natural logarithm of pi, \ln(\pi)
1191 Euler's constant, \gamma
1195 File: gsl-ref.info, Node: Infinities and Not-a-number, Next: Elementary Functions, Prev: Mathematical Constants, Up: Mathematical Functions
1197 4.2 Infinities and Not-a-number
1198 ===============================
1200 -- Macro: GSL_POSINF
1201 This macro contains the IEEE representation of positive infinity,
1202 +\infty. It is computed from the expression `+1.0/0.0'.
1204 -- Macro: GSL_NEGINF
1205 This macro contains the IEEE representation of negative infinity,
1206 -\infty. It is computed from the expression `-1.0/0.0'.
1209 This macro contains the IEEE representation of the Not-a-Number
1210 symbol, `NaN'. It is computed from the ratio `0.0/0.0'.
1212 -- Function: int gsl_isnan (const double X)
1213 This function returns 1 if X is not-a-number.
1215 -- Function: int gsl_isinf (const double X)
1216 This function returns +1 if X is positive infinity, -1 if X is
1217 negative infinity and 0 otherwise.
1219 -- Function: int gsl_finite (const double X)
1220 This function returns 1 if X is a real number, and 0 if it is
1221 infinite or not-a-number.
1224 File: gsl-ref.info, Node: Elementary Functions, Next: Small integer powers, Prev: Infinities and Not-a-number, Up: Mathematical Functions
1226 4.3 Elementary Functions
1227 ========================
1229 The following routines provide portable implementations of functions
1230 found in the BSD math library. When native versions are not available
1231 the functions described here can be used instead. The substitution can
1232 be made automatically if you use `autoconf' to compile your application
1233 (*note Portability functions::).
1235 -- Function: double gsl_log1p (const double X)
1236 This function computes the value of \log(1+x) in a way that is
1237 accurate for small X. It provides an alternative to the BSD math
1238 function `log1p(x)'.
1240 -- Function: double gsl_expm1 (const double X)
1241 This function computes the value of \exp(x)-1 in a way that is
1242 accurate for small X. It provides an alternative to the BSD math
1243 function `expm1(x)'.
1245 -- Function: double gsl_hypot (const double X, const double Y)
1246 This function computes the value of \sqrt{x^2 + y^2} in a way that
1247 avoids overflow. It provides an alternative to the BSD math
1248 function `hypot(x,y)'.
1250 -- Function: double gsl_hypot3 (const double X, const double Y, const
1252 This function computes the value of \sqrt{x^2 + y^2 + z^2} in a
1253 way that avoids overflow.
1255 -- Function: double gsl_acosh (const double X)
1256 This function computes the value of \arccosh(x). It provides an
1257 alternative to the standard math function `acosh(x)'.
1259 -- Function: double gsl_asinh (const double X)
1260 This function computes the value of \arcsinh(x). It provides an
1261 alternative to the standard math function `asinh(x)'.
1263 -- Function: double gsl_atanh (const double X)
1264 This function computes the value of \arctanh(x). It provides an
1265 alternative to the standard math function `atanh(x)'.
1267 -- Function: double gsl_ldexp (double X, int E)
1268 This function computes the value of x * 2^e. It provides an
1269 alternative to the standard math function `ldexp(x,e)'.
1271 -- Function: double gsl_frexp (double X, int * E)
1272 This function splits the number x into its normalized fraction f
1273 and exponent e, such that x = f * 2^e and 0.5 <= f < 1. The
1274 function returns f and stores the exponent in e. If x is zero,
1275 both f and e are set to zero. This function provides an
1276 alternative to the standard math function `frexp(x, e)'.
1279 File: gsl-ref.info, Node: Small integer powers, Next: Testing the Sign of Numbers, Prev: Elementary Functions, Up: Mathematical Functions
1281 4.4 Small integer powers
1282 ========================
1284 A common complaint about the standard C library is its lack of a
1285 function for calculating (small) integer powers. GSL provides some
1286 simple functions to fill this gap. For reasons of efficiency, these
1287 functions do not check for overflow or underflow conditions.
1289 -- Function: double gsl_pow_int (double X, int N)
1290 This routine computes the power x^n for integer N. The power is
1291 computed efficiently--for example, x^8 is computed as ((x^2)^2)^2,
1292 requiring only 3 multiplications. A version of this function
1293 which also computes the numerical error in the result is available
1294 as `gsl_sf_pow_int_e'.
1296 -- Function: double gsl_pow_2 (const double X)
1297 -- Function: double gsl_pow_3 (const double X)
1298 -- Function: double gsl_pow_4 (const double X)
1299 -- Function: double gsl_pow_5 (const double X)
1300 -- Function: double gsl_pow_6 (const double X)
1301 -- Function: double gsl_pow_7 (const double X)
1302 -- Function: double gsl_pow_8 (const double X)
1303 -- Function: double gsl_pow_9 (const double X)
1304 These functions can be used to compute small integer powers x^2,
1305 x^3, etc. efficiently. The functions will be inlined when
1306 `HAVE_INLINE' is defined, so that use of these functions should be
1307 as efficient as explicitly writing the corresponding product
1310 #include <gsl/gsl_math.h>
1311 double y = gsl_pow_4 (3.141) /* compute 3.141**4 */
1314 File: gsl-ref.info, Node: Testing the Sign of Numbers, Next: Testing for Odd and Even Numbers, Prev: Small integer powers, Up: Mathematical Functions
1316 4.5 Testing the Sign of Numbers
1317 ===============================
1319 -- Macro: GSL_SIGN (x)
1320 This macro returns the sign of X. It is defined as `((x) >= 0 ? 1
1321 : -1)'. Note that with this definition the sign of zero is positive
1322 (regardless of its IEEE sign bit).
1325 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
1327 4.6 Testing for Odd and Even Numbers
1328 ====================================
1330 -- Macro: GSL_IS_ODD (n)
1331 This macro evaluates to 1 if N is odd and 0 if N is even. The
1332 argument N must be of integer type.
1334 -- Macro: GSL_IS_EVEN (n)
1335 This macro is the opposite of `GSL_IS_ODD(n)'. It evaluates to 1 if
1336 N is even and 0 if N is odd. The argument N must be of integer
1340 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
1342 4.7 Maximum and Minimum functions
1343 =================================
1345 -- Macro: GSL_MAX (a, b)
1346 This macro returns the maximum of A and B. It is defined as `((a)
1349 -- Macro: GSL_MIN (a, b)
1350 This macro returns the minimum of A and B. It is defined as `((a)
1353 -- Function: extern inline double GSL_MAX_DBL (double A, double B)
1354 This function returns the maximum of the double precision numbers
1355 A and B using an inline function. The use of a function allows for
1356 type checking of the arguments as an extra safety feature. On
1357 platforms where inline functions are not available the macro
1358 `GSL_MAX' will be automatically substituted.
1360 -- Function: extern inline double GSL_MIN_DBL (double A, double B)
1361 This function returns the minimum of the double precision numbers
1362 A and B using an inline function. The use of a function allows for
1363 type checking of the arguments as an extra safety feature. On
1364 platforms where inline functions are not available the macro
1365 `GSL_MIN' will be automatically substituted.
1367 -- Function: extern inline int GSL_MAX_INT (int A, int B)
1368 -- Function: extern inline int GSL_MIN_INT (int A, int B)
1369 These functions return the maximum or minimum of the integers A
1370 and B using an inline function. On platforms where inline
1371 functions are not available the macros `GSL_MAX' or `GSL_MIN' will
1372 be automatically substituted.
1374 -- Function: extern inline long double GSL_MAX_LDBL (long double A,
1376 -- Function: extern inline long double GSL_MIN_LDBL (long double A,
1378 These functions return the maximum or minimum of the long doubles A
1379 and B using an inline function. On platforms where inline
1380 functions are not available the macros `GSL_MAX' or `GSL_MIN' will
1381 be automatically substituted.
1384 File: gsl-ref.info, Node: Approximate Comparison of Floating Point Numbers, Prev: Maximum and Minimum functions, Up: Mathematical Functions
1386 4.8 Approximate Comparison of Floating Point Numbers
1387 ====================================================
1389 It is sometimes useful to be able to compare two floating point numbers
1390 approximately, to allow for rounding and truncation errors. The
1391 following function implements the approximate floating-point comparison
1392 algorithm proposed by D.E. Knuth in Section 4.2.2 of `Seminumerical
1393 Algorithms' (3rd edition).
1395 -- Function: int gsl_fcmp (double X, double Y, double EPSILON)
1396 This function determines whether x and y are approximately equal
1397 to a relative accuracy EPSILON.
1399 The relative accuracy is measured using an interval of size 2
1400 \delta, where \delta = 2^k \epsilon and k is the maximum base-2
1401 exponent of x and y as computed by the function `frexp'.
1403 If x and y lie within this interval, they are considered
1404 approximately equal and the function returns 0. Otherwise if x <
1405 y, the function returns -1, or if x > y, the function returns +1.
1407 Note that x and y are compared to relative accuracy, so this
1408 function is not suitable for testing whether a value is
1411 The implementation is based on the package `fcmp' by T.C. Belding.
1414 File: gsl-ref.info, Node: Complex Numbers, Next: Polynomials, Prev: Mathematical Functions, Up: Top
1419 The functions described in this chapter provide support for complex
1420 numbers. The algorithms take care to avoid unnecessary intermediate
1421 underflows and overflows, allowing the functions to be evaluated over
1422 as much of the complex plane as possible.
1424 For multiple-valued functions the branch cuts have been chosen to
1425 follow the conventions of Abramowitz and Stegun in the `Handbook of
1426 Mathematical Functions'. The functions return principal values which are
1427 the same as those in GNU Calc, which in turn are the same as those in
1428 `Common Lisp, The Language (Second Edition)'(1) and the HP-28/48 series
1431 The complex types are defined in the header file `gsl_complex.h',
1432 while the corresponding complex functions and arithmetic operations are
1433 defined in `gsl_complex_math.h'.
1438 * Properties of complex numbers::
1439 * Complex arithmetic operators::
1440 * Elementary Complex Functions::
1441 * Complex Trigonometric Functions::
1442 * Inverse Complex Trigonometric Functions::
1443 * Complex Hyperbolic Functions::
1444 * Inverse Complex Hyperbolic Functions::
1445 * Complex Number References and Further Reading::
1447 ---------- Footnotes ----------
1449 (1) Note that the first edition uses different definitions.
1452 File: gsl-ref.info, Node: Complex numbers, Next: Properties of complex numbers, Up: Complex Numbers
1457 Complex numbers are represented using the type `gsl_complex'. The
1458 internal representation of this type may vary across platforms and
1459 should not be accessed directly. The functions and macros described
1460 below allow complex numbers to be manipulated in a portable way.
1462 For reference, the default form of the `gsl_complex' type is given
1463 by the following struct,
1470 The real and imaginary part are stored in contiguous elements of a two
1471 element array. This eliminates any padding between the real and
1472 imaginary parts, `dat[0]' and `dat[1]', allowing the struct to be
1473 mapped correctly onto packed complex arrays.
1475 -- Function: gsl_complex gsl_complex_rect (double X, double Y)
1476 This function uses the rectangular cartesian components (X,Y) to
1477 return the complex number z = x + i y. An inline version of this
1478 function is used when `HAVE_INLINE' is defined.
1480 -- Function: gsl_complex gsl_complex_polar (double R, double THETA)
1481 This function returns the complex number z = r \exp(i \theta) = r
1482 (\cos(\theta) + i \sin(\theta)) from the polar representation
1485 -- Macro: GSL_REAL (Z)
1486 -- Macro: GSL_IMAG (Z)
1487 These macros return the real and imaginary parts of the complex
1490 -- Macro: GSL_SET_COMPLEX (ZP, X, Y)
1491 This macro uses the cartesian components (X,Y) to set the real and
1492 imaginary parts of the complex number pointed to by ZP. For
1495 GSL_SET_COMPLEX(&z, 3, 4)
1497 sets Z to be 3 + 4i.
1499 -- Macro: GSL_SET_REAL (ZP,X)
1500 -- Macro: GSL_SET_IMAG (ZP,Y)
1501 These macros allow the real and imaginary parts of the complex
1502 number pointed to by ZP to be set independently.
1505 File: gsl-ref.info, Node: Properties of complex numbers, Next: Complex arithmetic operators, Prev: Complex numbers, Up: Complex Numbers
1507 5.2 Properties of complex numbers
1508 =================================
1510 -- Function: double gsl_complex_arg (gsl_complex Z)
1511 This function returns the argument of the complex number Z,
1512 \arg(z), where -\pi < \arg(z) <= \pi.
1514 -- Function: double gsl_complex_abs (gsl_complex Z)
1515 This function returns the magnitude of the complex number Z, |z|.
1517 -- Function: double gsl_complex_abs2 (gsl_complex Z)
1518 This function returns the squared magnitude of the complex number
1521 -- Function: double gsl_complex_logabs (gsl_complex Z)
1522 This function returns the natural logarithm of the magnitude of the
1523 complex number Z, \log|z|. It allows an accurate evaluation of
1524 \log|z| when |z| is close to one. The direct evaluation of
1525 `log(gsl_complex_abs(z))' would lead to a loss of precision in
1529 File: gsl-ref.info, Node: Complex arithmetic operators, Next: Elementary Complex Functions, Prev: Properties of complex numbers, Up: Complex Numbers
1531 5.3 Complex arithmetic operators
1532 ================================
1534 -- Function: gsl_complex gsl_complex_add (gsl_complex A, gsl_complex B)
1535 This function returns the sum of the complex numbers A and B,
1538 -- Function: gsl_complex gsl_complex_sub (gsl_complex A, gsl_complex B)
1539 This function returns the difference of the complex numbers A and
1542 -- Function: gsl_complex gsl_complex_mul (gsl_complex A, gsl_complex B)
1543 This function returns the product of the complex numbers A and B,
1546 -- Function: gsl_complex gsl_complex_div (gsl_complex A, gsl_complex B)
1547 This function returns the quotient of the complex numbers A and B,
1550 -- Function: gsl_complex gsl_complex_add_real (gsl_complex A, double X)
1551 This function returns the sum of the complex number A and the real
1554 -- Function: gsl_complex gsl_complex_sub_real (gsl_complex A, double X)
1555 This function returns the difference of the complex number A and
1556 the real number X, z=a-x.
1558 -- Function: gsl_complex gsl_complex_mul_real (gsl_complex A, double X)
1559 This function returns the product of the complex number A and the
1560 real number X, z=ax.
1562 -- Function: gsl_complex gsl_complex_div_real (gsl_complex A, double X)
1563 This function returns the quotient of the complex number A and the
1564 real number X, z=a/x.
1566 -- Function: gsl_complex gsl_complex_add_imag (gsl_complex A, double Y)
1567 This function returns the sum of the complex number A and the
1568 imaginary number iY, z=a+iy.
1570 -- Function: gsl_complex gsl_complex_sub_imag (gsl_complex A, double Y)
1571 This function returns the difference of the complex number A and
1572 the imaginary number iY, z=a-iy.
1574 -- Function: gsl_complex gsl_complex_mul_imag (gsl_complex A, double Y)
1575 This function returns the product of the complex number A and the
1576 imaginary number iY, z=a*(iy).
1578 -- Function: gsl_complex gsl_complex_div_imag (gsl_complex A, double Y)
1579 This function returns the quotient of the complex number A and the
1580 imaginary number iY, z=a/(iy).
1582 -- Function: gsl_complex gsl_complex_conjugate (gsl_complex Z)
1583 This function returns the complex conjugate of the complex number
1586 -- Function: gsl_complex gsl_complex_inverse (gsl_complex Z)
1587 This function returns the inverse, or reciprocal, of the complex
1588 number Z, 1/z = (x - i y)/(x^2 + y^2).
1590 -- Function: gsl_complex gsl_complex_negative (gsl_complex Z)
1591 This function returns the negative of the complex number Z, -z =
1595 File: gsl-ref.info, Node: Elementary Complex Functions, Next: Complex Trigonometric Functions, Prev: Complex arithmetic operators, Up: Complex Numbers
1597 5.4 Elementary Complex Functions
1598 ================================
1600 -- Function: gsl_complex gsl_complex_sqrt (gsl_complex Z)
1601 This function returns the square root of the complex number Z,
1602 \sqrt z. The branch cut is the negative real axis. The result
1603 always lies in the right half of the complex plane.
1605 -- Function: gsl_complex gsl_complex_sqrt_real (double X)
1606 This function returns the complex square root of the real number
1607 X, where X may be negative.
1609 -- Function: gsl_complex gsl_complex_pow (gsl_complex Z, gsl_complex A)
1610 The function returns the complex number Z raised to the complex
1611 power A, z^a. This is computed as \exp(\log(z)*a) using complex
1612 logarithms and complex exponentials.
1614 -- Function: gsl_complex gsl_complex_pow_real (gsl_complex Z, double X)
1615 This function returns the complex number Z raised to the real
1618 -- Function: gsl_complex gsl_complex_exp (gsl_complex Z)
1619 This function returns the complex exponential of the complex number
1622 -- Function: gsl_complex gsl_complex_log (gsl_complex Z)
1623 This function returns the complex natural logarithm (base e) of
1624 the complex number Z, \log(z). The branch cut is the negative
1627 -- Function: gsl_complex gsl_complex_log10 (gsl_complex Z)
1628 This function returns the complex base-10 logarithm of the complex
1629 number Z, \log_10 (z).
1631 -- Function: gsl_complex gsl_complex_log_b (gsl_complex Z, gsl_complex
1633 This function returns the complex base-B logarithm of the complex
1634 number Z, \log_b(z). This quantity is computed as the ratio
1638 File: gsl-ref.info, Node: Complex Trigonometric Functions, Next: Inverse Complex Trigonometric Functions, Prev: Elementary Complex Functions, Up: Complex Numbers
1640 5.5 Complex Trigonometric Functions
1641 ===================================
1643 -- Function: gsl_complex gsl_complex_sin (gsl_complex Z)
1644 This function returns the complex sine of the complex number Z,
1645 \sin(z) = (\exp(iz) - \exp(-iz))/(2i).
1647 -- Function: gsl_complex gsl_complex_cos (gsl_complex Z)
1648 This function returns the complex cosine of the complex number Z,
1649 \cos(z) = (\exp(iz) + \exp(-iz))/2.
1651 -- Function: gsl_complex gsl_complex_tan (gsl_complex Z)
1652 This function returns the complex tangent of the complex number Z,
1653 \tan(z) = \sin(z)/\cos(z).
1655 -- Function: gsl_complex gsl_complex_sec (gsl_complex Z)
1656 This function returns the complex secant of the complex number Z,
1657 \sec(z) = 1/\cos(z).
1659 -- Function: gsl_complex gsl_complex_csc (gsl_complex Z)
1660 This function returns the complex cosecant of the complex number Z,
1661 \csc(z) = 1/\sin(z).
1663 -- Function: gsl_complex gsl_complex_cot (gsl_complex Z)
1664 This function returns the complex cotangent of the complex number
1665 Z, \cot(z) = 1/\tan(z).
1668 File: gsl-ref.info, Node: Inverse Complex Trigonometric Functions, Next: Complex Hyperbolic Functions, Prev: Complex Trigonometric Functions, Up: Complex Numbers
1670 5.6 Inverse Complex Trigonometric Functions
1671 ===========================================
1673 -- Function: gsl_complex gsl_complex_arcsin (gsl_complex Z)
1674 This function returns the complex arcsine of the complex number Z,
1675 \arcsin(z). The branch cuts are on the real axis, less than -1 and
1678 -- Function: gsl_complex gsl_complex_arcsin_real (double Z)
1679 This function returns the complex arcsine of the real number Z,
1680 \arcsin(z). For z between -1 and 1, the function returns a real
1681 value in the range [-\pi/2,\pi/2]. For z less than -1 the result
1682 has a real part of -\pi/2 and a positive imaginary part. For z
1683 greater than 1 the result has a real part of \pi/2 and a negative
1686 -- Function: gsl_complex gsl_complex_arccos (gsl_complex Z)
1687 This function returns the complex arccosine of the complex number
1688 Z, \arccos(z). The branch cuts are on the real axis, less than -1
1691 -- Function: gsl_complex gsl_complex_arccos_real (double Z)
1692 This function returns the complex arccosine of the real number Z,
1693 \arccos(z). For z between -1 and 1, the function returns a real
1694 value in the range [0,\pi]. For z less than -1 the result has a
1695 real part of \pi and a negative imaginary part. For z greater
1696 than 1 the result is purely imaginary and positive.
1698 -- Function: gsl_complex gsl_complex_arctan (gsl_complex Z)
1699 This function returns the complex arctangent of the complex number
1700 Z, \arctan(z). The branch cuts are on the imaginary axis, below -i
1703 -- Function: gsl_complex gsl_complex_arcsec (gsl_complex Z)
1704 This function returns the complex arcsecant of the complex number
1705 Z, \arcsec(z) = \arccos(1/z).
1707 -- Function: gsl_complex gsl_complex_arcsec_real (double Z)
1708 This function returns the complex arcsecant of the real number Z,
1709 \arcsec(z) = \arccos(1/z).
1711 -- Function: gsl_complex gsl_complex_arccsc (gsl_complex Z)
1712 This function returns the complex arccosecant of the complex
1713 number Z, \arccsc(z) = \arcsin(1/z).
1715 -- Function: gsl_complex gsl_complex_arccsc_real (double Z)
1716 This function returns the complex arccosecant of the real number Z,
1717 \arccsc(z) = \arcsin(1/z).
1719 -- Function: gsl_complex gsl_complex_arccot (gsl_complex Z)
1720 This function returns the complex arccotangent of the complex
1721 number Z, \arccot(z) = \arctan(1/z).
1724 File: gsl-ref.info, Node: Complex Hyperbolic Functions, Next: Inverse Complex Hyperbolic Functions, Prev: Inverse Complex Trigonometric Functions, Up: Complex Numbers
1726 5.7 Complex Hyperbolic Functions
1727 ================================
1729 -- Function: gsl_complex gsl_complex_sinh (gsl_complex Z)
1730 This function returns the complex hyperbolic sine of the complex
1731 number Z, \sinh(z) = (\exp(z) - \exp(-z))/2.
1733 -- Function: gsl_complex gsl_complex_cosh (gsl_complex Z)
1734 This function returns the complex hyperbolic cosine of the complex
1735 number Z, \cosh(z) = (\exp(z) + \exp(-z))/2.
1737 -- Function: gsl_complex gsl_complex_tanh (gsl_complex Z)
1738 This function returns the complex hyperbolic tangent of the
1739 complex number Z, \tanh(z) = \sinh(z)/\cosh(z).
1741 -- Function: gsl_complex gsl_complex_sech (gsl_complex Z)
1742 This function returns the complex hyperbolic secant of the complex
1743 number Z, \sech(z) = 1/\cosh(z).
1745 -- Function: gsl_complex gsl_complex_csch (gsl_complex Z)
1746 This function returns the complex hyperbolic cosecant of the
1747 complex number Z, \csch(z) = 1/\sinh(z).
1749 -- Function: gsl_complex gsl_complex_coth (gsl_complex Z)
1750 This function returns the complex hyperbolic cotangent of the
1751 complex number Z, \coth(z) = 1/\tanh(z).
1754 File: gsl-ref.info, Node: Inverse Complex Hyperbolic Functions, Next: Complex Number References and Further Reading, Prev: Complex Hyperbolic Functions, Up: Complex Numbers
1756 5.8 Inverse Complex Hyperbolic Functions
1757 ========================================
1759 -- Function: gsl_complex gsl_complex_arcsinh (gsl_complex Z)
1760 This function returns the complex hyperbolic arcsine of the
1761 complex number Z, \arcsinh(z). The branch cuts are on the
1762 imaginary axis, below -i and above i.
1764 -- Function: gsl_complex gsl_complex_arccosh (gsl_complex Z)
1765 This function returns the complex hyperbolic arccosine of the
1766 complex number Z, \arccosh(z). The branch cut is on the real
1767 axis, less than 1. Note that in this case we use the negative
1768 square root in formula 4.6.21 of Abramowitz & Stegun giving
1769 \arccosh(z)=\log(z-\sqrt{z^2-1}).
1771 -- Function: gsl_complex gsl_complex_arccosh_real (double Z)
1772 This function returns the complex hyperbolic arccosine of the real
1773 number Z, \arccosh(z).
1775 -- Function: gsl_complex gsl_complex_arctanh (gsl_complex Z)
1776 This function returns the complex hyperbolic arctangent of the
1777 complex number Z, \arctanh(z). The branch cuts are on the real
1778 axis, less than -1 and greater than 1.
1780 -- Function: gsl_complex gsl_complex_arctanh_real (double Z)
1781 This function returns the complex hyperbolic arctangent of the real
1782 number Z, \arctanh(z).
1784 -- Function: gsl_complex gsl_complex_arcsech (gsl_complex Z)
1785 This function returns the complex hyperbolic arcsecant of the
1786 complex number Z, \arcsech(z) = \arccosh(1/z).
1788 -- Function: gsl_complex gsl_complex_arccsch (gsl_complex Z)
1789 This function returns the complex hyperbolic arccosecant of the
1790 complex number Z, \arccsch(z) = \arcsin(1/z).
1792 -- Function: gsl_complex gsl_complex_arccoth (gsl_complex Z)
1793 This function returns the complex hyperbolic arccotangent of the
1794 complex number Z, \arccoth(z) = \arctanh(1/z).
1797 File: gsl-ref.info, Node: Complex Number References and Further Reading, Prev: Inverse Complex Hyperbolic Functions, Up: Complex Numbers
1799 5.9 References and Further Reading
1800 ==================================
1802 The implementations of the elementary and trigonometric functions are
1803 based on the following papers,
1805 T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang,
1806 "Implementing Complex Elementary Functions Using Exception
1807 Handling", `ACM Transactions on Mathematical Software', Volume 20
1808 (1994), pp 215-244, Corrigenda, p553
1810 T. E. Hull, Thomas F. Fairgrieve, Ping Tak Peter Tang,
1811 "Implementing the complex arcsin and arccosine functions using
1812 exception handling", `ACM Transactions on Mathematical Software',
1813 Volume 23 (1997) pp 299-335
1815 The general formulas and details of branch cuts can be found in the
1818 Abramowitz and Stegun, `Handbook of Mathematical Functions',
1819 "Circular Functions in Terms of Real and Imaginary Parts", Formulas
1820 4.3.55-58, "Inverse Circular Functions in Terms of Real and
1821 Imaginary Parts", Formulas 4.4.37-39, "Hyperbolic Functions in
1822 Terms of Real and Imaginary Parts", Formulas 4.5.49-52, "Inverse
1823 Hyperbolic Functions--relation to Inverse Circular Functions",
1826 Dave Gillespie, `Calc Manual', Free Software Foundation, ISBN
1830 File: gsl-ref.info, Node: Polynomials, Next: Special Functions, Prev: Complex Numbers, Up: Top
1835 This chapter describes functions for evaluating and solving polynomials.
1836 There are routines for finding real and complex roots of quadratic and
1837 cubic equations using analytic methods. An iterative polynomial solver
1838 is also available for finding the roots of general polynomials with real
1839 coefficients (of any order). The functions are declared in the header
1844 * Polynomial Evaluation::
1845 * Divided Difference Representation of Polynomials::
1846 * Quadratic Equations::
1848 * General Polynomial Equations::
1849 * Roots of Polynomials Examples::
1850 * Roots of Polynomials References and Further Reading::
1853 File: gsl-ref.info, Node: Polynomial Evaluation, Next: Divided Difference Representation of Polynomials, Up: Polynomials
1855 6.1 Polynomial Evaluation
1856 =========================
1858 The functions described here evaluate the polynomial c[0] + c[1] x +
1859 c[2] x^2 + \dots + c[len-1] x^{len-1} using Horner's method for
1860 stability. Inline versions of these functions are used when
1861 `HAVE_INLINE' is defined.
1863 -- Function: double gsl_poly_eval (const double C[], const int LEN,
1865 This function evaluates a polynomial with real coefficients for
1866 the real variable X.
1868 -- Function: gsl_complex gsl_poly_complex_eval (const double C[],
1869 const int LEN, const gsl_complex Z)
1870 This function evaluates a polynomial with real coefficients for
1871 the complex variable Z.
1873 -- Function: gsl_complex gsl_complex_poly_complex_eval (const
1874 gsl_complex C[], const int LEN, const gsl_complex Z)
1875 This function evaluates a polynomial with complex coefficients for
1876 the complex variable Z.
1879 File: gsl-ref.info, Node: Divided Difference Representation of Polynomials, Next: Quadratic Equations, Prev: Polynomial Evaluation, Up: Polynomials
1881 6.2 Divided Difference Representation of Polynomials
1882 ====================================================
1884 The functions described here manipulate polynomials stored in Newton's
1885 divided-difference representation. The use of divided-differences is
1886 described in Abramowitz & Stegun sections 25.1.4 and 25.2.26.
1888 -- Function: int gsl_poly_dd_init (double DD[], const double XA[],
1889 const double YA[], size_t SIZE)
1890 This function computes a divided-difference representation of the
1891 interpolating polynomial for the points (XA, YA) stored in the
1892 arrays XA and YA of length SIZE. On output the
1893 divided-differences of (XA,YA) are stored in the array DD, also of
1896 -- Function: double gsl_poly_dd_eval (const double DD[], const double
1897 XA[], const size_t SIZE, const double X)
1898 This function evaluates the polynomial stored in
1899 divided-difference form in the arrays DD and XA of length SIZE at
1900 the point X. An inline version of this function is used when
1901 `HAVE_INLINE' is defined.
1903 -- Function: int gsl_poly_dd_taylor (double C[], double XP, const
1904 double DD[], const double XA[], size_t SIZE, double W[])
1905 This function converts the divided-difference representation of a
1906 polynomial to a Taylor expansion. The divided-difference
1907 representation is supplied in the arrays DD and XA of length SIZE.
1908 On output the Taylor coefficients of the polynomial expanded about
1909 the point XP are stored in the array C also of length SIZE. A
1910 workspace of length SIZE must be provided in the array W.
1913 File: gsl-ref.info, Node: Quadratic Equations, Next: Cubic Equations, Prev: Divided Difference Representation of Polynomials, Up: Polynomials
1915 6.3 Quadratic Equations
1916 =======================
1918 -- Function: int gsl_poly_solve_quadratic (double A, double B, double
1919 C, double * X0, double * X1)
1920 This function finds the real roots of the quadratic equation,
1924 The number of real roots (either zero, one or two) is returned, and
1925 their locations are stored in X0 and X1. If no real roots are
1926 found then X0 and X1 are not modified. If one real root is found
1927 (i.e. if a=0) then it is stored in X0. When two real roots are
1928 found they are stored in X0 and X1 in ascending order. The case
1929 of coincident roots is not considered special. For example
1930 (x-1)^2=0 will have two roots, which happen to have exactly equal
1933 The number of roots found depends on the sign of the discriminant
1934 b^2 - 4 a c. This will be subject to rounding and cancellation
1935 errors when computed in double precision, and will also be subject
1936 to errors if the coefficients of the polynomial are inexact.
1937 These errors may cause a discrete change in the number of roots.
1938 However, for polynomials with small integer coefficients the
1939 discriminant can always be computed exactly.
1942 -- Function: int gsl_poly_complex_solve_quadratic (double A, double B,
1943 double C, gsl_complex * Z0, gsl_complex * Z1)
1944 This function finds the complex roots of the quadratic equation,
1948 The number of complex roots is returned (either one or two) and the
1949 locations of the roots are stored in Z0 and Z1. The roots are
1950 returned in ascending order, sorted first by their real components
1951 and then by their imaginary components. If only one real root is
1952 found (i.e. if a=0) then it is stored in Z0.
1956 File: gsl-ref.info, Node: Cubic Equations, Next: General Polynomial Equations, Prev: Quadratic Equations, Up: Polynomials
1961 -- Function: int gsl_poly_solve_cubic (double A, double B, double C,
1962 double * X0, double * X1, double * X2)
1963 This function finds the real roots of the cubic equation,
1965 x^3 + a x^2 + b x + c = 0
1967 with a leading coefficient of unity. The number of real roots
1968 (either one or three) is returned, and their locations are stored
1969 in X0, X1 and X2. If one real root is found then only X0 is
1970 modified. When three real roots are found they are stored in X0,
1971 X1 and X2 in ascending order. The case of coincident roots is not
1972 considered special. For example, the equation (x-1)^3=0 will have
1973 three roots with exactly equal values. As in the quadratic case,
1974 finite precision may cause equal or closely-spaced real roots to
1975 move off the real axis into the complex plane, leading to a
1976 discrete change in the number of real roots.
1978 -- Function: int gsl_poly_complex_solve_cubic (double A, double B,
1979 double C, gsl_complex * Z0, gsl_complex * Z1, gsl_complex *
1981 This function finds the complex roots of the cubic equation,
1983 z^3 + a z^2 + b z + c = 0
1985 The number of complex roots is returned (always three) and the
1986 locations of the roots are stored in Z0, Z1 and Z2. The roots are
1987 returned in ascending order, sorted first by their real components
1988 and then by their imaginary components.
1992 File: gsl-ref.info, Node: General Polynomial Equations, Next: Roots of Polynomials Examples, Prev: Cubic Equations, Up: Polynomials
1994 6.5 General Polynomial Equations
1995 ================================
1997 The roots of polynomial equations cannot be found analytically beyond
1998 the special cases of the quadratic, cubic and quartic equation. The
1999 algorithm described in this section uses an iterative method to find the
2000 approximate locations of roots of higher order polynomials.
2002 -- Function: gsl_poly_complex_workspace *
2003 gsl_poly_complex_workspace_alloc (size_t N)
2004 This function allocates space for a `gsl_poly_complex_workspace'
2005 struct and a workspace suitable for solving a polynomial with N
2006 coefficients using the routine `gsl_poly_complex_solve'.
2008 The function returns a pointer to the newly allocated
2009 `gsl_poly_complex_workspace' if no errors were detected, and a null
2010 pointer in the case of error.
2012 -- Function: void gsl_poly_complex_workspace_free
2013 (gsl_poly_complex_workspace * W)
2014 This function frees all the memory associated with the workspace W.
2016 -- Function: int gsl_poly_complex_solve (const double * A, size_t N,
2017 gsl_poly_complex_workspace * W, gsl_complex_packed_ptr Z)
2018 This function computes the roots of the general polynomial P(x) =
2019 a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} using balanced-QR
2020 reduction of the companion matrix. The parameter N specifies the
2021 length of the coefficient array. The coefficient of the highest
2022 order term must be non-zero. The function requires a workspace W
2023 of the appropriate size. The n-1 roots are returned in the packed
2024 complex array Z of length 2(n-1), alternating real and imaginary
2027 The function returns `GSL_SUCCESS' if all the roots are found. If
2028 the QR reduction does not converge, the error handler is invoked
2029 with an error code of `GSL_EFAILED'. Note that due to finite
2030 precision, roots of higher multiplicity are returned as a cluster
2031 of simple roots with reduced accuracy. The solution of
2032 polynomials with higher-order roots requires specialized
2033 algorithms that take the multiplicity structure into account (see
2034 e.g. Z. Zeng, Algorithm 835, ACM Transactions on Mathematical
2035 Software, Volume 30, Issue 2 (2004), pp 218-236).
2038 File: gsl-ref.info, Node: Roots of Polynomials Examples, Next: Roots of Polynomials References and Further Reading, Prev: General Polynomial Equations, Up: Polynomials
2043 To demonstrate the use of the general polynomial solver we will take the
2044 polynomial P(x) = x^5 - 1 which has the following roots,
2046 1, e^{2\pi i /5}, e^{4\pi i /5}, e^{6\pi i /5}, e^{8\pi i /5}
2048 The following program will find these roots.
2051 #include <gsl/gsl_poly.h>
2057 /* coefficients of P(x) = -1 + x^5 */
2058 double a[6] = { -1, 0, 0, 0, 0, 1 };
2061 gsl_poly_complex_workspace * w
2062 = gsl_poly_complex_workspace_alloc (6);
2064 gsl_poly_complex_solve (a, 6, w, z);
2066 gsl_poly_complex_workspace_free (w);
2068 for (i = 0; i < 5; i++)
2070 printf ("z%d = %+.18f %+.18f\n",
2071 i, z[2*i], z[2*i+1]);
2077 The output of the program is,
2080 z0 = -0.809016994374947451 +0.587785252292473137
2081 z1 = -0.809016994374947451 -0.587785252292473137
2082 z2 = +0.309016994374947451 +0.951056516295153642
2083 z3 = +0.309016994374947451 -0.951056516295153642
2084 z4 = +1.000000000000000000 +0.000000000000000000
2086 which agrees with the analytic result, z_n = \exp(2 \pi n i/5).
2089 File: gsl-ref.info, Node: Roots of Polynomials References and Further Reading, Prev: Roots of Polynomials Examples, Up: Polynomials
2091 6.7 References and Further Reading
2092 ==================================
2094 The balanced-QR method and its error analysis are described in the
2097 R.S. Martin, G. Peters and J.H. Wilkinson, "The QR Algorithm for
2098 Real Hessenberg Matrices", `Numerische Mathematik', 14 (1970),
2101 B.N. Parlett and C. Reinsch, "Balancing a Matrix for Calculation of
2102 Eigenvalues and Eigenvectors", `Numerische Mathematik', 13 (1969),
2105 A. Edelman and H. Murakami, "Polynomial roots from companion matrix
2106 eigenvalues", `Mathematics of Computation', Vol. 64, No. 210
2109 The formulas for divided differences are given in Abramowitz and Stegun,
2111 Abramowitz and Stegun, `Handbook of Mathematical Functions',
2112 Sections 25.1.4 and 25.2.26.
2115 File: gsl-ref.info, Node: Special Functions, Next: Vectors and Matrices, Prev: Polynomials, Up: Top
2120 This chapter describes the GSL special function library. The library
2121 includes routines for calculating the values of Airy functions, Bessel
2122 functions, Clausen functions, Coulomb wave functions, Coupling
2123 coefficients, the Dawson function, Debye functions, Dilogarithms,
2124 Elliptic integrals, Jacobi elliptic functions, Error functions,
2125 Exponential integrals, Fermi-Dirac functions, Gamma functions,
2126 Gegenbauer functions, Hypergeometric functions, Laguerre functions,
2127 Legendre functions and Spherical Harmonics, the Psi (Digamma) Function,
2128 Synchrotron functions, Transport functions, Trigonometric functions and
2129 Zeta functions. Each routine also computes an estimate of the numerical
2130 error in the calculated value of the function.
2132 The functions in this chapter are declared in individual header
2133 files, such as `gsl_sf_airy.h', `gsl_sf_bessel.h', etc. The complete
2134 set of header files can be included using the file `gsl_sf.h'.
2138 * Special Function Usage::
2139 * The gsl_sf_result struct::
2140 * Special Function Modes::
2141 * Airy Functions and Derivatives::
2142 * Bessel Functions::
2143 * Clausen Functions::
2144 * Coulomb Functions::
2145 * Coupling Coefficients::
2149 * Elementary Operations::
2150 * Elliptic Integrals::
2151 * Elliptic Functions (Jacobi)::
2153 * Exponential Functions::
2154 * Exponential Integrals::
2155 * Fermi-Dirac Function::
2156 * Gamma and Beta Functions::
2157 * Gegenbauer Functions::
2158 * Hypergeometric Functions::
2159 * Laguerre Functions::
2160 * Lambert W Functions::
2161 * Legendre Functions and Spherical Harmonics::
2162 * Logarithm and Related Functions::
2163 * Mathieu Functions::
2165 * Psi (Digamma) Function::
2166 * Synchrotron Functions::
2167 * Transport Functions::
2168 * Trigonometric Functions::
2170 * Special Functions Examples::
2171 * Special Functions References and Further Reading::
2174 File: gsl-ref.info, Node: Special Function Usage, Next: The gsl_sf_result struct, Up: Special Functions
2179 The special functions are available in two calling conventions, a
2180 "natural form" which returns the numerical value of the function and an
2181 "error-handling form" which returns an error code. The two types of
2182 function provide alternative ways of accessing the same underlying code.
2184 The "natural form" returns only the value of the function and can be
2185 used directly in mathematical expressions. For example, the following
2186 function call will compute the value of the Bessel function J_0(x),
2188 double y = gsl_sf_bessel_J0 (x);
2190 There is no way to access an error code or to estimate the error using
2191 this method. To allow access to this information the alternative
2192 error-handling form stores the value and error in a modifiable argument,
2194 gsl_sf_result result;
2195 int status = gsl_sf_bessel_J0_e (x, &result);
2197 The error-handling functions have the suffix `_e'. The returned status
2198 value indicates error conditions such as overflow, underflow or loss of
2199 precision. If there are no errors the error-handling functions return
2203 File: gsl-ref.info, Node: The gsl_sf_result struct, Next: Special Function Modes, Prev: Special Function Usage, Up: Special Functions
2205 7.2 The gsl_sf_result struct
2206 ============================
2208 The error handling form of the special functions always calculate an
2209 error estimate along with the value of the result. Therefore,
2210 structures are provided for amalgamating a value and error estimate.
2211 These structures are declared in the header file `gsl_sf_result.h'.
2213 The `gsl_sf_result' struct contains value and error fields.
2221 The field VAL contains the value and the field ERR contains an estimate
2222 of the absolute error in the value.
2224 In some cases, an overflow or underflow can be detected and handled
2225 by a function. In this case, it may be possible to return a scaling
2226 exponent as well as an error/value pair in order to save the result from
2227 exceeding the dynamic range of the built-in types. The
2228 `gsl_sf_result_e10' struct contains value and error fields as well as
2229 an exponent field such that the actual result is obtained as `result *
2237 } gsl_sf_result_e10;
2240 File: gsl-ref.info, Node: Special Function Modes, Next: Airy Functions and Derivatives, Prev: The gsl_sf_result struct, Up: Special Functions
2245 The goal of the library is to achieve double precision accuracy wherever
2246 possible. However the cost of evaluating some special functions to
2247 double precision can be significant, particularly where very high order
2248 terms are required. In these cases a `mode' argument allows the
2249 accuracy of the function to be reduced in order to improve performance.
2250 The following precision levels are available for the mode argument,
2253 Double-precision, a relative accuracy of approximately 2 * 10^-16.
2256 Single-precision, a relative accuracy of approximately 10^-7.
2259 Approximate values, a relative accuracy of approximately 5 * 10^-4.
2261 The approximate mode provides the fastest evaluation at the lowest
2265 File: gsl-ref.info, Node: Airy Functions and Derivatives, Next: Bessel Functions, Prev: Special Function Modes, Up: Special Functions
2267 7.4 Airy Functions and Derivatives
2268 ==================================
2270 The Airy functions Ai(x) and Bi(x) are defined by the integral
2273 Ai(x) = (1/\pi) \int_0^\infty \cos((1/3) t^3 + xt) dt
2274 Bi(x) = (1/\pi) \int_0^\infty (e^(-(1/3) t^3) + \sin((1/3) t^3 + xt)) dt
2276 For further information see Abramowitz & Stegun, Section 10.4. The Airy
2277 functions are defined in the header file `gsl_sf_airy.h'.
2282 * Derivatives of Airy Functions::
2283 * Zeros of Airy Functions::
2284 * Zeros of Derivatives of Airy Functions::
2287 File: gsl-ref.info, Node: Airy Functions, Next: Derivatives of Airy Functions, Up: Airy Functions and Derivatives
2289 7.4.1 Airy Functions
2290 --------------------
2292 -- Function: double gsl_sf_airy_Ai (double X, gsl_mode_t MODE)
2293 -- Function: int gsl_sf_airy_Ai_e (double X, gsl_mode_t MODE,
2294 gsl_sf_result * RESULT)
2295 These routines compute the Airy function Ai(x) with an accuracy
2298 -- Function: double gsl_sf_airy_Bi (double X, gsl_mode_t MODE)
2299 -- Function: int gsl_sf_airy_Bi_e (double X, gsl_mode_t MODE,
2300 gsl_sf_result * RESULT)
2301 These routines compute the Airy function Bi(x) with an accuracy
2304 -- Function: double gsl_sf_airy_Ai_scaled (double X, gsl_mode_t MODE)
2305 -- Function: int gsl_sf_airy_Ai_scaled_e (double X, gsl_mode_t MODE,
2306 gsl_sf_result * RESULT)
2307 These routines compute a scaled version of the Airy function
2308 S_A(x) Ai(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3)
2309 x^(3/2)), and is 1 for x<0.
2311 -- Function: double gsl_sf_airy_Bi_scaled (double X, gsl_mode_t MODE)
2312 -- Function: int gsl_sf_airy_Bi_scaled_e (double X, gsl_mode_t MODE,
2313 gsl_sf_result * RESULT)
2314 These routines compute a scaled version of the Airy function
2315 S_B(x) Bi(x). For x>0 the scaling factor S_B(x) is exp(-(2/3)
2316 x^(3/2)), and is 1 for x<0.
2319 File: gsl-ref.info, Node: Derivatives of Airy Functions, Next: Zeros of Airy Functions, Prev: Airy Functions, Up: Airy Functions and Derivatives
2321 7.4.2 Derivatives of Airy Functions
2322 -----------------------------------
2324 -- Function: double gsl_sf_airy_Ai_deriv (double X, gsl_mode_t MODE)
2325 -- Function: int gsl_sf_airy_Ai_deriv_e (double X, gsl_mode_t MODE,
2326 gsl_sf_result * RESULT)
2327 These routines compute the Airy function derivative Ai'(x) with an
2328 accuracy specified by MODE.
2330 -- Function: double gsl_sf_airy_Bi_deriv (double X, gsl_mode_t MODE)
2331 -- Function: int gsl_sf_airy_Bi_deriv_e (double X, gsl_mode_t MODE,
2332 gsl_sf_result * RESULT)
2333 These routines compute the Airy function derivative Bi'(x) with an
2334 accuracy specified by MODE.
2336 -- Function: double gsl_sf_airy_Ai_deriv_scaled (double X, gsl_mode_t
2338 -- Function: int gsl_sf_airy_Ai_deriv_scaled_e (double X, gsl_mode_t
2339 MODE, gsl_sf_result * RESULT)
2340 These routines compute the scaled Airy function derivative S_A(x)
2341 Ai'(x). For x>0 the scaling factor S_A(x) is \exp(+(2/3)
2342 x^(3/2)), and is 1 for x<0.
2344 -- Function: double gsl_sf_airy_Bi_deriv_scaled (double X, gsl_mode_t
2346 -- Function: int gsl_sf_airy_Bi_deriv_scaled_e (double X, gsl_mode_t
2347 MODE, gsl_sf_result * RESULT)
2348 These routines compute the scaled Airy function derivative S_B(x)
2349 Bi'(x). For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)),
2353 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
2355 7.4.3 Zeros of Airy Functions
2356 -----------------------------
2358 -- Function: double gsl_sf_airy_zero_Ai (unsigned int S)
2359 -- Function: int gsl_sf_airy_zero_Ai_e (unsigned int S, gsl_sf_result
2361 These routines compute the location of the S-th zero of the Airy
2364 -- Function: double gsl_sf_airy_zero_Bi (unsigned int S)
2365 -- Function: int gsl_sf_airy_zero_Bi_e (unsigned int S, gsl_sf_result
2367 These routines compute the location of the S-th zero of the Airy
2371 File: gsl-ref.info, Node: Zeros of Derivatives of Airy Functions, Prev: Zeros of Airy Functions, Up: Airy Functions and Derivatives
2373 7.4.4 Zeros of Derivatives of Airy Functions
2374 --------------------------------------------
2376 -- Function: double gsl_sf_airy_zero_Ai_deriv (unsigned int S)
2377 -- Function: int gsl_sf_airy_zero_Ai_deriv_e (unsigned int S,
2378 gsl_sf_result * RESULT)
2379 These routines compute the location of the S-th zero of the Airy
2380 function derivative Ai'(x).
2382 -- Function: double gsl_sf_airy_zero_Bi_deriv (unsigned int S)
2383 -- Function: int gsl_sf_airy_zero_Bi_deriv_e (unsigned int S,
2384 gsl_sf_result * RESULT)
2385 These routines compute the location of the S-th zero of the Airy
2386 function derivative Bi'(x).
2389 File: gsl-ref.info, Node: Bessel Functions, Next: Clausen Functions, Prev: Airy Functions and Derivatives, Up: Special Functions
2391 7.5 Bessel Functions
2392 ====================
2394 The routines described in this section compute the Cylindrical Bessel
2395 functions J_n(x), Y_n(x), Modified cylindrical Bessel functions I_n(x),
2396 K_n(x), Spherical Bessel functions j_l(x), y_l(x), and Modified
2397 Spherical Bessel functions i_l(x), k_l(x). For more information see
2398 Abramowitz & Stegun, Chapters 9 and 10. The Bessel functions are
2399 defined in the header file `gsl_sf_bessel.h'.
2403 * Regular Cylindrical Bessel Functions::
2404 * Irregular Cylindrical Bessel Functions::
2405 * Regular Modified Cylindrical Bessel Functions::
2406 * Irregular Modified Cylindrical Bessel Functions::
2407 * Regular Spherical Bessel Functions::
2408 * Irregular Spherical Bessel Functions::
2409 * Regular Modified Spherical Bessel Functions::
2410 * Irregular Modified Spherical Bessel Functions::
2411 * Regular Bessel Function - Fractional Order::
2412 * Irregular Bessel Functions - Fractional Order::
2413 * Regular Modified Bessel Functions - Fractional Order::
2414 * Irregular Modified Bessel Functions - Fractional Order::
2415 * Zeros of Regular Bessel Functions::
2418 File: gsl-ref.info, Node: Regular Cylindrical Bessel Functions, Next: Irregular Cylindrical Bessel Functions, Up: Bessel Functions
2420 7.5.1 Regular Cylindrical Bessel Functions
2421 ------------------------------------------
2423 -- Function: double gsl_sf_bessel_J0 (double X)
2424 -- Function: int gsl_sf_bessel_J0_e (double X, gsl_sf_result * RESULT)
2425 These routines compute the regular cylindrical Bessel function of
2426 zeroth order, J_0(x).
2428 -- Function: double gsl_sf_bessel_J1 (double X)
2429 -- Function: int gsl_sf_bessel_J1_e (double X, gsl_sf_result * RESULT)
2430 These routines compute the regular cylindrical Bessel function of
2431 first order, J_1(x).
2433 -- Function: double gsl_sf_bessel_Jn (int N, double X)
2434 -- Function: int gsl_sf_bessel_Jn_e (int N, double X, gsl_sf_result *
2436 These routines compute the regular cylindrical Bessel function of
2439 -- Function: int gsl_sf_bessel_Jn_array (int NMIN, int NMAX, double X,
2440 double RESULT_ARRAY[])
2441 This routine computes the values of the regular cylindrical Bessel
2442 functions J_n(x) for n from NMIN to NMAX inclusive, storing the
2443 results in the array RESULT_ARRAY. The values are computed using
2444 recurrence relations for efficiency, and therefore may differ
2445 slightly from the exact values.
2448 File: gsl-ref.info, Node: Irregular Cylindrical Bessel Functions, Next: Regular Modified Cylindrical Bessel Functions, Prev: Regular Cylindrical Bessel Functions, Up: Bessel Functions
2450 7.5.2 Irregular Cylindrical Bessel Functions
2451 --------------------------------------------
2453 -- Function: double gsl_sf_bessel_Y0 (double X)
2454 -- Function: int gsl_sf_bessel_Y0_e (double X, gsl_sf_result * RESULT)
2455 These routines compute the irregular cylindrical Bessel function
2456 of zeroth order, Y_0(x), for x>0.
2458 -- Function: double gsl_sf_bessel_Y1 (double X)
2459 -- Function: int gsl_sf_bessel_Y1_e (double X, gsl_sf_result * RESULT)
2460 These routines compute the irregular cylindrical Bessel function
2461 of first order, Y_1(x), for x>0.
2463 -- Function: double gsl_sf_bessel_Yn (int N,double X)
2464 -- Function: int gsl_sf_bessel_Yn_e (int N,double X, gsl_sf_result *
2466 These routines compute the irregular cylindrical Bessel function of
2467 order N, Y_n(x), for x>0.
2469 -- Function: int gsl_sf_bessel_Yn_array (int NMIN, int NMAX, double X,
2470 double RESULT_ARRAY[])
2471 This routine computes the values of the irregular cylindrical
2472 Bessel functions Y_n(x) for n from NMIN to NMAX inclusive, storing
2473 the results in the array RESULT_ARRAY. The domain of the function
2474 is x>0. The values are computed using recurrence relations for
2475 efficiency, and therefore may differ slightly from the exact
2479 File: gsl-ref.info, Node: Regular Modified Cylindrical Bessel Functions, Next: Irregular Modified Cylindrical Bessel Functions, Prev: Irregular Cylindrical Bessel Functions, Up: Bessel Functions
2481 7.5.3 Regular Modified Cylindrical Bessel Functions
2482 ---------------------------------------------------
2484 -- Function: double gsl_sf_bessel_I0 (double X)
2485 -- Function: int gsl_sf_bessel_I0_e (double X, gsl_sf_result * RESULT)
2486 These routines compute the regular modified cylindrical Bessel
2487 function of zeroth order, I_0(x).
2489 -- Function: double gsl_sf_bessel_I1 (double X)
2490 -- Function: int gsl_sf_bessel_I1_e (double X, gsl_sf_result * RESULT)
2491 These routines compute the regular modified cylindrical Bessel
2492 function of first order, I_1(x).
2494 -- Function: double gsl_sf_bessel_In (int N, double X)
2495 -- Function: int gsl_sf_bessel_In_e (int N, double X, gsl_sf_result *
2497 These routines compute the regular modified cylindrical Bessel
2498 function of order N, I_n(x).
2500 -- Function: int gsl_sf_bessel_In_array (int NMIN, int NMAX, double X,
2501 double RESULT_ARRAY[])
2502 This routine computes the values of the regular modified
2503 cylindrical Bessel functions I_n(x) for n from NMIN to NMAX
2504 inclusive, storing the results in the array RESULT_ARRAY. The
2505 start of the range NMIN must be positive or zero. The values are
2506 computed using recurrence relations for efficiency, and therefore
2507 may differ slightly from the exact values.
2509 -- Function: double gsl_sf_bessel_I0_scaled (double X)
2510 -- Function: int gsl_sf_bessel_I0_scaled_e (double X, gsl_sf_result *
2512 These routines compute the scaled regular modified cylindrical
2513 Bessel function of zeroth order \exp(-|x|) I_0(x).
2515 -- Function: double gsl_sf_bessel_I1_scaled (double X)
2516 -- Function: int gsl_sf_bessel_I1_scaled_e (double X, gsl_sf_result *
2518 These routines compute the scaled regular modified cylindrical
2519 Bessel function of first order \exp(-|x|) I_1(x).
2521 -- Function: double gsl_sf_bessel_In_scaled (int N, double X)
2522 -- Function: int gsl_sf_bessel_In_scaled_e (int N, double X,
2523 gsl_sf_result * RESULT)
2524 These routines compute the scaled regular modified cylindrical
2525 Bessel function of order N, \exp(-|x|) I_n(x)
2527 -- Function: int gsl_sf_bessel_In_scaled_array (int NMIN, int NMAX,
2528 double X, double RESULT_ARRAY[])
2529 This routine computes the values of the scaled regular cylindrical
2530 Bessel functions \exp(-|x|) I_n(x) for n from NMIN to NMAX
2531 inclusive, storing the results in the array RESULT_ARRAY. The
2532 start of the range NMIN must be positive or zero. The values are
2533 computed using recurrence relations for efficiency, and therefore
2534 may differ slightly from the exact values.
2537 File: gsl-ref.info, Node: Irregular Modified Cylindrical Bessel Functions, Next: Regular Spherical Bessel Functions, Prev: Regular Modified Cylindrical Bessel Functions, Up: Bessel Functions
2539 7.5.4 Irregular Modified Cylindrical Bessel Functions
2540 -----------------------------------------------------
2542 -- Function: double gsl_sf_bessel_K0 (double X)
2543 -- Function: int gsl_sf_bessel_K0_e (double X, gsl_sf_result * RESULT)
2544 These routines compute the irregular modified cylindrical Bessel
2545 function of zeroth order, K_0(x), for x > 0.
2547 -- Function: double gsl_sf_bessel_K1 (double X)
2548 -- Function: int gsl_sf_bessel_K1_e (double X, gsl_sf_result * RESULT)
2549 These routines compute the irregular modified cylindrical Bessel
2550 function of first order, K_1(x), for x > 0.
2552 -- Function: double gsl_sf_bessel_Kn (int N, double X)
2553 -- Function: int gsl_sf_bessel_Kn_e (int N, double X, gsl_sf_result *
2555 These routines compute the irregular modified cylindrical Bessel
2556 function of order N, K_n(x), for x > 0.
2558 -- Function: int gsl_sf_bessel_Kn_array (int NMIN, int NMAX, double X,
2559 double RESULT_ARRAY[])
2560 This routine computes the values of the irregular modified
2561 cylindrical Bessel functions K_n(x) for n from NMIN to NMAX
2562 inclusive, storing the results in the array RESULT_ARRAY. The
2563 start of the range NMIN must be positive or zero. The domain of
2564 the function is x>0. The values are computed using recurrence
2565 relations for efficiency, and therefore may differ slightly from
2568 -- Function: double gsl_sf_bessel_K0_scaled (double X)
2569 -- Function: int gsl_sf_bessel_K0_scaled_e (double X, gsl_sf_result *
2571 These routines compute the scaled irregular modified cylindrical
2572 Bessel function of zeroth order \exp(x) K_0(x) for x>0.
2574 -- Function: double gsl_sf_bessel_K1_scaled (double X)
2575 -- Function: int gsl_sf_bessel_K1_scaled_e (double X, gsl_sf_result *
2577 These routines compute the scaled irregular modified cylindrical
2578 Bessel function of first order \exp(x) K_1(x) for x>0.
2580 -- Function: double gsl_sf_bessel_Kn_scaled (int N, double X)
2581 -- Function: int gsl_sf_bessel_Kn_scaled_e (int N, double X,
2582 gsl_sf_result * RESULT)
2583 These routines compute the scaled irregular modified cylindrical
2584 Bessel function of order N, \exp(x) K_n(x), for x>0.
2586 -- Function: int gsl_sf_bessel_Kn_scaled_array (int NMIN, int NMAX,
2587 double X, double RESULT_ARRAY[])
2588 This routine computes the values of the scaled irregular
2589 cylindrical Bessel functions \exp(x) K_n(x) for n from NMIN to
2590 NMAX inclusive, storing the results in the array RESULT_ARRAY. The
2591 start of the range NMIN must be positive or zero. The domain of
2592 the function is x>0. The values are computed using recurrence
2593 relations for efficiency, and therefore may differ slightly from
2597 File: gsl-ref.info, Node: Regular Spherical Bessel Functions, Next: Irregular Spherical Bessel Functions, Prev: Irregular Modified Cylindrical Bessel Functions, Up: Bessel Functions
2599 7.5.5 Regular Spherical Bessel Functions
2600 ----------------------------------------
2602 -- Function: double gsl_sf_bessel_j0 (double X)
2603 -- Function: int gsl_sf_bessel_j0_e (double X, gsl_sf_result * RESULT)
2604 These routines compute the regular spherical Bessel function of
2605 zeroth order, j_0(x) = \sin(x)/x.
2607 -- Function: double gsl_sf_bessel_j1 (double X)
2608 -- Function: int gsl_sf_bessel_j1_e (double X, gsl_sf_result * RESULT)
2609 These routines compute the regular spherical Bessel function of
2610 first order, j_1(x) = (\sin(x)/x - \cos(x))/x.
2612 -- Function: double gsl_sf_bessel_j2 (double X)
2613 -- Function: int gsl_sf_bessel_j2_e (double X, gsl_sf_result * RESULT)
2614 These routines compute the regular spherical Bessel function of
2615 second order, j_2(x) = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x.
2617 -- Function: double gsl_sf_bessel_jl (int L, double X)
2618 -- Function: int gsl_sf_bessel_jl_e (int L, double X, gsl_sf_result *
2620 These routines compute the regular spherical Bessel function of
2621 order L, j_l(x), for l >= 0 and x >= 0.
2623 -- Function: int gsl_sf_bessel_jl_array (int LMAX, double X, double
2625 This routine computes the values of the regular spherical Bessel
2626 functions j_l(x) for l from 0 to LMAX inclusive for lmax >= 0 and
2627 x >= 0, storing the results in the array RESULT_ARRAY. The values
2628 are computed using recurrence relations for efficiency, and
2629 therefore may differ slightly from the exact values.
2631 -- Function: int gsl_sf_bessel_jl_steed_array (int LMAX, double X,
2632 double * JL_X_ARRAY)
2633 This routine uses Steed's method to compute the values of the
2634 regular spherical Bessel functions j_l(x) for l from 0 to LMAX
2635 inclusive for lmax >= 0 and x >= 0, storing the results in the
2636 array RESULT_ARRAY. The Steed/Barnett algorithm is described in
2637 `Comp. Phys. Comm.' 21, 297 (1981). Steed's method is more stable
2638 than the recurrence used in the other functions but is also slower.
2641 File: gsl-ref.info, Node: Irregular Spherical Bessel Functions, Next: Regular Modified Spherical Bessel Functions, Prev: Regular Spherical Bessel Functions, Up: Bessel Functions
2643 7.5.6 Irregular Spherical Bessel Functions
2644 ------------------------------------------
2646 -- Function: double gsl_sf_bessel_y0 (double X)
2647 -- Function: int gsl_sf_bessel_y0_e (double X, gsl_sf_result * RESULT)
2648 These routines compute the irregular spherical Bessel function of
2649 zeroth order, y_0(x) = -\cos(x)/x.
2651 -- Function: double gsl_sf_bessel_y1 (double X)
2652 -- Function: int gsl_sf_bessel_y1_e (double X, gsl_sf_result * RESULT)
2653 These routines compute the irregular spherical Bessel function of
2654 first order, y_1(x) = -(\cos(x)/x + \sin(x))/x.
2656 -- Function: double gsl_sf_bessel_y2 (double X)
2657 -- Function: int gsl_sf_bessel_y2_e (double X, gsl_sf_result * RESULT)
2658 These routines compute the irregular spherical Bessel function of
2659 second order, y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x).
2661 -- Function: double gsl_sf_bessel_yl (int L, double X)
2662 -- Function: int gsl_sf_bessel_yl_e (int L, double X, gsl_sf_result *
2664 These routines compute the irregular spherical Bessel function of
2665 order L, y_l(x), for l >= 0.
2667 -- Function: int gsl_sf_bessel_yl_array (int LMAX, double X, double
2669 This routine computes the values of the irregular spherical Bessel
2670 functions y_l(x) for l from 0 to LMAX inclusive for lmax >= 0,
2671 storing the results in the array RESULT_ARRAY. The values are
2672 computed using recurrence relations for efficiency, and therefore
2673 may differ slightly from the exact values.
2676 File: gsl-ref.info, Node: Regular Modified Spherical Bessel Functions, Next: Irregular Modified Spherical Bessel Functions, Prev: Irregular Spherical Bessel Functions, Up: Bessel Functions
2678 7.5.7 Regular Modified Spherical Bessel Functions
2679 -------------------------------------------------
2681 The regular modified spherical Bessel functions i_l(x) are related to
2682 the modified Bessel functions of fractional order, i_l(x) =
2683 \sqrt{\pi/(2x)} I_{l+1/2}(x)
2685 -- Function: double gsl_sf_bessel_i0_scaled (double X)
2686 -- Function: int gsl_sf_bessel_i0_scaled_e (double X, gsl_sf_result *
2688 These routines compute the scaled regular modified spherical Bessel
2689 function of zeroth order, \exp(-|x|) i_0(x).
2691 -- Function: double gsl_sf_bessel_i1_scaled (double X)
2692 -- Function: int gsl_sf_bessel_i1_scaled_e (double X, gsl_sf_result *
2694 These routines compute the scaled regular modified spherical Bessel
2695 function of first order, \exp(-|x|) i_1(x).
2697 -- Function: double gsl_sf_bessel_i2_scaled (double X)
2698 -- Function: int gsl_sf_bessel_i2_scaled_e (double X, gsl_sf_result *
2700 These routines compute the scaled regular modified spherical Bessel
2701 function of second order, \exp(-|x|) i_2(x)
2703 -- Function: double gsl_sf_bessel_il_scaled (int L, double X)
2704 -- Function: int gsl_sf_bessel_il_scaled_e (int L, double X,
2705 gsl_sf_result * RESULT)
2706 These routines compute the scaled regular modified spherical Bessel
2707 function of order L, \exp(-|x|) i_l(x)
2709 -- Function: int gsl_sf_bessel_il_scaled_array (int LMAX, double X,
2710 double RESULT_ARRAY[])
2711 This routine computes the values of the scaled regular modified
2712 cylindrical Bessel functions \exp(-|x|) i_l(x) for l from 0 to
2713 LMAX inclusive for lmax >= 0, storing the results in the array
2714 RESULT_ARRAY. The values are computed using recurrence relations
2715 for efficiency, and therefore may differ slightly from the exact
2719 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
2721 7.5.8 Irregular Modified Spherical Bessel Functions
2722 ---------------------------------------------------
2724 The irregular modified spherical Bessel functions k_l(x) are related to
2725 the irregular modified Bessel functions of fractional order, k_l(x) =
2726 \sqrt{\pi/(2x)} K_{l+1/2}(x).
2728 -- Function: double gsl_sf_bessel_k0_scaled (double X)
2729 -- Function: int gsl_sf_bessel_k0_scaled_e (double X, gsl_sf_result *
2731 These routines compute the scaled irregular modified spherical
2732 Bessel function of zeroth order, \exp(x) k_0(x), for x>0.
2734 -- Function: double gsl_sf_bessel_k1_scaled (double X)
2735 -- Function: int gsl_sf_bessel_k1_scaled_e (double X, gsl_sf_result *
2737 These routines compute the scaled irregular modified spherical
2738 Bessel function of first order, \exp(x) k_1(x), for x>0.
2740 -- Function: double gsl_sf_bessel_k2_scaled (double X)
2741 -- Function: int gsl_sf_bessel_k2_scaled_e (double X, gsl_sf_result *
2743 These routines compute the scaled irregular modified spherical
2744 Bessel function of second order, \exp(x) k_2(x), for x>0.
2746 -- Function: double gsl_sf_bessel_kl_scaled (int L, double X)
2747 -- Function: int gsl_sf_bessel_kl_scaled_e (int L, double X,
2748 gsl_sf_result * RESULT)
2749 These routines compute the scaled irregular modified spherical
2750 Bessel function of order L, \exp(x) k_l(x), for x>0.
2752 -- Function: int gsl_sf_bessel_kl_scaled_array (int LMAX, double X,
2753 double RESULT_ARRAY[])
2754 This routine computes the values of the scaled irregular modified
2755 spherical Bessel functions \exp(x) k_l(x) for l from 0 to LMAX
2756 inclusive for lmax >= 0 and x>0, storing the results in the array
2757 RESULT_ARRAY. The values are computed using recurrence relations
2758 for efficiency, and therefore may differ slightly from the exact
2762 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
2764 7.5.9 Regular Bessel Function--Fractional Order
2765 -----------------------------------------------
2767 -- Function: double gsl_sf_bessel_Jnu (double NU, double X)
2768 -- Function: int gsl_sf_bessel_Jnu_e (double NU, double X,
2769 gsl_sf_result * RESULT)
2770 These routines compute the regular cylindrical Bessel function of
2771 fractional order \nu, J_\nu(x).
2773 -- Function: int gsl_sf_bessel_sequence_Jnu_e (double NU, gsl_mode_t
2774 MODE, size_t SIZE, double V[])
2775 This function computes the regular cylindrical Bessel function of
2776 fractional order \nu, J_\nu(x), evaluated at a series of x values.
2777 The array V of length SIZE contains the x values. They are
2778 assumed to be strictly ordered and positive. The array is
2779 over-written with the values of J_\nu(x_i).
2782 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
2784 7.5.10 Irregular Bessel Functions--Fractional Order
2785 ---------------------------------------------------
2787 -- Function: double gsl_sf_bessel_Ynu (double NU, double X)
2788 -- Function: int gsl_sf_bessel_Ynu_e (double NU, double X,
2789 gsl_sf_result * RESULT)
2790 These routines compute the irregular cylindrical Bessel function of
2791 fractional order \nu, Y_\nu(x).
2794 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
2796 7.5.11 Regular Modified Bessel Functions--Fractional Order
2797 ----------------------------------------------------------
2799 -- Function: double gsl_sf_bessel_Inu (double NU, double X)
2800 -- Function: int gsl_sf_bessel_Inu_e (double NU, double X,
2801 gsl_sf_result * RESULT)
2802 These routines compute the regular modified Bessel function of
2803 fractional order \nu, I_\nu(x) for x>0, \nu>0.
2805 -- Function: double gsl_sf_bessel_Inu_scaled (double NU, double X)
2806 -- Function: int gsl_sf_bessel_Inu_scaled_e (double NU, double X,
2807 gsl_sf_result * RESULT)
2808 These routines compute the scaled regular modified Bessel function
2809 of fractional order \nu, \exp(-|x|)I_\nu(x) for x>0, \nu>0.
2812 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
2814 7.5.12 Irregular Modified Bessel Functions--Fractional Order
2815 ------------------------------------------------------------
2817 -- Function: double gsl_sf_bessel_Knu (double NU, double X)
2818 -- Function: int gsl_sf_bessel_Knu_e (double NU, double X,
2819 gsl_sf_result * RESULT)
2820 These routines compute the irregular modified Bessel function of
2821 fractional order \nu, K_\nu(x) for x>0, \nu>0.
2823 -- Function: double gsl_sf_bessel_lnKnu (double NU, double X)
2824 -- Function: int gsl_sf_bessel_lnKnu_e (double NU, double X,
2825 gsl_sf_result * RESULT)
2826 These routines compute the logarithm of the irregular modified
2827 Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0,
2830 -- Function: double gsl_sf_bessel_Knu_scaled (double NU, double X)
2831 -- Function: int gsl_sf_bessel_Knu_scaled_e (double NU, double X,
2832 gsl_sf_result * RESULT)
2833 These routines compute the scaled irregular modified Bessel
2834 function of fractional order \nu, \exp(+|x|) K_\nu(x) for x>0,
2838 File: gsl-ref.info, Node: Zeros of Regular Bessel Functions, Prev: Irregular Modified Bessel Functions - Fractional Order, Up: Bessel Functions
2840 7.5.13 Zeros of Regular Bessel Functions
2841 ----------------------------------------
2843 -- Function: double gsl_sf_bessel_zero_J0 (unsigned int S)
2844 -- Function: int gsl_sf_bessel_zero_J0_e (unsigned int S,
2845 gsl_sf_result * RESULT)
2846 These routines compute the location of the S-th positive zero of
2847 the Bessel function J_0(x).
2849 -- Function: double gsl_sf_bessel_zero_J1 (unsigned int S)
2850 -- Function: int gsl_sf_bessel_zero_J1_e (unsigned int S,
2851 gsl_sf_result * RESULT)
2852 These routines compute the location of the S-th positive zero of
2853 the Bessel function J_1(x).
2855 -- Function: double gsl_sf_bessel_zero_Jnu (double NU, unsigned int S)
2856 -- Function: int gsl_sf_bessel_zero_Jnu_e (double NU, unsigned int S,
2857 gsl_sf_result * RESULT)
2858 These routines compute the location of the S-th positive zero of
2859 the Bessel function J_\nu(x). The current implementation does not
2860 support negative values of NU.
2863 File: gsl-ref.info, Node: Clausen Functions, Next: Coulomb Functions, Prev: Bessel Functions, Up: Special Functions
2865 7.6 Clausen Functions
2866 =====================
2868 The Clausen function is defined by the following integral,
2870 Cl_2(x) = - \int_0^x dt \log(2 \sin(t/2))
2872 It is related to the dilogarithm by Cl_2(\theta) = \Im
2873 Li_2(\exp(i\theta)). The Clausen functions are declared in the header
2874 file `gsl_sf_clausen.h'.
2876 -- Function: double gsl_sf_clausen (double X)
2877 -- Function: int gsl_sf_clausen_e (double X, gsl_sf_result * RESULT)
2878 These routines compute the Clausen integral Cl_2(x).
2881 File: gsl-ref.info, Node: Coulomb Functions, Next: Coupling Coefficients, Prev: Clausen Functions, Up: Special Functions
2883 7.7 Coulomb Functions
2884 =====================
2886 The prototypes of the Coulomb functions are declared in the header file
2887 `gsl_sf_coulomb.h'. Both bound state and scattering solutions are
2892 * Normalized Hydrogenic Bound States::
2893 * Coulomb Wave Functions::
2894 * Coulomb Wave Function Normalization Constant::
2897 File: gsl-ref.info, Node: Normalized Hydrogenic Bound States, Next: Coulomb Wave Functions, Up: Coulomb Functions
2899 7.7.1 Normalized Hydrogenic Bound States
2900 ----------------------------------------
2902 -- Function: double gsl_sf_hydrogenicR_1 (double Z, double R)
2903 -- Function: int gsl_sf_hydrogenicR_1_e (double Z, double R,
2904 gsl_sf_result * RESULT)
2905 These routines compute the lowest-order normalized hydrogenic bound
2906 state radial wavefunction R_1 := 2Z \sqrt{Z} \exp(-Z r).
2908 -- Function: double gsl_sf_hydrogenicR (int N, int L, double Z, double
2910 -- Function: int gsl_sf_hydrogenicR_e (int N, int L, double Z, double
2911 R, gsl_sf_result * RESULT)
2912 These routines compute the N-th normalized hydrogenic bound state
2913 radial wavefunction,
2915 R_n := 2 (Z^{3/2}/n^2) \sqrt{(n-l-1)!/(n+l)!} \exp(-Z r/n) (2Zr/n)^l
2916 L^{2l+1}_{n-l-1}(2Zr/n).
2918 where L^a_b(x) is the generalized Laguerre polynomial (*note
2919 Laguerre Functions::). The normalization is chosen such that the
2920 wavefunction \psi is given by \psi(n,l,r) = R_n Y_{lm}.
2923 File: gsl-ref.info, Node: Coulomb Wave Functions, Next: Coulomb Wave Function Normalization Constant, Prev: Normalized Hydrogenic Bound States, Up: Coulomb Functions
2925 7.7.2 Coulomb Wave Functions
2926 ----------------------------
2928 The Coulomb wave functions F_L(\eta,x), G_L(\eta,x) are described in
2929 Abramowitz & Stegun, Chapter 14. Because there can be a large dynamic
2930 range of values for these functions, overflows are handled gracefully.
2931 If an overflow occurs, `GSL_EOVRFLW' is signalled and exponent(s) are
2932 returned through the modifiable parameters EXP_F, EXP_G. The full
2933 solution can be reconstructed from the following relations,
2935 F_L(eta,x) = fc[k_L] * exp(exp_F)
2936 G_L(eta,x) = gc[k_L] * exp(exp_G)
2938 F_L'(eta,x) = fcp[k_L] * exp(exp_F)
2939 G_L'(eta,x) = gcp[k_L] * exp(exp_G)
2942 -- Function: int gsl_sf_coulomb_wave_FG_e (double ETA, double X,
2943 double L_F, int K, gsl_sf_result * F, gsl_sf_result * FP,
2944 gsl_sf_result * G, gsl_sf_result * GP, double * EXP_F, double
2946 This function computes the Coulomb wave functions F_L(\eta,x),
2947 G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x),
2948 G'_{L-k}(\eta,x) with respect to x. The parameters are restricted
2949 to L, L-k > -1/2, x > 0 and integer k. Note that L itself is not
2950 restricted to being an integer. The results are stored in the
2951 parameters F, G for the function values and FP, GP for the
2952 derivative values. If an overflow occurs, `GSL_EOVRFLW' is
2953 returned and scaling exponents are stored in the modifiable
2954 parameters EXP_F, EXP_G.
2956 -- Function: int gsl_sf_coulomb_wave_F_array (double L_MIN, int KMAX,
2957 double ETA, double X, double FC_ARRAY[], double * F_EXPONENT)
2958 This function computes the Coulomb wave function F_L(\eta,x) for L
2959 = Lmin \dots Lmin + kmax, storing the results in FC_ARRAY. In the
2960 case of overflow the exponent is stored in F_EXPONENT.
2962 -- Function: int gsl_sf_coulomb_wave_FG_array (double L_MIN, int KMAX,
2963 double ETA, double X, double FC_ARRAY[], double GC_ARRAY[],
2964 double * F_EXPONENT, double * G_EXPONENT)
2965 This function computes the functions F_L(\eta,x), G_L(\eta,x) for
2966 L = Lmin \dots Lmin + kmax storing the results in FC_ARRAY and
2967 GC_ARRAY. In the case of overflow the exponents are stored in
2968 F_EXPONENT and G_EXPONENT.
2970 -- Function: int gsl_sf_coulomb_wave_FGp_array (double L_MIN, int
2971 KMAX, double ETA, double X, double FC_ARRAY[], double
2972 FCP_ARRAY[], double GC_ARRAY[], double GCP_ARRAY[], double *
2973 F_EXPONENT, double * G_EXPONENT)
2974 This function computes the functions F_L(\eta,x), G_L(\eta,x) and
2975 their derivatives F'_L(\eta,x), G'_L(\eta,x) for L = Lmin \dots
2976 Lmin + kmax storing the results in FC_ARRAY, GC_ARRAY, FCP_ARRAY
2977 and GCP_ARRAY. In the case of overflow the exponents are stored
2978 in F_EXPONENT and G_EXPONENT.
2980 -- Function: int gsl_sf_coulomb_wave_sphF_array (double L_MIN, int
2981 KMAX, double ETA, double X, double FC_ARRAY[], double
2983 This function computes the Coulomb wave function divided by the
2984 argument F_L(\eta, x)/x for L = Lmin \dots Lmin + kmax, storing the
2985 results in FC_ARRAY. In the case of overflow the exponent is
2986 stored in F_EXPONENT. This function reduces to spherical Bessel
2987 functions in the limit \eta \to 0.
2990 File: gsl-ref.info, Node: Coulomb Wave Function Normalization Constant, Prev: Coulomb Wave Functions, Up: Coulomb Functions
2992 7.7.3 Coulomb Wave Function Normalization Constant
2993 --------------------------------------------------
2995 The Coulomb wave function normalization constant is defined in
2998 -- Function: int gsl_sf_coulomb_CL_e (double L, double ETA,
2999 gsl_sf_result * RESULT)
3000 This function computes the Coulomb wave function normalization
3001 constant C_L(\eta) for L > -1.
3003 -- Function: int gsl_sf_coulomb_CL_array (double LMIN, int KMAX,
3004 double ETA, double CL[])
3005 This function computes the Coulomb wave function normalization
3006 constant C_L(\eta) for L = Lmin \dots Lmin + kmax, Lmin > -1.
3009 File: gsl-ref.info, Node: Coupling Coefficients, Next: Dawson Function, Prev: Coulomb Functions, Up: Special Functions
3011 7.8 Coupling Coefficients
3012 =========================
3014 The Wigner 3-j, 6-j and 9-j symbols give the coupling coefficients for
3015 combined angular momentum vectors. Since the arguments of the standard
3016 coupling coefficient functions are integer or half-integer, the
3017 arguments of the following functions are, by convention, integers equal
3018 to twice the actual spin value. For information on the 3-j coefficients
3019 see Abramowitz & Stegun, Section 27.9. The functions described in this
3020 section are declared in the header file `gsl_sf_coupling.h'.
3029 File: gsl-ref.info, Node: 3-j Symbols, Next: 6-j Symbols, Up: Coupling Coefficients
3034 -- Function: double gsl_sf_coupling_3j (int TWO_JA, int TWO_JB, int
3035 TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC)
3036 -- Function: int gsl_sf_coupling_3j_e (int TWO_JA, int TWO_JB, int
3037 TWO_JC, int TWO_MA, int TWO_MB, int TWO_MC, gsl_sf_result *
3039 These routines compute the Wigner 3-j coefficient,
3044 where the arguments are given in half-integer units, ja =
3045 TWO_JA/2, ma = TWO_MA/2, etc.
3048 File: gsl-ref.info, Node: 6-j Symbols, Next: 9-j Symbols, Prev: 3-j Symbols, Up: Coupling Coefficients
3053 -- Function: double gsl_sf_coupling_6j (int TWO_JA, int TWO_JB, int
3054 TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF)
3055 -- Function: int gsl_sf_coupling_6j_e (int TWO_JA, int TWO_JB, int
3056 TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, gsl_sf_result *
3058 These routines compute the Wigner 6-j coefficient,
3063 where the arguments are given in half-integer units, ja =
3064 TWO_JA/2, ma = TWO_MA/2, etc.
3067 File: gsl-ref.info, Node: 9-j Symbols, Prev: 6-j Symbols, Up: Coupling Coefficients
3072 -- Function: double gsl_sf_coupling_9j (int TWO_JA, int TWO_JB, int
3073 TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int
3075 -- Function: int gsl_sf_coupling_9j_e (int TWO_JA, int TWO_JB, int
3076 TWO_JC, int TWO_JD, int TWO_JE, int TWO_JF, int TWO_JG, int
3077 TWO_JH, int TWO_JI, gsl_sf_result * RESULT)
3078 These routines compute the Wigner 9-j coefficient,
3084 where the arguments are given in half-integer units, ja =
3085 TWO_JA/2, ma = TWO_MA/2, etc.
3088 File: gsl-ref.info, Node: Dawson Function, Next: Debye Functions, Prev: Coupling Coefficients, Up: Special Functions
3093 The Dawson integral is defined by \exp(-x^2) \int_0^x dt \exp(t^2). A
3094 table of Dawson's integral can be found in Abramowitz & Stegun, Table
3095 7.5. The Dawson functions are declared in the header file
3098 -- Function: double gsl_sf_dawson (double X)
3099 -- Function: int gsl_sf_dawson_e (double X, gsl_sf_result * RESULT)
3100 These routines compute the value of Dawson's integral for X.
3103 File: gsl-ref.info, Node: Debye Functions, Next: Dilogarithm, Prev: Dawson Function, Up: Special Functions
3105 7.10 Debye Functions
3106 ====================
3108 The Debye functions D_n(x) are defined by the following integral,
3110 D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1))
3112 For further information see Abramowitz & Stegun, Section 27.1. The
3113 Debye functions are declared in the header file `gsl_sf_debye.h'.
3115 -- Function: double gsl_sf_debye_1 (double X)
3116 -- Function: int gsl_sf_debye_1_e (double X, gsl_sf_result * RESULT)
3117 These routines compute the first-order Debye function D_1(x) =
3118 (1/x) \int_0^x dt (t/(e^t - 1)).
3120 -- Function: double gsl_sf_debye_2 (double X)
3121 -- Function: int gsl_sf_debye_2_e (double X, gsl_sf_result * RESULT)
3122 These routines compute the second-order Debye function D_2(x) =
3123 (2/x^2) \int_0^x dt (t^2/(e^t - 1)).
3125 -- Function: double gsl_sf_debye_3 (double X)
3126 -- Function: int gsl_sf_debye_3_e (double X, gsl_sf_result * RESULT)
3127 These routines compute the third-order Debye function D_3(x) =
3128 (3/x^3) \int_0^x dt (t^3/(e^t - 1)).
3130 -- Function: double gsl_sf_debye_4 (double X)
3131 -- Function: int gsl_sf_debye_4_e (double X, gsl_sf_result * RESULT)
3132 These routines compute the fourth-order Debye function D_4(x) =
3133 (4/x^4) \int_0^x dt (t^4/(e^t - 1)).
3135 -- Function: double gsl_sf_debye_5 (double X)
3136 -- Function: int gsl_sf_debye_5_e (double X, gsl_sf_result * RESULT)
3137 These routines compute the fifth-order Debye function D_5(x) =
3138 (5/x^5) \int_0^x dt (t^5/(e^t - 1)).
3140 -- Function: double gsl_sf_debye_6 (double X)
3141 -- Function: int gsl_sf_debye_6_e (double X, gsl_sf_result * RESULT)
3142 These routines compute the sixth-order Debye function D_6(x) =
3143 (6/x^6) \int_0^x dt (t^6/(e^t - 1)).
3146 File: gsl-ref.info, Node: Dilogarithm, Next: Elementary Operations, Prev: Debye Functions, Up: Special Functions
3151 The functions described in this section are declared in the header file
3157 * Complex Argument::
3160 File: gsl-ref.info, Node: Real Argument, Next: Complex Argument, Up: Dilogarithm
3162 7.11.1 Real Argument
3163 --------------------
3165 -- Function: double gsl_sf_dilog (double X)
3166 -- Function: int gsl_sf_dilog_e (double X, gsl_sf_result * RESULT)
3167 These routines compute the dilogarithm for a real argument. In
3168 Lewin's notation this is Li_2(x), the real part of the dilogarithm
3169 of a real x. It is defined by the integral representation Li_2(x)
3170 = - \Re \int_0^x ds \log(1-s) / s. Note that \Im(Li_2(x)) = 0 for
3171 x <= 1, and -\pi\log(x) for x > 1.
3173 Note that Abramowitz & Stegun refer to the Spence integral
3174 S(x)=Li_2(1-x) as the dilogarithm rather than Li_2(x).
3177 File: gsl-ref.info, Node: Complex Argument, Prev: Real Argument, Up: Dilogarithm
3179 7.11.2 Complex Argument
3180 -----------------------
3182 -- Function: int gsl_sf_complex_dilog_e (double R, double THETA,
3183 gsl_sf_result * RESULT_RE, gsl_sf_result * RESULT_IM)
3184 This function computes the full complex-valued dilogarithm for the
3185 complex argument z = r \exp(i \theta). The real and imaginary
3186 parts of the result are returned in RESULT_RE, RESULT_IM.
3189 File: gsl-ref.info, Node: Elementary Operations, Next: Elliptic Integrals, Prev: Dilogarithm, Up: Special Functions
3191 7.12 Elementary Operations
3192 ==========================
3194 The following functions allow for the propagation of errors when
3195 combining quantities by multiplication. The functions are declared in
3196 the header file `gsl_sf_elementary.h'.
3198 -- Function: int gsl_sf_multiply_e (double X, double Y, gsl_sf_result
3200 This function multiplies X and Y storing the product and its
3201 associated error in RESULT.
3203 -- Function: int gsl_sf_multiply_err_e (double X, double DX, double Y,
3204 double DY, gsl_sf_result * RESULT)
3205 This function multiplies X and Y with associated absolute errors
3206 DX and DY. The product xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2) is
3210 File: gsl-ref.info, Node: Elliptic Integrals, Next: Elliptic Functions (Jacobi), Prev: Elementary Operations, Up: Special Functions
3212 7.13 Elliptic Integrals
3213 =======================
3215 The functions described in this section are declared in the header file
3216 `gsl_sf_ellint.h'. Further information about the elliptic integrals
3217 can be found in Abramowitz & Stegun, Chapter 17.
3221 * Definition of Legendre Forms::
3222 * Definition of Carlson Forms::
3223 * Legendre Form of Complete Elliptic Integrals::
3224 * Legendre Form of Incomplete Elliptic Integrals::
3228 File: gsl-ref.info, Node: Definition of Legendre Forms, Next: Definition of Carlson Forms, Up: Elliptic Integrals
3230 7.13.1 Definition of Legendre Forms
3231 -----------------------------------
3233 The Legendre forms of elliptic integrals F(\phi,k), E(\phi,k) and
3234 \Pi(\phi,k,n) are defined by,
3236 F(\phi,k) = \int_0^\phi dt 1/\sqrt((1 - k^2 \sin^2(t)))
3238 E(\phi,k) = \int_0^\phi dt \sqrt((1 - k^2 \sin^2(t)))
3240 Pi(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))
3242 The complete Legendre forms are denoted by K(k) = F(\pi/2, k) and E(k)
3245 The notation used here is based on Carlson, `Numerische Mathematik'
3246 33 (1979) 1 and differs slightly from that used by Abramowitz & Stegun,
3247 where the functions are given in terms of the parameter m = k^2 and n
3251 File: gsl-ref.info, Node: Definition of Carlson Forms, Next: Legendre Form of Complete Elliptic Integrals, Prev: Definition of Legendre Forms, Up: Elliptic Integrals
3253 7.13.2 Definition of Carlson Forms
3254 ----------------------------------
3256 The Carlson symmetric forms of elliptical integrals RC(x,y), RD(x,y,z),
3257 RF(x,y,z) and RJ(x,y,z,p) are defined by,
3259 RC(x,y) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1)
3261 RD(x,y,z) = 3/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-3/2)
3263 RF(x,y,z) = 1/2 \int_0^\infty dt (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2)
3265 RJ(x,y,z,p) = 3/2 \int_0^\infty dt
3266 (t+x)^(-1/2) (t+y)^(-1/2) (t+z)^(-1/2) (t+p)^(-1)
3269 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
3271 7.13.3 Legendre Form of Complete Elliptic Integrals
3272 ---------------------------------------------------
3274 -- Function: double gsl_sf_ellint_Kcomp (double K, gsl_mode_t MODE)
3275 -- Function: int gsl_sf_ellint_Kcomp_e (double K, gsl_mode_t MODE,
3276 gsl_sf_result * RESULT)
3277 These routines compute the complete elliptic integral K(k) to the
3278 accuracy specified by the mode variable MODE. Note that
3279 Abramowitz & Stegun define this function in terms of the parameter
3282 -- Function: double gsl_sf_ellint_Ecomp (double K, gsl_mode_t MODE)
3283 -- Function: int gsl_sf_ellint_Ecomp_e (double K, gsl_mode_t MODE,
3284 gsl_sf_result * RESULT)
3285 These routines compute the complete elliptic integral E(k) to the
3286 accuracy specified by the mode variable MODE. Note that
3287 Abramowitz & Stegun define this function in terms of the parameter
3290 -- Function: double gsl_sf_ellint_Pcomp (double K, double N,
3292 -- Function: int gsl_sf_ellint_Pcomp_e (double K, double N, gsl_mode_t
3293 MODE, gsl_sf_result * RESULT)
3294 These routines compute the complete elliptic integral \Pi(k,n) to
3295 the accuracy specified by the mode variable MODE. Note that
3296 Abramowitz & Stegun define this function in terms of the
3297 parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of
3301 File: gsl-ref.info, Node: Legendre Form of Incomplete Elliptic Integrals, Next: Carlson Forms, Prev: Legendre Form of Complete Elliptic Integrals, Up: Elliptic Integrals
3303 7.13.4 Legendre Form of Incomplete Elliptic Integrals
3304 -----------------------------------------------------
3306 -- Function: double gsl_sf_ellint_F (double PHI, double K, gsl_mode_t
3308 -- Function: int gsl_sf_ellint_F_e (double PHI, double K, gsl_mode_t
3309 MODE, gsl_sf_result * RESULT)
3310 These routines compute the incomplete elliptic integral F(\phi,k)
3311 to the accuracy specified by the mode variable MODE. Note that
3312 Abramowitz & Stegun define this function in terms of the parameter
3315 -- Function: double gsl_sf_ellint_E (double PHI, double K, gsl_mode_t
3317 -- Function: int gsl_sf_ellint_E_e (double PHI, double K, gsl_mode_t
3318 MODE, gsl_sf_result * RESULT)
3319 These routines compute the incomplete elliptic integral E(\phi,k)
3320 to the accuracy specified by the mode variable MODE. Note that
3321 Abramowitz & Stegun define this function in terms of the parameter
3324 -- Function: double gsl_sf_ellint_P (double PHI, double K, double N,
3326 -- Function: int gsl_sf_ellint_P_e (double PHI, double K, double N,
3327 gsl_mode_t MODE, gsl_sf_result * RESULT)
3328 These routines compute the incomplete elliptic integral
3329 \Pi(\phi,k,n) to the accuracy specified by the mode variable MODE.
3330 Note that Abramowitz & Stegun define this function in terms of the
3331 parameters m = k^2 and \sin^2(\alpha) = k^2, with the change of
3334 -- Function: double gsl_sf_ellint_D (double PHI, double K, double N,
3336 -- Function: int gsl_sf_ellint_D_e (double PHI, double K, double N,
3337 gsl_mode_t MODE, gsl_sf_result * RESULT)
3338 These functions compute the incomplete elliptic integral D(\phi,k)
3339 which is defined through the Carlson form RD(x,y,z) by the
3342 D(\phi,k,n) = (1/3)(\sin(\phi))^3 RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).
3343 The argument N is not used and will be removed in a future release.
3347 File: gsl-ref.info, Node: Carlson Forms, Prev: Legendre Form of Incomplete Elliptic Integrals, Up: Elliptic Integrals
3349 7.13.5 Carlson Forms
3350 --------------------
3352 -- Function: double gsl_sf_ellint_RC (double X, double Y, gsl_mode_t
3354 -- Function: int gsl_sf_ellint_RC_e (double X, double Y, gsl_mode_t
3355 MODE, gsl_sf_result * RESULT)
3356 These routines compute the incomplete elliptic integral RC(x,y) to
3357 the accuracy specified by the mode variable MODE.
3359 -- Function: double gsl_sf_ellint_RD (double X, double Y, double Z,
3361 -- Function: int gsl_sf_ellint_RD_e (double X, double Y, double Z,
3362 gsl_mode_t MODE, gsl_sf_result * RESULT)
3363 These routines compute the incomplete elliptic integral RD(x,y,z)
3364 to the accuracy specified by the mode variable MODE.
3366 -- Function: double gsl_sf_ellint_RF (double X, double Y, double Z,
3368 -- Function: int gsl_sf_ellint_RF_e (double X, double Y, double Z,
3369 gsl_mode_t MODE, gsl_sf_result * RESULT)
3370 These routines compute the incomplete elliptic integral RF(x,y,z)
3371 to the accuracy specified by the mode variable MODE.
3373 -- Function: double gsl_sf_ellint_RJ (double X, double Y, double Z,
3374 double P, gsl_mode_t MODE)
3375 -- Function: int gsl_sf_ellint_RJ_e (double X, double Y, double Z,
3376 double P, gsl_mode_t MODE, gsl_sf_result * RESULT)
3377 These routines compute the incomplete elliptic integral RJ(x,y,z,p)
3378 to the accuracy specified by the mode variable MODE.
3381 File: gsl-ref.info, Node: Elliptic Functions (Jacobi), Next: Error Functions, Prev: Elliptic Integrals, Up: Special Functions
3383 7.14 Elliptic Functions (Jacobi)
3384 ================================
3386 The Jacobian Elliptic functions are defined in Abramowitz & Stegun,
3387 Chapter 16. The functions are declared in the header file
3390 -- Function: int gsl_sf_elljac_e (double U, double M, double * SN,
3391 double * CN, double * DN)
3392 This function computes the Jacobian elliptic functions sn(u|m),
3393 cn(u|m), dn(u|m) by descending Landen transformations.
3396 File: gsl-ref.info, Node: Error Functions, Next: Exponential Functions, Prev: Elliptic Functions (Jacobi), Up: Special Functions
3398 7.15 Error Functions
3399 ====================
3401 The error function is described in Abramowitz & Stegun, Chapter 7. The
3402 functions in this section are declared in the header file
3408 * Complementary Error Function::
3409 * Log Complementary Error Function::
3410 * Probability functions::
3413 File: gsl-ref.info, Node: Error Function, Next: Complementary Error Function, Up: Error Functions
3415 7.15.1 Error Function
3416 ---------------------
3418 -- Function: double gsl_sf_erf (double X)
3419 -- Function: int gsl_sf_erf_e (double X, gsl_sf_result * RESULT)
3420 These routines compute the error function erf(x), where erf(x) =
3421 (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2).
3424 File: gsl-ref.info, Node: Complementary Error Function, Next: Log Complementary Error Function, Prev: Error Function, Up: Error Functions
3426 7.15.2 Complementary Error Function
3427 -----------------------------------
3429 -- Function: double gsl_sf_erfc (double X)
3430 -- Function: int gsl_sf_erfc_e (double X, gsl_sf_result * RESULT)
3431 These routines compute the complementary error function erfc(x) =
3432 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2).
3435 File: gsl-ref.info, Node: Log Complementary Error Function, Next: Probability functions, Prev: Complementary Error Function, Up: Error Functions
3437 7.15.3 Log Complementary Error Function
3438 ---------------------------------------
3440 -- Function: double gsl_sf_log_erfc (double X)
3441 -- Function: int gsl_sf_log_erfc_e (double X, gsl_sf_result * RESULT)
3442 These routines compute the logarithm of the complementary error
3443 function \log(\erfc(x)).
3446 File: gsl-ref.info, Node: Probability functions, Prev: Log Complementary Error Function, Up: Error Functions
3448 7.15.4 Probability functions
3449 ----------------------------
3451 The probability functions for the Normal or Gaussian distribution are
3452 described in Abramowitz & Stegun, Section 26.2.
3454 -- Function: double gsl_sf_erf_Z (double X)
3455 -- Function: int gsl_sf_erf_Z_e (double X, gsl_sf_result * RESULT)
3456 These routines compute the Gaussian probability density function
3457 Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2).
3459 -- Function: double gsl_sf_erf_Q (double X)
3460 -- Function: int gsl_sf_erf_Q_e (double X, gsl_sf_result * RESULT)
3461 These routines compute the upper tail of the Gaussian probability
3462 function Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2).
3464 The "hazard function" for the normal distribution, also known as the
3465 inverse Mill's ratio, is defined as,
3467 h(x) = Z(x)/Q(x) = \sqrt{2/\pi} \exp(-x^2 / 2) / \erfc(x/\sqrt 2)
3469 It decreases rapidly as x approaches -\infty and asymptotes to h(x)
3470 \sim x as x approaches +\infty.
3472 -- Function: double gsl_sf_hazard (double X)
3473 -- Function: int gsl_sf_hazard_e (double X, gsl_sf_result * RESULT)
3474 These routines compute the hazard function for the normal
3478 File: gsl-ref.info, Node: Exponential Functions, Next: Exponential Integrals, Prev: Error Functions, Up: Special Functions
3480 7.16 Exponential Functions
3481 ==========================
3483 The functions described in this section are declared in the header file
3488 * Exponential Function::
3489 * Relative Exponential Functions::
3490 * Exponentiation With Error Estimate::
3493 File: gsl-ref.info, Node: Exponential Function, Next: Relative Exponential Functions, Up: Exponential Functions
3495 7.16.1 Exponential Function
3496 ---------------------------
3498 -- Function: double gsl_sf_exp (double X)
3499 -- Function: int gsl_sf_exp_e (double X, gsl_sf_result * RESULT)
3500 These routines provide an exponential function \exp(x) using GSL
3501 semantics and error checking.
3503 -- Function: int gsl_sf_exp_e10_e (double X, gsl_sf_result_e10 *
3505 This function computes the exponential \exp(x) using the
3506 `gsl_sf_result_e10' type to return a result with extended range.
3507 This function may be useful if the value of \exp(x) would overflow
3508 the numeric range of `double'.
3510 -- Function: double gsl_sf_exp_mult (double X, double Y)
3511 -- Function: int gsl_sf_exp_mult_e (double X, double Y, gsl_sf_result
3513 These routines exponentiate X and multiply by the factor Y to
3514 return the product y \exp(x).
3516 -- Function: int gsl_sf_exp_mult_e10_e (const double X, const double
3517 Y, gsl_sf_result_e10 * RESULT)
3518 This function computes the product y \exp(x) using the
3519 `gsl_sf_result_e10' type to return a result with extended numeric
3523 File: gsl-ref.info, Node: Relative Exponential Functions, Next: Exponentiation With Error Estimate, Prev: Exponential Function, Up: Exponential Functions
3525 7.16.2 Relative Exponential Functions
3526 -------------------------------------
3528 -- Function: double gsl_sf_expm1 (double X)
3529 -- Function: int gsl_sf_expm1_e (double X, gsl_sf_result * RESULT)
3530 These routines compute the quantity \exp(x)-1 using an algorithm
3531 that is accurate for small x.
3533 -- Function: double gsl_sf_exprel (double X)
3534 -- Function: int gsl_sf_exprel_e (double X, gsl_sf_result * RESULT)
3535 These routines compute the quantity (\exp(x)-1)/x using an
3536 algorithm that is accurate for small x. For small x the algorithm
3537 is based on the expansion (\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) +
3538 x^3/(2*3*4) + \dots.
3540 -- Function: double gsl_sf_exprel_2 (double X)
3541 -- Function: int gsl_sf_exprel_2_e (double X, gsl_sf_result * RESULT)
3542 These routines compute the quantity 2(\exp(x)-1-x)/x^2 using an
3543 algorithm that is accurate for small x. For small x the algorithm
3544 is based on the expansion 2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4)
3545 + x^3/(3*4*5) + \dots.
3547 -- Function: double gsl_sf_exprel_n (int N, double X)
3548 -- Function: int gsl_sf_exprel_n_e (int N, double X, gsl_sf_result *
3550 These routines compute the N-relative exponential, which is the
3551 N-th generalization of the functions `gsl_sf_exprel' and
3552 `gsl_sf_exprel2'. The N-relative exponential is given by,
3554 exprel_N(x) = N!/x^N (\exp(x) - \sum_{k=0}^{N-1} x^k/k!)
3555 = 1 + x/(N+1) + x^2/((N+1)(N+2)) + ...
3559 File: gsl-ref.info, Node: Exponentiation With Error Estimate, Prev: Relative Exponential Functions, Up: Exponential Functions
3561 7.16.3 Exponentiation With Error Estimate
3562 -----------------------------------------
3564 -- Function: int gsl_sf_exp_err_e (double X, double DX, gsl_sf_result
3566 This function exponentiates X with an associated absolute error DX.
3568 -- Function: int gsl_sf_exp_err_e10_e (double X, double DX,
3569 gsl_sf_result_e10 * RESULT)
3570 This function exponentiates a quantity X with an associated
3571 absolute error DX using the `gsl_sf_result_e10' type to return a
3572 result with extended range.
3574 -- Function: int gsl_sf_exp_mult_err_e (double X, double DX, double Y,
3575 double DY, gsl_sf_result * RESULT)
3576 This routine computes the product y \exp(x) for the quantities X,
3577 Y with associated absolute errors DX, DY.
3579 -- Function: int gsl_sf_exp_mult_err_e10_e (double X, double DX,
3580 double Y, double DY, gsl_sf_result_e10 * RESULT)
3581 This routine computes the product y \exp(x) for the quantities X,
3582 Y with associated absolute errors DX, DY using the
3583 `gsl_sf_result_e10' type to return a result with extended range.
3586 File: gsl-ref.info, Node: Exponential Integrals, Next: Fermi-Dirac Function, Prev: Exponential Functions, Up: Special Functions
3588 7.17 Exponential Integrals
3589 ==========================
3591 Information on the exponential integrals can be found in Abramowitz &
3592 Stegun, Chapter 5. These functions are declared in the header file
3597 * Exponential Integral::
3599 * Hyperbolic Integrals::
3601 * Trigonometric Integrals::
3602 * Arctangent Integral::
3605 File: gsl-ref.info, Node: Exponential Integral, Next: Ei(x), Up: Exponential Integrals
3607 7.17.1 Exponential Integral
3608 ---------------------------
3610 -- Function: double gsl_sf_expint_E1 (double X)
3611 -- Function: int gsl_sf_expint_E1_e (double X, gsl_sf_result * RESULT)
3612 These routines compute the exponential integral E_1(x),
3614 E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t.
3618 -- Function: double gsl_sf_expint_E2 (double X)
3619 -- Function: int gsl_sf_expint_E2_e (double X, gsl_sf_result * RESULT)
3620 These routines compute the second-order exponential integral
3623 E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.
3627 -- Function: double gsl_sf_expint_En (int N, double X)
3628 -- Function: int gsl_sf_expint_En_e (int N, double X, gsl_sf_result *
3630 These routines compute the exponential integral E_n(x) of order n,
3632 E_n(x) := \Re \int_1^\infty dt \exp(-xt)/t^n.
3637 File: gsl-ref.info, Node: Ei(x), Next: Hyperbolic Integrals, Prev: Exponential Integral, Up: Exponential Integrals
3642 -- Function: double gsl_sf_expint_Ei (double X)
3643 -- Function: int gsl_sf_expint_Ei_e (double X, gsl_sf_result * RESULT)
3644 These routines compute the exponential integral Ei(x),
3646 Ei(x) := - PV(\int_{-x}^\infty dt \exp(-t)/t)
3648 where PV denotes the principal value of the integral.
3651 File: gsl-ref.info, Node: Hyperbolic Integrals, Next: Ei_3(x), Prev: Ei(x), Up: Exponential Integrals
3653 7.17.3 Hyperbolic Integrals
3654 ---------------------------
3656 -- Function: double gsl_sf_Shi (double X)
3657 -- Function: int gsl_sf_Shi_e (double X, gsl_sf_result * RESULT)
3658 These routines compute the integral Shi(x) = \int_0^x dt
3661 -- Function: double gsl_sf_Chi (double X)
3662 -- Function: int gsl_sf_Chi_e (double X, gsl_sf_result * RESULT)
3663 These routines compute the integral Chi(x) := \Re[ \gamma_E +
3664 \log(x) + \int_0^x dt (\cosh[t]-1)/t] , where \gamma_E is the
3665 Euler constant (available as the macro `M_EULER').
3668 File: gsl-ref.info, Node: Ei_3(x), Next: Trigonometric Integrals, Prev: Hyperbolic Integrals, Up: Exponential Integrals
3673 -- Function: double gsl_sf_expint_3 (double X)
3674 -- Function: int gsl_sf_expint_3_e (double X, gsl_sf_result * RESULT)
3675 These routines compute the third-order exponential integral
3676 Ei_3(x) = \int_0^xdt \exp(-t^3) for x >= 0.
3679 File: gsl-ref.info, Node: Trigonometric Integrals, Next: Arctangent Integral, Prev: Ei_3(x), Up: Exponential Integrals
3681 7.17.5 Trigonometric Integrals
3682 ------------------------------
3684 -- Function: double gsl_sf_Si (const double X)
3685 -- Function: int gsl_sf_Si_e (double X, gsl_sf_result * RESULT)
3686 These routines compute the Sine integral Si(x) = \int_0^x dt
3689 -- Function: double gsl_sf_Ci (const double X)
3690 -- Function: int gsl_sf_Ci_e (double X, gsl_sf_result * RESULT)
3691 These routines compute the Cosine integral Ci(x) = -\int_x^\infty
3692 dt \cos(t)/t for x > 0.
3695 File: gsl-ref.info, Node: Arctangent Integral, Prev: Trigonometric Integrals, Up: Exponential Integrals
3697 7.17.6 Arctangent Integral
3698 --------------------------
3700 -- Function: double gsl_sf_atanint (double X)
3701 -- Function: int gsl_sf_atanint_e (double X, gsl_sf_result * RESULT)
3702 These routines compute the Arctangent integral, which is defined
3703 as AtanInt(x) = \int_0^x dt \arctan(t)/t.
3706 File: gsl-ref.info, Node: Fermi-Dirac Function, Next: Gamma and Beta Functions, Prev: Exponential Integrals, Up: Special Functions
3708 7.18 Fermi-Dirac Function
3709 =========================
3711 The functions described in this section are declared in the header file
3712 `gsl_sf_fermi_dirac.h'.
3716 * Complete Fermi-Dirac Integrals::
3717 * Incomplete Fermi-Dirac Integrals::
3720 File: gsl-ref.info, Node: Complete Fermi-Dirac Integrals, Next: Incomplete Fermi-Dirac Integrals, Up: Fermi-Dirac Function
3722 7.18.1 Complete Fermi-Dirac Integrals
3723 -------------------------------------
3725 The complete Fermi-Dirac integral F_j(x) is given by,
3727 F_j(x) := (1/\Gamma(j+1)) \int_0^\infty dt (t^j / (\exp(t-x) + 1))
3728 Note that the Fermi-Dirac integral is sometimes defined without the
3729 normalisation factor in other texts.
3731 -- Function: double gsl_sf_fermi_dirac_m1 (double X)
3732 -- Function: int gsl_sf_fermi_dirac_m1_e (double X, gsl_sf_result *
3734 These routines compute the complete Fermi-Dirac integral with an
3735 index of -1. This integral is given by F_{-1}(x) = e^x / (1 +
3738 -- Function: double gsl_sf_fermi_dirac_0 (double X)
3739 -- Function: int gsl_sf_fermi_dirac_0_e (double X, gsl_sf_result *
3741 These routines compute the complete Fermi-Dirac integral with an
3742 index of 0. This integral is given by F_0(x) = \ln(1 + e^x).
3744 -- Function: double gsl_sf_fermi_dirac_1 (double X)
3745 -- Function: int gsl_sf_fermi_dirac_1_e (double X, gsl_sf_result *
3747 These routines compute the complete Fermi-Dirac integral with an
3748 index of 1, F_1(x) = \int_0^\infty dt (t /(\exp(t-x)+1)).
3750 -- Function: double gsl_sf_fermi_dirac_2 (double X)
3751 -- Function: int gsl_sf_fermi_dirac_2_e (double X, gsl_sf_result *
3753 These routines compute the complete Fermi-Dirac integral with an
3754 index of 2, F_2(x) = (1/2) \int_0^\infty dt (t^2 /(\exp(t-x)+1)).
3756 -- Function: double gsl_sf_fermi_dirac_int (int J, double X)
3757 -- Function: int gsl_sf_fermi_dirac_int_e (int J, double X,
3758 gsl_sf_result * RESULT)
3759 These routines compute the complete Fermi-Dirac integral with an
3760 integer index of j, F_j(x) = (1/\Gamma(j+1)) \int_0^\infty dt (t^j
3763 -- Function: double gsl_sf_fermi_dirac_mhalf (double X)
3764 -- Function: int gsl_sf_fermi_dirac_mhalf_e (double X, gsl_sf_result *
3766 These routines compute the complete Fermi-Dirac integral
3769 -- Function: double gsl_sf_fermi_dirac_half (double X)
3770 -- Function: int gsl_sf_fermi_dirac_half_e (double X, gsl_sf_result *
3772 These routines compute the complete Fermi-Dirac integral
3775 -- Function: double gsl_sf_fermi_dirac_3half (double X)
3776 -- Function: int gsl_sf_fermi_dirac_3half_e (double X, gsl_sf_result *
3778 These routines compute the complete Fermi-Dirac integral
3782 File: gsl-ref.info, Node: Incomplete Fermi-Dirac Integrals, Prev: Complete Fermi-Dirac Integrals, Up: Fermi-Dirac Function
3784 7.18.2 Incomplete Fermi-Dirac Integrals
3785 ---------------------------------------
3787 The incomplete Fermi-Dirac integral F_j(x,b) is given by,
3789 F_j(x,b) := (1/\Gamma(j+1)) \int_b^\infty dt (t^j / (\Exp(t-x) + 1))
3791 -- Function: double gsl_sf_fermi_dirac_inc_0 (double X, double B)
3792 -- Function: int gsl_sf_fermi_dirac_inc_0_e (double X, double B,
3793 gsl_sf_result * RESULT)
3794 These routines compute the incomplete Fermi-Dirac integral with an
3795 index of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x).
3798 File: gsl-ref.info, Node: Gamma and Beta Functions, Next: Gegenbauer Functions, Prev: Fermi-Dirac Function, Up: Special Functions
3800 7.19 Gamma and Beta Functions
3801 =============================
3803 The functions described in this section are declared in the header
3804 file `gsl_sf_gamma.h'.
3810 * Pochhammer Symbol::
3811 * Incomplete Gamma Functions::
3813 * Incomplete Beta Function::
3816 File: gsl-ref.info, Node: Gamma Functions, Next: Factorials, Up: Gamma and Beta Functions
3818 7.19.1 Gamma Functions
3819 ----------------------
3821 The Gamma function is defined by the following integral,
3823 \Gamma(x) = \int_0^\infty dt t^{x-1} \exp(-t)
3825 It is related to the factorial function by \Gamma(n)=(n-1)! for
3826 positive integer n. Further information on the Gamma function can be
3827 found in Abramowitz & Stegun, Chapter 6. The functions described in
3828 this section are declared in the header file `gsl_sf_gamma.h'.
3830 -- Function: double gsl_sf_gamma (double X)
3831 -- Function: int gsl_sf_gamma_e (double X, gsl_sf_result * RESULT)
3832 These routines compute the Gamma function \Gamma(x), subject to x
3833 not being a negative integer or zero. The function is computed
3834 using the real Lanczos method. The maximum value of x such that
3835 \Gamma(x) is not considered an overflow is given by the macro
3836 `GSL_SF_GAMMA_XMAX' and is 171.0.
3838 -- Function: double gsl_sf_lngamma (double X)
3839 -- Function: int gsl_sf_lngamma_e (double X, gsl_sf_result * RESULT)
3840 These routines compute the logarithm of the Gamma function,
3841 \log(\Gamma(x)), subject to x not being a negative integer or
3842 zero. For x<0 the real part of \log(\Gamma(x)) is returned, which
3843 is equivalent to \log(|\Gamma(x)|). The function is computed
3844 using the real Lanczos method.
3846 -- Function: int gsl_sf_lngamma_sgn_e (double X, gsl_sf_result *
3847 RESULT_LG, double * SGN)
3848 This routine computes the sign of the gamma function and the
3849 logarithm of its magnitude, subject to x not being a negative
3850 integer or zero. The function is computed using the real Lanczos
3851 method. The value of the gamma function can be reconstructed
3852 using the relation \Gamma(x) = sgn * \exp(resultlg).
3854 -- Function: double gsl_sf_gammastar (double X)
3855 -- Function: int gsl_sf_gammastar_e (double X, gsl_sf_result * RESULT)
3856 These routines compute the regulated Gamma Function \Gamma^*(x)
3857 for x > 0. The regulated gamma function is given by,
3859 \Gamma^*(x) = \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x))
3860 = (1 + (1/12x) + ...) for x \to \infty
3861 and is a useful suggestion of Temme.
3863 -- Function: double gsl_sf_gammainv (double X)
3864 -- Function: int gsl_sf_gammainv_e (double X, gsl_sf_result * RESULT)
3865 These routines compute the reciprocal of the gamma function,
3866 1/\Gamma(x) using the real Lanczos method.
3868 -- Function: int gsl_sf_lngamma_complex_e (double ZR, double ZI,
3869 gsl_sf_result * LNR, gsl_sf_result * ARG)
3870 This routine computes \log(\Gamma(z)) for complex z=z_r+i z_i and
3871 z not a negative integer or zero, using the complex Lanczos
3872 method. The returned parameters are lnr = \log|\Gamma(z)| and arg
3873 = \arg(\Gamma(z)) in (-\pi,\pi]. Note that the phase part (ARG)
3874 is not well-determined when |z| is very large, due to inevitable
3875 roundoff in restricting to (-\pi,\pi]. This will result in a
3876 `GSL_ELOSS' error when it occurs. The absolute value part (LNR),
3877 however, never suffers from loss of precision.
3880 File: gsl-ref.info, Node: Factorials, Next: Pochhammer Symbol, Prev: Gamma Functions, Up: Gamma and Beta Functions
3885 Although factorials can be computed from the Gamma function, using the
3886 relation n! = \Gamma(n+1) for non-negative integer n, it is usually
3887 more efficient to call the functions in this section, particularly for
3888 small values of n, whose factorial values are maintained in hardcoded
3891 -- Function: double gsl_sf_fact (unsigned int N)
3892 -- Function: int gsl_sf_fact_e (unsigned int N, gsl_sf_result * RESULT)
3893 These routines compute the factorial n!. The factorial is related
3894 to the Gamma function by n! = \Gamma(n+1). The maximum value of n
3895 such that n! is not considered an overflow is given by the macro
3896 `GSL_SF_FACT_NMAX' and is 170.
3898 -- Function: double gsl_sf_doublefact (unsigned int N)
3899 -- Function: int gsl_sf_doublefact_e (unsigned int N, gsl_sf_result *
3901 These routines compute the double factorial n!! = n(n-2)(n-4)
3902 \dots. The maximum value of n such that n!! is not considered an
3903 overflow is given by the macro `GSL_SF_DOUBLEFACT_NMAX' and is 297.
3905 -- Function: double gsl_sf_lnfact (unsigned int N)
3906 -- Function: int gsl_sf_lnfact_e (unsigned int N, gsl_sf_result *
3908 These routines compute the logarithm of the factorial of N,
3909 \log(n!). The algorithm is faster than computing \ln(\Gamma(n+1))
3910 via `gsl_sf_lngamma' for n < 170, but defers for larger N.
3912 -- Function: double gsl_sf_lndoublefact (unsigned int N)
3913 -- Function: int gsl_sf_lndoublefact_e (unsigned int N, gsl_sf_result
3915 These routines compute the logarithm of the double factorial of N,
3918 -- Function: double gsl_sf_choose (unsigned int N, unsigned int M)
3919 -- Function: int gsl_sf_choose_e (unsigned int N, unsigned int M,
3920 gsl_sf_result * RESULT)
3921 These routines compute the combinatorial factor `n choose m' =
3924 -- Function: double gsl_sf_lnchoose (unsigned int N, unsigned int M)
3925 -- Function: int gsl_sf_lnchoose_e (unsigned int N, unsigned int M,
3926 gsl_sf_result * RESULT)
3927 These routines compute the logarithm of `n choose m'. This is
3928 equivalent to the sum \log(n!) - \log(m!) - \log((n-m)!).
3930 -- Function: double gsl_sf_taylorcoeff (int N, double X)
3931 -- Function: int gsl_sf_taylorcoeff_e (int N, double X, gsl_sf_result
3933 These routines compute the Taylor coefficient x^n / n! for x >= 0,
3937 File: gsl-ref.info, Node: Pochhammer Symbol, Next: Incomplete Gamma Functions, Prev: Factorials, Up: Gamma and Beta Functions
3939 7.19.3 Pochhammer Symbol
3940 ------------------------
3942 -- Function: double gsl_sf_poch (double A, double X)
3943 -- Function: int gsl_sf_poch_e (double A, double X, gsl_sf_result *
3945 These routines compute the Pochhammer symbol (a)_x = \Gamma(a +
3946 x)/\Gamma(a), subject to a and a+x not being negative integers or
3947 zero. The Pochhammer symbol is also known as the Apell symbol and
3948 sometimes written as (a,x).
3950 -- Function: double gsl_sf_lnpoch (double A, double X)
3951 -- Function: int gsl_sf_lnpoch_e (double A, double X, gsl_sf_result *
3953 These routines compute the logarithm of the Pochhammer symbol,
3954 \log((a)_x) = \log(\Gamma(a + x)/\Gamma(a)) for a > 0, a+x > 0.
3956 -- Function: int gsl_sf_lnpoch_sgn_e (double A, double X,
3957 gsl_sf_result * RESULT, double * SGN)
3958 These routines compute the sign of the Pochhammer symbol and the
3959 logarithm of its magnitude. The computed parameters are result =
3960 \log(|(a)_x|) and sgn = \sgn((a)_x) where (a)_x = \Gamma(a +
3961 x)/\Gamma(a), subject to a, a+x not being negative integers or
3964 -- Function: double gsl_sf_pochrel (double A, double X)
3965 -- Function: int gsl_sf_pochrel_e (double A, double X, gsl_sf_result *
3967 These routines compute the relative Pochhammer symbol ((a)_x -
3968 1)/x where (a)_x = \Gamma(a + x)/\Gamma(a).
3971 File: gsl-ref.info, Node: Incomplete Gamma Functions, Next: Beta Functions, Prev: Pochhammer Symbol, Up: Gamma and Beta Functions
3973 7.19.4 Incomplete Gamma Functions
3974 ---------------------------------
3976 -- Function: double gsl_sf_gamma_inc (double A, double X)
3977 -- Function: int gsl_sf_gamma_inc_e (double A, double X, gsl_sf_result
3979 These functions compute the unnormalized incomplete Gamma Function
3980 \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x
3983 -- Function: double gsl_sf_gamma_inc_Q (double A, double X)
3984 -- Function: int gsl_sf_gamma_inc_Q_e (double A, double X,
3985 gsl_sf_result * RESULT)
3986 These routines compute the normalized incomplete Gamma Function
3987 Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0,
3990 -- Function: double gsl_sf_gamma_inc_P (double A, double X)
3991 -- Function: int gsl_sf_gamma_inc_P_e (double A, double X,
3992 gsl_sf_result * RESULT)
3993 These routines compute the complementary normalized incomplete
3994 Gamma Function P(a,x) = 1 - Q(a,x) = 1/\Gamma(a) \int_0^x dt
3995 t^{a-1} \exp(-t) for a > 0, x >= 0.
3997 Note that Abramowitz & Stegun call P(a,x) the incomplete gamma
3998 function (section 6.5).
4001 File: gsl-ref.info, Node: Beta Functions, Next: Incomplete Beta Function, Prev: Incomplete Gamma Functions, Up: Gamma and Beta Functions
4003 7.19.5 Beta Functions
4004 ---------------------
4006 -- Function: double gsl_sf_beta (double A, double B)
4007 -- Function: int gsl_sf_beta_e (double A, double B, gsl_sf_result *
4009 These routines compute the Beta Function, B(a,b) =
4010 \Gamma(a)\Gamma(b)/\Gamma(a+b) subject to a and b not being
4013 -- Function: double gsl_sf_lnbeta (double A, double B)
4014 -- Function: int gsl_sf_lnbeta_e (double A, double B, gsl_sf_result *
4016 These routines compute the logarithm of the Beta Function,
4017 \log(B(a,b)) subject to a and b not being negative integers.
4020 File: gsl-ref.info, Node: Incomplete Beta Function, Prev: Beta Functions, Up: Gamma and Beta Functions
4022 7.19.6 Incomplete Beta Function
4023 -------------------------------
4025 -- Function: double gsl_sf_beta_inc (double A, double B, double X)
4026 -- Function: int gsl_sf_beta_inc_e (double A, double B, double X,
4027 gsl_sf_result * RESULT)
4028 These routines compute the normalized incomplete Beta function
4029 I_x(a,b)=B_x(a,b)/B(a,b) where B_x(a,b) = \int_0^x t^{a-1}
4030 (1-t)^{b-1} dt for 0 <= x <= 1. For a > 0, b > 0 the value is
4031 computed using a continued fraction expansion. For all other
4032 values it is computed using the relation I_x(a,b,x) = (1/a) x^a
4033 2F1(a,1-b,a+1,x)/B(a,b).
4036 File: gsl-ref.info, Node: Gegenbauer Functions, Next: Hypergeometric Functions, Prev: Gamma and Beta Functions, Up: Special Functions
4038 7.20 Gegenbauer Functions
4039 =========================
4041 The Gegenbauer polynomials are defined in Abramowitz & Stegun, Chapter
4042 22, where they are known as Ultraspherical polynomials. The functions
4043 described in this section are declared in the header file
4044 `gsl_sf_gegenbauer.h'.
4046 -- Function: double gsl_sf_gegenpoly_1 (double LAMBDA, double X)
4047 -- Function: double gsl_sf_gegenpoly_2 (double LAMBDA, double X)
4048 -- Function: double gsl_sf_gegenpoly_3 (double LAMBDA, double X)
4049 -- Function: int gsl_sf_gegenpoly_1_e (double LAMBDA, double X,
4050 gsl_sf_result * RESULT)
4051 -- Function: int gsl_sf_gegenpoly_2_e (double LAMBDA, double X,
4052 gsl_sf_result * RESULT)
4053 -- Function: int gsl_sf_gegenpoly_3_e (double LAMBDA, double X,
4054 gsl_sf_result * RESULT)
4055 These functions evaluate the Gegenbauer polynomials
4056 C^{(\lambda)}_n(x) using explicit representations for n =1, 2, 3.
4058 -- Function: double gsl_sf_gegenpoly_n (int N, double LAMBDA, double X)
4059 -- Function: int gsl_sf_gegenpoly_n_e (int N, double LAMBDA, double X,
4060 gsl_sf_result * RESULT)
4061 These functions evaluate the Gegenbauer polynomial
4062 C^{(\lambda)}_n(x) for a specific value of N, LAMBDA, X subject to
4063 \lambda > -1/2, n >= 0.
4065 -- Function: int gsl_sf_gegenpoly_array (int NMAX, double LAMBDA,
4066 double X, double RESULT_ARRAY[])
4067 This function computes an array of Gegenbauer polynomials
4068 C^{(\lambda)}_n(x) for n = 0, 1, 2, \dots, nmax, subject to
4069 \lambda > -1/2, nmax >= 0.
4072 File: gsl-ref.info, Node: Hypergeometric Functions, Next: Laguerre Functions, Prev: Gegenbauer Functions, Up: Special Functions
4074 7.21 Hypergeometric Functions
4075 =============================
4077 Hypergeometric functions are described in Abramowitz & Stegun, Chapters
4078 13 and 15. These functions are declared in the header file
4081 -- Function: double gsl_sf_hyperg_0F1 (double C, double X)
4082 -- Function: int gsl_sf_hyperg_0F1_e (double C, double X,
4083 gsl_sf_result * RESULT)
4084 These routines compute the hypergeometric function 0F1(c,x).
4086 -- Function: double gsl_sf_hyperg_1F1_int (int M, int N, double X)
4087 -- Function: int gsl_sf_hyperg_1F1_int_e (int M, int N, double X,
4088 gsl_sf_result * RESULT)
4089 These routines compute the confluent hypergeometric function
4090 1F1(m,n,x) = M(m,n,x) for integer parameters M, N.
4092 -- Function: double gsl_sf_hyperg_1F1 (double A, double B, double X)
4093 -- Function: int gsl_sf_hyperg_1F1_e (double A, double B, double X,
4094 gsl_sf_result * RESULT)
4095 These routines compute the confluent hypergeometric function
4096 1F1(a,b,x) = M(a,b,x) for general parameters A, B.
4098 -- Function: double gsl_sf_hyperg_U_int (int M, int N, double X)
4099 -- Function: int gsl_sf_hyperg_U_int_e (int M, int N, double X,
4100 gsl_sf_result * RESULT)
4101 These routines compute the confluent hypergeometric function
4102 U(m,n,x) for integer parameters M, N.
4104 -- Function: int gsl_sf_hyperg_U_int_e10_e (int M, int N, double X,
4105 gsl_sf_result_e10 * RESULT)
4106 This routine computes the confluent hypergeometric function
4107 U(m,n,x) for integer parameters M, N using the `gsl_sf_result_e10'
4108 type to return a result with extended range.
4110 -- Function: double gsl_sf_hyperg_U (double A, double B, double X)
4111 -- Function: int gsl_sf_hyperg_U_e (double A, double B, double X,
4112 gsl_sf_result * RESULT)
4113 These routines compute the confluent hypergeometric function
4116 -- Function: int gsl_sf_hyperg_U_e10_e (double A, double B, double X,
4117 gsl_sf_result_e10 * RESULT)
4118 This routine computes the confluent hypergeometric function
4119 U(a,b,x) using the `gsl_sf_result_e10' type to return a result
4120 with extended range.
4122 -- Function: double gsl_sf_hyperg_2F1 (double A, double B, double C,
4124 -- Function: int gsl_sf_hyperg_2F1_e (double A, double B, double C,
4125 double X, gsl_sf_result * RESULT)
4126 These routines compute the Gauss hypergeometric function
4127 2F1(a,b,c,x) = F(a,b,c,x) for |x| < 1.
4129 If the arguments (a,b,c,x) are too close to a singularity then the
4130 function can return the error code `GSL_EMAXITER' when the series
4131 approximation converges too slowly. This occurs in the region of
4132 x=1, c - a - b = m for integer m.
4134 -- Function: double gsl_sf_hyperg_2F1_conj (double AR, double AI,
4136 -- Function: int gsl_sf_hyperg_2F1_conj_e (double AR, double AI,
4137 double C, double X, gsl_sf_result * RESULT)
4138 These routines compute the Gauss hypergeometric function 2F1(a_R +
4139 i a_I, a_R - i a_I, c, x) with complex parameters for |x| < 1.
4142 -- Function: double gsl_sf_hyperg_2F1_renorm (double A, double B,
4144 -- Function: int gsl_sf_hyperg_2F1_renorm_e (double A, double B,
4145 double C, double X, gsl_sf_result * RESULT)
4146 These routines compute the renormalized Gauss hypergeometric
4147 function 2F1(a,b,c,x) / \Gamma(c) for |x| < 1.
4149 -- Function: double gsl_sf_hyperg_2F1_conj_renorm (double AR, double
4150 AI, double C, double X)
4151 -- Function: int gsl_sf_hyperg_2F1_conj_renorm_e (double AR, double
4152 AI, double C, double X, gsl_sf_result * RESULT)
4153 These routines compute the renormalized Gauss hypergeometric
4154 function 2F1(a_R + i a_I, a_R - i a_I, c, x) / \Gamma(c) for |x| <
4157 -- Function: double gsl_sf_hyperg_2F0 (double A, double B, double X)
4158 -- Function: int gsl_sf_hyperg_2F0_e (double A, double B, double X,
4159 gsl_sf_result * RESULT)
4160 These routines compute the hypergeometric function 2F0(a,b,x).
4161 The series representation is a divergent hypergeometric series.
4162 However, for x < 0 we have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)
4165 File: gsl-ref.info, Node: Laguerre Functions, Next: Lambert W Functions, Prev: Hypergeometric Functions, Up: Special Functions
4167 7.22 Laguerre Functions
4168 =======================
4170 The generalized Laguerre polynomials are defined in terms of confluent
4171 hypergeometric functions as L^a_n(x) = ((a+1)_n / n!) 1F1(-n,a+1,x),
4172 and are sometimes referred to as the associated Laguerre polynomials.
4173 They are related to the plain Laguerre polynomials L_n(x) by L^0_n(x) =
4174 L_n(x) and L^k_n(x) = (-1)^k (d^k/dx^k) L_(n+k)(x). For more
4175 information see Abramowitz & Stegun, Chapter 22.
4177 The functions described in this section are declared in the header
4178 file `gsl_sf_laguerre.h'.
4180 -- Function: double gsl_sf_laguerre_1 (double A, double X)
4181 -- Function: double gsl_sf_laguerre_2 (double A, double X)
4182 -- Function: double gsl_sf_laguerre_3 (double A, double X)
4183 -- Function: int gsl_sf_laguerre_1_e (double A, double X,
4184 gsl_sf_result * RESULT)
4185 -- Function: int gsl_sf_laguerre_2_e (double A, double X,
4186 gsl_sf_result * RESULT)
4187 -- Function: int gsl_sf_laguerre_3_e (double A, double X,
4188 gsl_sf_result * RESULT)
4189 These routines evaluate the generalized Laguerre polynomials
4190 L^a_1(x), L^a_2(x), L^a_3(x) using explicit representations.
4192 -- Function: double gsl_sf_laguerre_n (const int N, const double A,
4194 -- Function: int gsl_sf_laguerre_n_e (int N, double A, double X,
4195 gsl_sf_result * RESULT)
4196 These routines evaluate the generalized Laguerre polynomials
4197 L^a_n(x) for a > -1, n >= 0.
4201 File: gsl-ref.info, Node: Lambert W Functions, Next: Legendre Functions and Spherical Harmonics, Prev: Laguerre Functions, Up: Special Functions
4203 7.23 Lambert W Functions
4204 ========================
4206 Lambert's W functions, W(x), are defined to be solutions of the
4207 equation W(x) \exp(W(x)) = x. This function has multiple branches for x
4208 < 0; however, it has only two real-valued branches. We define W_0(x) to
4209 be the principal branch, where W > -1 for x < 0, and W_{-1}(x) to be
4210 the other real branch, where W < -1 for x < 0. The Lambert functions
4211 are declared in the header file `gsl_sf_lambert.h'.
4213 -- Function: double gsl_sf_lambert_W0 (double X)
4214 -- Function: int gsl_sf_lambert_W0_e (double X, gsl_sf_result * RESULT)
4215 These compute the principal branch of the Lambert W function,
4218 -- Function: double gsl_sf_lambert_Wm1 (double X)
4219 -- Function: int gsl_sf_lambert_Wm1_e (double X, gsl_sf_result *
4221 These compute the secondary real-valued branch of the Lambert W
4222 function, W_{-1}(x).
4225 File: gsl-ref.info, Node: Legendre Functions and Spherical Harmonics, Next: Logarithm and Related Functions, Prev: Lambert W Functions, Up: Special Functions
4227 7.24 Legendre Functions and Spherical Harmonics
4228 ===============================================
4230 The Legendre Functions and Legendre Polynomials are described in
4231 Abramowitz & Stegun, Chapter 8. These functions are declared in the
4232 header file `gsl_sf_legendre.h'.
4236 * Legendre Polynomials::
4237 * Associated Legendre Polynomials and Spherical Harmonics::
4238 * Conical Functions::
4239 * Radial Functions for Hyperbolic Space::
4242 File: gsl-ref.info, Node: Legendre Polynomials, Next: Associated Legendre Polynomials and Spherical Harmonics, Up: Legendre Functions and Spherical Harmonics
4244 7.24.1 Legendre Polynomials
4245 ---------------------------
4247 -- Function: double gsl_sf_legendre_P1 (double X)
4248 -- Function: double gsl_sf_legendre_P2 (double X)
4249 -- Function: double gsl_sf_legendre_P3 (double X)
4250 -- Function: int gsl_sf_legendre_P1_e (double X, gsl_sf_result *
4252 -- Function: int gsl_sf_legendre_P2_e (double X, gsl_sf_result *
4254 -- Function: int gsl_sf_legendre_P3_e (double X, gsl_sf_result *
4256 These functions evaluate the Legendre polynomials P_l(x) using
4257 explicit representations for l=1, 2, 3.
4259 -- Function: double gsl_sf_legendre_Pl (int L, double X)
4260 -- Function: int gsl_sf_legendre_Pl_e (int L, double X, gsl_sf_result
4262 These functions evaluate the Legendre polynomial P_l(x) for a
4263 specific value of L, X subject to l >= 0, |x| <= 1
4265 -- Function: int gsl_sf_legendre_Pl_array (int LMAX, double X, double
4267 -- Function: int gsl_sf_legendre_Pl_deriv_array (int LMAX, double X,
4268 double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
4269 These functions compute an array of Legendre polynomials P_l(x),
4270 and optionally their derivatives dP_l(x)/dx, for l = 0, \dots,
4273 -- Function: double gsl_sf_legendre_Q0 (double X)
4274 -- Function: int gsl_sf_legendre_Q0_e (double X, gsl_sf_result *
4276 These routines compute the Legendre function Q_0(x) for x > -1, x
4279 -- Function: double gsl_sf_legendre_Q1 (double X)
4280 -- Function: int gsl_sf_legendre_Q1_e (double X, gsl_sf_result *
4282 These routines compute the Legendre function Q_1(x) for x > -1, x
4285 -- Function: double gsl_sf_legendre_Ql (int L, double X)
4286 -- Function: int gsl_sf_legendre_Ql_e (int L, double X, gsl_sf_result
4288 These routines compute the Legendre function Q_l(x) for x > -1, x
4292 File: gsl-ref.info, Node: Associated Legendre Polynomials and Spherical Harmonics, Next: Conical Functions, Prev: Legendre Polynomials, Up: Legendre Functions and Spherical Harmonics
4294 7.24.2 Associated Legendre Polynomials and Spherical Harmonics
4295 --------------------------------------------------------------
4297 The following functions compute the associated Legendre Polynomials
4298 P_l^m(x). Note that this function grows combinatorially with l and can
4299 overflow for l larger than about 150. There is no trouble for small m,
4300 but overflow occurs when m and l are both large. Rather than allow
4301 overflows, these functions refuse to calculate P_l^m(x) and return
4302 `GSL_EOVRFLW' when they can sense that l and m are too big.
4304 If you want to calculate a spherical harmonic, then _do not_ use
4305 these functions. Instead use `gsl_sf_legendre_sphPlm' below, which
4306 uses a similar recursion, but with the normalized functions.
4308 -- Function: double gsl_sf_legendre_Plm (int L, int M, double X)
4309 -- Function: int gsl_sf_legendre_Plm_e (int L, int M, double X,
4310 gsl_sf_result * RESULT)
4311 These routines compute the associated Legendre polynomial P_l^m(x)
4312 for m >= 0, l >= m, |x| <= 1.
4314 -- Function: int gsl_sf_legendre_Plm_array (int LMAX, int M, double X,
4315 double RESULT_ARRAY[])
4316 -- Function: int gsl_sf_legendre_Plm_deriv_array (int LMAX, int M,
4317 double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
4318 These functions compute an array of Legendre polynomials P_l^m(x),
4319 and optionally their derivatives dP_l^m(x)/dx, for m >= 0, l =
4320 |m|, ..., lmax, |x| <= 1.
4322 -- Function: double gsl_sf_legendre_sphPlm (int L, int M, double X)
4323 -- Function: int gsl_sf_legendre_sphPlm_e (int L, int M, double X,
4324 gsl_sf_result * RESULT)
4325 These routines compute the normalized associated Legendre
4326 polynomial $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$
4327 suitable for use in spherical harmonics. The parameters must
4328 satisfy m >= 0, l >= m, |x| <= 1. Theses routines avoid the
4329 overflows that occur for the standard normalization of P_l^m(x).
4331 -- Function: int gsl_sf_legendre_sphPlm_array (int LMAX, int M, double
4332 X, double RESULT_ARRAY[])
4333 -- Function: int gsl_sf_legendre_sphPlm_deriv_array (int LMAX, int M,
4334 double X, double RESULT_ARRAY[], double RESULT_DERIV_ARRAY[])
4335 These functions compute an array of normalized associated Legendre
4336 functions $\sqrt{(2l+1)/(4\pi)} \sqrt{(l-m)!/(l+m)!} P_l^m(x)$,
4337 and optionally their derivatives, for m >= 0, l = |m|, ..., lmax,
4340 -- Function: int gsl_sf_legendre_array_size (const int LMAX, const int
4342 This function returns the size of RESULT_ARRAY[] needed for the
4343 array versions of P_l^m(x), LMAX - M + 1. An inline version of
4344 this function is used when `HAVE_INLINE' is defined.
4347 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
4349 7.24.3 Conical Functions
4350 ------------------------
4352 The Conical Functions P^\mu_{-(1/2)+i\lambda}(x) and
4353 Q^\mu_{-(1/2)+i\lambda} are described in Abramowitz & Stegun, Section
4356 -- Function: double gsl_sf_conicalP_half (double LAMBDA, double X)
4357 -- Function: int gsl_sf_conicalP_half_e (double LAMBDA, double X,
4358 gsl_sf_result * RESULT)
4359 These routines compute the irregular Spherical Conical Function
4360 P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.
4362 -- Function: double gsl_sf_conicalP_mhalf (double LAMBDA, double X)
4363 -- Function: int gsl_sf_conicalP_mhalf_e (double LAMBDA, double X,
4364 gsl_sf_result * RESULT)
4365 These routines compute the regular Spherical Conical Function
4366 P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.
4368 -- Function: double gsl_sf_conicalP_0 (double LAMBDA, double X)
4369 -- Function: int gsl_sf_conicalP_0_e (double LAMBDA, double X,
4370 gsl_sf_result * RESULT)
4371 These routines compute the conical function P^0_{-1/2 + i
4372 \lambda}(x) for x > -1.
4374 -- Function: double gsl_sf_conicalP_1 (double LAMBDA, double X)
4375 -- Function: int gsl_sf_conicalP_1_e (double LAMBDA, double X,
4376 gsl_sf_result * RESULT)
4377 These routines compute the conical function P^1_{-1/2 + i
4378 \lambda}(x) for x > -1.
4380 -- Function: double gsl_sf_conicalP_sph_reg (int L, double LAMBDA,
4382 -- Function: int gsl_sf_conicalP_sph_reg_e (int L, double LAMBDA,
4383 double X, gsl_sf_result * RESULT)
4384 These routines compute the Regular Spherical Conical Function
4385 P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.
4387 -- Function: double gsl_sf_conicalP_cyl_reg (int M, double LAMBDA,
4389 -- Function: int gsl_sf_conicalP_cyl_reg_e (int M, double LAMBDA,
4390 double X, gsl_sf_result * RESULT)
4391 These routines compute the Regular Cylindrical Conical Function
4392 P^{-m}_{-1/2 + i \lambda}(x) for x > -1, m >= -1.
4395 File: gsl-ref.info, Node: Radial Functions for Hyperbolic Space, Prev: Conical Functions, Up: Legendre Functions and Spherical Harmonics
4397 7.24.4 Radial Functions for Hyperbolic Space
4398 --------------------------------------------
4400 The following spherical functions are specializations of Legendre
4401 functions which give the regular eigenfunctions of the Laplacian on a
4402 3-dimensional hyperbolic space H3d. Of particular interest is the flat
4403 limit, \lambda \to \infty, \eta \to 0, \lambda\eta fixed.
4405 -- Function: double gsl_sf_legendre_H3d_0 (double LAMBDA, double ETA)
4406 -- Function: int gsl_sf_legendre_H3d_0_e (double LAMBDA, double ETA,
4407 gsl_sf_result * RESULT)
4408 These routines compute the zeroth radial eigenfunction of the
4409 Laplacian on the 3-dimensional hyperbolic space,
4410 L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta))
4411 for \eta >= 0. In the flat limit this takes the form
4412 L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta).
4414 -- Function: double gsl_sf_legendre_H3d_1 (double LAMBDA, double ETA)
4415 -- Function: int gsl_sf_legendre_H3d_1_e (double LAMBDA, double ETA,
4416 gsl_sf_result * RESULT)
4417 These routines compute the first radial eigenfunction of the
4418 Laplacian on the 3-dimensional hyperbolic space,
4419 L^{H3d}_1(\lambda,\eta) := 1/\sqrt{\lambda^2 + 1} \sin(\lambda
4420 \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda
4421 \cot(\lambda\eta)) for \eta >= 0. In the flat limit this takes
4422 the form L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta).
4424 -- Function: double gsl_sf_legendre_H3d (int L, double LAMBDA, double
4426 -- Function: int gsl_sf_legendre_H3d_e (int L, double LAMBDA, double
4427 ETA, gsl_sf_result * RESULT)
4428 These routines compute the L-th radial eigenfunction of the
4429 Laplacian on the 3-dimensional hyperbolic space \eta >= 0, l >= 0.
4430 In the flat limit this takes the form L^{H3d}_l(\lambda,\eta) =
4433 -- Function: int gsl_sf_legendre_H3d_array (int LMAX, double LAMBDA,
4434 double ETA, double RESULT_ARRAY[])
4435 This function computes an array of radial eigenfunctions
4436 L^{H3d}_l(\lambda, \eta) for 0 <= l <= lmax.
4439 File: gsl-ref.info, Node: Logarithm and Related Functions, Next: Mathieu Functions, Prev: Legendre Functions and Spherical Harmonics, Up: Special Functions
4441 7.25 Logarithm and Related Functions
4442 ====================================
4444 Information on the properties of the Logarithm function can be found in
4445 Abramowitz & Stegun, Chapter 4. The functions described in this section
4446 are declared in the header file `gsl_sf_log.h'.
4448 -- Function: double gsl_sf_log (double X)
4449 -- Function: int gsl_sf_log_e (double X, gsl_sf_result * RESULT)
4450 These routines compute the logarithm of X, \log(x), for x > 0.
4452 -- Function: double gsl_sf_log_abs (double X)
4453 -- Function: int gsl_sf_log_abs_e (double X, gsl_sf_result * RESULT)
4454 These routines compute the logarithm of the magnitude of X,
4455 \log(|x|), for x \ne 0.
4457 -- Function: int gsl_sf_complex_log_e (double ZR, double ZI,
4458 gsl_sf_result * LNR, gsl_sf_result * THETA)
4459 This routine computes the complex logarithm of z = z_r + i z_i.
4460 The results are returned as LNR, THETA such that \exp(lnr + i
4461 \theta) = z_r + i z_i, where \theta lies in the range [-\pi,\pi].
4463 -- Function: double gsl_sf_log_1plusx (double X)
4464 -- Function: int gsl_sf_log_1plusx_e (double X, gsl_sf_result * RESULT)
4465 These routines compute \log(1 + x) for x > -1 using an algorithm
4466 that is accurate for small x.
4468 -- Function: double gsl_sf_log_1plusx_mx (double X)
4469 -- Function: int gsl_sf_log_1plusx_mx_e (double X, gsl_sf_result *
4471 These routines compute \log(1 + x) - x for x > -1 using an
4472 algorithm that is accurate for small x.
4475 File: gsl-ref.info, Node: Mathieu Functions, Next: Power Function, Prev: Logarithm and Related Functions, Up: Special Functions
4477 7.26 Mathieu Functions
4478 ======================
4480 The routines described in this section compute the angular and radial
4481 Mathieu functions, and their characteristic values. Mathieu functions
4482 are the solutions of the following two differential equations:
4484 d^2y/dv^2 + (a - 2q\cos 2v)y = 0
4485 d^2f/du^2 - (a - 2q\cosh 2u)f = 0
4487 The angular Mathieu functions ce_r(x,q), se_r(x,q) are the even and odd
4488 periodic solutions of the first equation, which is known as Mathieu's
4489 equation. These exist only for the discrete sequence of characteristic
4490 values a=a_r(q) (even-periodic) and a=b_r(q) (odd-periodic).
4492 The radial Mathieu functions Mc^{(j)}_{r}(z,q), Ms^{(j)}_{r}(z,q)
4493 are the solutions of the second equation, which is referred to as
4494 Mathieu's modified equation. The radial Mathieu functions of the
4495 first, second, third and fourth kind are denoted by the parameter j,
4496 which takes the value 1, 2, 3 or 4.
4498 For more information on the Mathieu functions, see Abramowitz and
4499 Stegun, Chapter 20. These functions are defined in the header file
4504 * Mathieu Function Workspace::
4505 * Mathieu Function Characteristic Values::
4506 * Angular Mathieu Functions::
4507 * Radial Mathieu Functions::
4510 File: gsl-ref.info, Node: Mathieu Function Workspace, Next: Mathieu Function Characteristic Values, Up: Mathieu Functions
4512 7.26.1 Mathieu Function Workspace
4513 ---------------------------------
4515 The Mathieu functions can be computed for a single order or for
4516 multiple orders, using array-based routines. The array-based routines
4517 require a preallocated workspace.
4519 -- Function: gsl_sf_mathieu_workspace * gsl_sf_mathieu_alloc (size_t
4521 This function returns a workspace for the array versions of the
4522 Mathieu routines. The arguments N and QMAX specify the maximum
4523 order and q-value of Mathieu functions which can be computed with
4527 -- Function: void gsl_sf_mathieu_free (gsl_sf_mathieu_workspace *WORK)
4528 This function frees the workspace WORK.
4531 File: gsl-ref.info, Node: Mathieu Function Characteristic Values, Next: Angular Mathieu Functions, Prev: Mathieu Function Workspace, Up: Mathieu Functions
4533 7.26.2 Mathieu Function Characteristic Values
4534 ---------------------------------------------
4536 -- Function: int gsl_sf_mathieu_a (int N, double Q, gsl_sf_result
4538 -- Function: int gsl_sf_mathieu_b (int N, double Q, gsl_sf_result
4540 These routines compute the characteristic values a_n(q), b_n(q) of
4541 the Mathieu functions ce_n(q,x) and se_n(q,x), respectively.
4543 -- Function: int gsl_sf_mathieu_a_array (int ORDER_MIN, int ORDER_MAX,
4544 double Q, gsl_sf_mathieu_workspace *WORK, double
4546 -- Function: int gsl_sf_mathieu_b_array (int ORDER_MIN, int ORDER_MAX,
4547 double Q, gsl_sf_mathieu_workspace *WORK, double
4549 These routines compute a series of Mathieu characteristic values
4550 a_n(q), b_n(q) for n from ORDER_MIN to ORDER_MAX inclusive,
4551 storing the results in the array RESULT_ARRAY.
4554 File: gsl-ref.info, Node: Angular Mathieu Functions, Next: Radial Mathieu Functions, Prev: Mathieu Function Characteristic Values, Up: Mathieu Functions
4556 7.26.3 Angular Mathieu Functions
4557 --------------------------------
4559 -- Function: int gsl_sf_mathieu_ce (int N, double Q, double X,
4560 gsl_sf_result *RESULT)
4561 -- Function: int gsl_sf_mathieu_se (int N, double Q, double X,
4562 gsl_sf_result *RESULT)
4563 These routines compute the angular Mathieu functions ce_n(q,x) and
4564 se_n(q,x), respectively.
4566 -- Function: int gsl_sf_mathieu_ce_array (int NMIN, int NMAX, double
4567 Q, double X, gsl_sf_mathieu_workspace *WORK, double
4569 -- Function: int gsl_sf_mathieu_se_array (int NMIN, int NMAX, double
4570 Q, double X, gsl_sf_mathieu_workspace *WORK, double
4572 These routines compute a series of the angular Mathieu functions
4573 ce_n(q,x) and se_n(q,x) of order n from NMIN to NMAX inclusive,
4574 storing the results in the array RESULT_ARRAY.
4577 File: gsl-ref.info, Node: Radial Mathieu Functions, Prev: Angular Mathieu Functions, Up: Mathieu Functions
4579 7.26.4 Radial Mathieu Functions
4580 -------------------------------
4582 -- Function: int gsl_sf_mathieu_Mc (int J, int N, double Q, double X,
4583 gsl_sf_result *RESULT)
4584 -- Function: int gsl_sf_mathieu_Ms (int J, int N, double Q, double X,
4585 gsl_sf_result *RESULT)
4586 These routines compute the radial J-th kind Mathieu functions
4587 Mc_n^{(j)}(q,x) and Ms_n^{(j)}(q,x) of order N.
4589 The allowed values of J are 1 and 2. The functions for j = 3,4
4590 can be computed as M_n^{(3)} = M_n^{(1)} + iM_n^{(2)} and
4591 M_n^{(4)} = M_n^{(1)} - iM_n^{(2)}, where M_n^{(j)} = Mc_n^{(j)} or
4594 -- Function: int gsl_sf_mathieu_Mc_array (int J, int NMIN, int NMAX,
4595 double Q, double X, gsl_sf_mathieu_workspace *WORK, double
4597 -- Function: int gsl_sf_mathieu_Ms_array (int J, int NMIN, int NMAX,
4598 double Q, double X, gsl_sf_mathieu_workspace *WORK, double
4600 These routines compute a series of the radial Mathieu functions of
4601 kind J, with order from NMIN to NMAX inclusive, storing the
4602 results in the array RESULT_ARRAY.
4605 File: gsl-ref.info, Node: Power Function, Next: Psi (Digamma) Function, Prev: Mathieu Functions, Up: Special Functions
4610 The following functions are equivalent to the function `gsl_pow_int'
4611 (*note Small integer powers::) with an error estimate. These functions
4612 are declared in the header file `gsl_sf_pow_int.h'.
4614 -- Function: double gsl_sf_pow_int (double X, int N)
4615 -- Function: int gsl_sf_pow_int_e (double X, int N, gsl_sf_result *
4617 These routines compute the power x^n for integer N. The power is
4618 computed using the minimum number of multiplications. For example,
4619 x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications.
4620 For reasons of efficiency, these functions do not check for
4621 overflow or underflow conditions.
4623 #include <gsl/gsl_sf_pow_int.h>
4624 /* compute 3.0**12 */
4625 double y = gsl_sf_pow_int(3.0, 12);
4628 File: gsl-ref.info, Node: Psi (Digamma) Function, Next: Synchrotron Functions, Prev: Power Function, Up: Special Functions
4630 7.28 Psi (Digamma) Function
4631 ===========================
4633 The polygamma functions of order n are defined by
4635 \psi^{(n)}(x) = (d/dx)^n \psi(x) = (d/dx)^{n+1} \log(\Gamma(x))
4637 where \psi(x) = \Gamma'(x)/\Gamma(x) is known as the digamma function.
4638 These functions are declared in the header file `gsl_sf_psi.h'.
4642 * Digamma Function::
4643 * Trigamma Function::
4644 * Polygamma Function::
4647 File: gsl-ref.info, Node: Digamma Function, Next: Trigamma Function, Up: Psi (Digamma) Function
4649 7.28.1 Digamma Function
4650 -----------------------
4652 -- Function: double gsl_sf_psi_int (int N)
4653 -- Function: int gsl_sf_psi_int_e (int N, gsl_sf_result * RESULT)
4654 These routines compute the digamma function \psi(n) for positive
4655 integer N. The digamma function is also called the Psi function.
4657 -- Function: double gsl_sf_psi (double X)
4658 -- Function: int gsl_sf_psi_e (double X, gsl_sf_result * RESULT)
4659 These routines compute the digamma function \psi(x) for general x,
4662 -- Function: double gsl_sf_psi_1piy (double Y)
4663 -- Function: int gsl_sf_psi_1piy_e (double Y, gsl_sf_result * RESULT)
4664 These routines compute the real part of the digamma function on
4665 the line 1+i y, \Re[\psi(1 + i y)].
4668 File: gsl-ref.info, Node: Trigamma Function, Next: Polygamma Function, Prev: Digamma Function, Up: Psi (Digamma) Function
4670 7.28.2 Trigamma Function
4671 ------------------------
4673 -- Function: double gsl_sf_psi_1_int (int N)
4674 -- Function: int gsl_sf_psi_1_int_e (int N, gsl_sf_result * RESULT)
4675 These routines compute the Trigamma function \psi'(n) for positive
4678 -- Function: double gsl_sf_psi_1 (double X)
4679 -- Function: int gsl_sf_psi_1_e (double X, gsl_sf_result * RESULT)
4680 These routines compute the Trigamma function \psi'(x) for general
4684 File: gsl-ref.info, Node: Polygamma Function, Prev: Trigamma Function, Up: Psi (Digamma) Function
4686 7.28.3 Polygamma Function
4687 -------------------------
4689 -- Function: double gsl_sf_psi_n (int N, double X)
4690 -- Function: int gsl_sf_psi_n_e (int N, double X, gsl_sf_result *
4692 These routines compute the polygamma function \psi^{(n)}(x) for n
4696 File: gsl-ref.info, Node: Synchrotron Functions, Next: Transport Functions, Prev: Psi (Digamma) Function, Up: Special Functions
4698 7.29 Synchrotron Functions
4699 ==========================
4701 The functions described in this section are declared in the header file
4702 `gsl_sf_synchrotron.h'.
4704 -- Function: double gsl_sf_synchrotron_1 (double X)
4705 -- Function: int gsl_sf_synchrotron_1_e (double X, gsl_sf_result *
4707 These routines compute the first synchrotron function x
4708 \int_x^\infty dt K_{5/3}(t) for x >= 0.
4710 -- Function: double gsl_sf_synchrotron_2 (double X)
4711 -- Function: int gsl_sf_synchrotron_2_e (double X, gsl_sf_result *
4713 These routines compute the second synchrotron function x
4714 K_{2/3}(x) for x >= 0.
4717 File: gsl-ref.info, Node: Transport Functions, Next: Trigonometric Functions, Prev: Synchrotron Functions, Up: Special Functions
4719 7.30 Transport Functions
4720 ========================
4722 The transport functions J(n,x) are defined by the integral
4723 representations J(n,x) := \int_0^x dt t^n e^t /(e^t - 1)^2. They are
4724 declared in the header file `gsl_sf_transport.h'.
4726 -- Function: double gsl_sf_transport_2 (double X)
4727 -- Function: int gsl_sf_transport_2_e (double X, gsl_sf_result *
4729 These routines compute the transport function J(2,x).
4731 -- Function: double gsl_sf_transport_3 (double X)
4732 -- Function: int gsl_sf_transport_3_e (double X, gsl_sf_result *
4734 These routines compute the transport function J(3,x).
4736 -- Function: double gsl_sf_transport_4 (double X)
4737 -- Function: int gsl_sf_transport_4_e (double X, gsl_sf_result *
4739 These routines compute the transport function J(4,x).
4741 -- Function: double gsl_sf_transport_5 (double X)
4742 -- Function: int gsl_sf_transport_5_e (double X, gsl_sf_result *
4744 These routines compute the transport function J(5,x).
4747 File: gsl-ref.info, Node: Trigonometric Functions, Next: Zeta Functions, Prev: Transport Functions, Up: Special Functions
4749 7.31 Trigonometric Functions
4750 ============================
4752 The library includes its own trigonometric functions in order to provide
4753 consistency across platforms and reliable error estimates. These
4754 functions are declared in the header file `gsl_sf_trig.h'.
4758 * Circular Trigonometric Functions::
4759 * Trigonometric Functions for Complex Arguments::
4760 * Hyperbolic Trigonometric Functions::
4761 * Conversion Functions::
4762 * Restriction Functions::
4763 * Trigonometric Functions With Error Estimates::
4766 File: gsl-ref.info, Node: Circular Trigonometric Functions, Next: Trigonometric Functions for Complex Arguments, Up: Trigonometric Functions
4768 7.31.1 Circular Trigonometric Functions
4769 ---------------------------------------
4771 -- Function: double gsl_sf_sin (double X)
4772 -- Function: int gsl_sf_sin_e (double X, gsl_sf_result * RESULT)
4773 These routines compute the sine function \sin(x).
4775 -- Function: double gsl_sf_cos (double X)
4776 -- Function: int gsl_sf_cos_e (double X, gsl_sf_result * RESULT)
4777 These routines compute the cosine function \cos(x).
4779 -- Function: double gsl_sf_hypot (double X, double Y)
4780 -- Function: int gsl_sf_hypot_e (double X, double Y, gsl_sf_result *
4782 These routines compute the hypotenuse function \sqrt{x^2 + y^2}
4783 avoiding overflow and underflow.
4785 -- Function: double gsl_sf_sinc (double X)
4786 -- Function: int gsl_sf_sinc_e (double X, gsl_sf_result * RESULT)
4787 These routines compute \sinc(x) = \sin(\pi x) / (\pi x) for any
4791 File: gsl-ref.info, Node: Trigonometric Functions for Complex Arguments, Next: Hyperbolic Trigonometric Functions, Prev: Circular Trigonometric Functions, Up: Trigonometric Functions
4793 7.31.2 Trigonometric Functions for Complex Arguments
4794 ----------------------------------------------------
4796 -- Function: int gsl_sf_complex_sin_e (double ZR, double ZI,
4797 gsl_sf_result * SZR, gsl_sf_result * SZI)
4798 This function computes the complex sine, \sin(z_r + i z_i) storing
4799 the real and imaginary parts in SZR, SZI.
4801 -- Function: int gsl_sf_complex_cos_e (double ZR, double ZI,
4802 gsl_sf_result * CZR, gsl_sf_result * CZI)
4803 This function computes the complex cosine, \cos(z_r + i z_i)
4804 storing the real and imaginary parts in SZR, SZI.
4806 -- Function: int gsl_sf_complex_logsin_e (double ZR, double ZI,
4807 gsl_sf_result * LSZR, gsl_sf_result * LSZI)
4808 This function computes the logarithm of the complex sine,
4809 \log(\sin(z_r + i z_i)) storing the real and imaginary parts in
4813 File: gsl-ref.info, Node: Hyperbolic Trigonometric Functions, Next: Conversion Functions, Prev: Trigonometric Functions for Complex Arguments, Up: Trigonometric Functions
4815 7.31.3 Hyperbolic Trigonometric Functions
4816 -----------------------------------------
4818 -- Function: double gsl_sf_lnsinh (double X)
4819 -- Function: int gsl_sf_lnsinh_e (double X, gsl_sf_result * RESULT)
4820 These routines compute \log(\sinh(x)) for x > 0.
4822 -- Function: double gsl_sf_lncosh (double X)
4823 -- Function: int gsl_sf_lncosh_e (double X, gsl_sf_result * RESULT)
4824 These routines compute \log(\cosh(x)) for any X.
4827 File: gsl-ref.info, Node: Conversion Functions, Next: Restriction Functions, Prev: Hyperbolic Trigonometric Functions, Up: Trigonometric Functions
4829 7.31.4 Conversion Functions
4830 ---------------------------
4832 -- Function: int gsl_sf_polar_to_rect (double R, double THETA,
4833 gsl_sf_result * X, gsl_sf_result * Y);
4834 This function converts the polar coordinates (R,THETA) to
4835 rectilinear coordinates (X,Y), x = r\cos(\theta), y =
4838 -- Function: int gsl_sf_rect_to_polar (double X, double Y,
4839 gsl_sf_result * R, gsl_sf_result * THETA)
4840 This function converts the rectilinear coordinates (X,Y) to polar
4841 coordinates (R,THETA), such that x = r\cos(\theta), y =
4842 r\sin(\theta). The argument THETA lies in the range [-\pi, \pi].
4845 File: gsl-ref.info, Node: Restriction Functions, Next: Trigonometric Functions With Error Estimates, Prev: Conversion Functions, Up: Trigonometric Functions
4847 7.31.5 Restriction Functions
4848 ----------------------------
4850 -- Function: double gsl_sf_angle_restrict_symm (double THETA)
4851 -- Function: int gsl_sf_angle_restrict_symm_e (double * THETA)
4852 These routines force the angle THETA to lie in the range
4855 Note that the mathematical value of \pi is slightly greater than
4856 `M_PI', so the machine numbers `M_PI' and `-M_PI' are included in
4859 -- Function: double gsl_sf_angle_restrict_pos (double THETA)
4860 -- Function: int gsl_sf_angle_restrict_pos_e (double * THETA)
4861 These routines force the angle THETA to lie in the range [0, 2\pi).
4863 Note that the mathematical value of 2\pi is slightly greater than
4864 `2*M_PI', so the machine number `2*M_PI' is included in the range.
4868 File: gsl-ref.info, Node: Trigonometric Functions With Error Estimates, Prev: Restriction Functions, Up: Trigonometric Functions
4870 7.31.6 Trigonometric Functions With Error Estimates
4871 ---------------------------------------------------
4873 -- Function: int gsl_sf_sin_err_e (double X, double DX, gsl_sf_result
4875 This routine computes the sine of an angle X with an associated
4876 absolute error DX, \sin(x \pm dx). Note that this function is
4877 provided in the error-handling form only since its purpose is to
4878 compute the propagated error.
4880 -- Function: int gsl_sf_cos_err_e (double X, double DX, gsl_sf_result
4882 This routine computes the cosine of an angle X with an associated
4883 absolute error DX, \cos(x \pm dx). Note that this function is
4884 provided in the error-handling form only since its purpose is to
4885 compute the propagated error.
4888 File: gsl-ref.info, Node: Zeta Functions, Next: Special Functions Examples, Prev: Trigonometric Functions, Up: Special Functions
4893 The Riemann zeta function is defined in Abramowitz & Stegun, Section
4894 23.2. The functions described in this section are declared in the
4895 header file `gsl_sf_zeta.h'.
4899 * Riemann Zeta Function::
4900 * Riemann Zeta Function Minus One::
4901 * Hurwitz Zeta Function::
4905 File: gsl-ref.info, Node: Riemann Zeta Function, Next: Riemann Zeta Function Minus One, Up: Zeta Functions
4907 7.32.1 Riemann Zeta Function
4908 ----------------------------
4910 The Riemann zeta function is defined by the infinite sum \zeta(s) =
4911 \sum_{k=1}^\infty k^{-s}.
4913 -- Function: double gsl_sf_zeta_int (int N)
4914 -- Function: int gsl_sf_zeta_int_e (int N, gsl_sf_result * RESULT)
4915 These routines compute the Riemann zeta function \zeta(n) for
4918 -- Function: double gsl_sf_zeta (double S)
4919 -- Function: int gsl_sf_zeta_e (double S, gsl_sf_result * RESULT)
4920 These routines compute the Riemann zeta function \zeta(s) for
4921 arbitrary S, s \ne 1.
4924 File: gsl-ref.info, Node: Riemann Zeta Function Minus One, Next: Hurwitz Zeta Function, Prev: Riemann Zeta Function, Up: Zeta Functions
4926 7.32.2 Riemann Zeta Function Minus One
4927 --------------------------------------
4929 For large positive argument, the Riemann zeta function approaches one.
4930 In this region the fractional part is interesting, and therefore we
4931 need a function to evaluate it explicitly.
4933 -- Function: double gsl_sf_zetam1_int (int N)
4934 -- Function: int gsl_sf_zetam1_int_e (int N, gsl_sf_result * RESULT)
4935 These routines compute \zeta(n) - 1 for integer N, n \ne 1.
4937 -- Function: double gsl_sf_zetam1 (double S)
4938 -- Function: int gsl_sf_zetam1_e (double S, gsl_sf_result * RESULT)
4939 These routines compute \zeta(s) - 1 for arbitrary S, s \ne 1.
4942 File: gsl-ref.info, Node: Hurwitz Zeta Function, Next: Eta Function, Prev: Riemann Zeta Function Minus One, Up: Zeta Functions
4944 7.32.3 Hurwitz Zeta Function
4945 ----------------------------
4947 The Hurwitz zeta function is defined by \zeta(s,q) = \sum_0^\infty
4950 -- Function: double gsl_sf_hzeta (double S, double Q)
4951 -- Function: int gsl_sf_hzeta_e (double S, double Q, gsl_sf_result *
4953 These routines compute the Hurwitz zeta function \zeta(s,q) for s
4957 File: gsl-ref.info, Node: Eta Function, Prev: Hurwitz Zeta Function, Up: Zeta Functions
4962 The eta function is defined by \eta(s) = (1-2^{1-s}) \zeta(s).
4964 -- Function: double gsl_sf_eta_int (int N)
4965 -- Function: int gsl_sf_eta_int_e (int N, gsl_sf_result * RESULT)
4966 These routines compute the eta function \eta(n) for integer N.
4968 -- Function: double gsl_sf_eta (double S)
4969 -- Function: int gsl_sf_eta_e (double S, gsl_sf_result * RESULT)
4970 These routines compute the eta function \eta(s) for arbitrary S.
4973 File: gsl-ref.info, Node: Special Functions Examples, Next: Special Functions References and Further Reading, Prev: Zeta Functions, Up: Special Functions
4978 The following example demonstrates the use of the error handling form of
4979 the special functions, in this case to compute the Bessel function
4983 #include <gsl/gsl_errno.h>
4984 #include <gsl/gsl_sf_bessel.h>
4990 gsl_sf_result result;
4992 double expected = -0.17759677131433830434739701;
4994 int status = gsl_sf_bessel_J0_e (x, &result);
4996 printf ("status = %s\n", gsl_strerror(status));
4997 printf ("J0(5.0) = %.18f\n"
4999 result.val, result.err);
5000 printf ("exact = %.18f\n", expected);
5004 Here are the results of running the program,
5008 J0(5.0) = -0.177596771314338292
5009 +/- 0.000000000000000193
5010 exact = -0.177596771314338292
5012 The next program computes the same quantity using the natural form of
5013 the function. In this case the error term RESULT.ERR and return status
5017 #include <gsl/gsl_sf_bessel.h>
5023 double expected = -0.17759677131433830434739701;
5025 double y = gsl_sf_bessel_J0 (x);
5027 printf ("J0(5.0) = %.18f\n", y);
5028 printf ("exact = %.18f\n", expected);
5032 The results of the function are the same,
5035 J0(5.0) = -0.177596771314338292
5036 exact = -0.177596771314338292
5039 File: gsl-ref.info, Node: Special Functions References and Further Reading, Prev: Special Functions Examples, Up: Special Functions
5041 7.34 References and Further Reading
5042 ===================================
5044 The library follows the conventions of `Abramowitz & Stegun' where
5046 Abramowitz & Stegun (eds.), `Handbook of Mathematical Functions'
5048 The following papers contain information on the algorithms used to
5049 compute the special functions,
5050 MISCFUN: A software package to compute uncommon special functions.
5051 `ACM Trans. Math. Soft.', vol. 22, 1996, 288-301
5053 G.N. Watson, A Treatise on the Theory of Bessel Functions, 2nd
5054 Edition (Cambridge University Press, 1944).
5056 G. Nemeth, Mathematical Approximations of Special Functions, Nova
5057 Science Publishers, ISBN 1-56072-052-2
5059 B.C. Carlson, Special Functions of Applied Mathematics (1977)
5061 W.J. Thompson, Atlas for Computing Mathematical Functions, John
5062 Wiley & Sons, New York (1997).
5064 Y.Y. Luke, Algorithms for the Computation of Mathematical
5065 Functions, Academic Press, New York (1977).
5069 File: gsl-ref.info, Node: Vectors and Matrices, Next: Permutations, Prev: Special Functions, Up: Top
5071 8 Vectors and Matrices
5072 **********************
5074 The functions described in this chapter provide a simple vector and
5075 matrix interface to ordinary C arrays. The memory management of these
5076 arrays is implemented using a single underlying type, known as a block.
5077 By writing your functions in terms of vectors and matrices you can pass
5078 a single structure containing both data and dimensions as an argument
5079 without needing additional function parameters. The structures are
5080 compatible with the vector and matrix formats used by BLAS routines.
5088 * Vector and Matrix References and Further Reading::
5091 File: gsl-ref.info, Node: Data types, Next: Blocks, Up: Vectors and Matrices
5096 All the functions are available for each of the standard data-types.
5097 The versions for `double' have the prefix `gsl_block', `gsl_vector' and
5098 `gsl_matrix'. Similarly the versions for single-precision `float'
5099 arrays have the prefix `gsl_block_float', `gsl_vector_float' and
5100 `gsl_matrix_float'. The full list of available types is given below,
5103 gsl_block_float float
5104 gsl_block_long_double long double
5106 gsl_block_uint unsigned int
5108 gsl_block_ulong unsigned long
5109 gsl_block_short short
5110 gsl_block_ushort unsigned short
5112 gsl_block_uchar unsigned char
5113 gsl_block_complex complex double
5114 gsl_block_complex_float complex float
5115 gsl_block_complex_long_double complex long double
5117 Corresponding types exist for the `gsl_vector' and `gsl_matrix'
5121 File: gsl-ref.info, Node: Blocks, Next: Vectors, Prev: Data types, Up: Vectors and Matrices
5126 For consistency all memory is allocated through a `gsl_block'
5127 structure. The structure contains two components, the size of an area
5128 of memory and a pointer to the memory. The `gsl_block' structure looks
5137 Vectors and matrices are made by "slicing" an underlying block. A slice
5138 is a set of elements formed from an initial offset and a combination of
5139 indices and step-sizes. In the case of a matrix the step-size for the
5140 column index represents the row-length. The step-size for a vector is
5141 known as the "stride".
5143 The functions for allocating and deallocating blocks are defined in
5148 * Block allocation::
5149 * Reading and writing blocks::
5150 * Example programs for blocks::
5153 File: gsl-ref.info, Node: Block allocation, Next: Reading and writing blocks, Up: Blocks
5155 8.2.1 Block allocation
5156 ----------------------
5158 The functions for allocating memory to a block follow the style of
5159 `malloc' and `free'. In addition they also perform their own error
5160 checking. If there is insufficient memory available to allocate a
5161 block then the functions call the GSL error handler (with an error
5162 number of `GSL_ENOMEM') in addition to returning a null pointer. Thus
5163 if you use the library error handler to abort your program then it
5164 isn't necessary to check every `alloc'.
5166 -- Function: gsl_block * gsl_block_alloc (size_t N)
5167 This function allocates memory for a block of N double-precision
5168 elements, returning a pointer to the block struct. The block is
5169 not initialized and so the values of its elements are undefined.
5170 Use the function `gsl_block_calloc' if you want to ensure that all
5171 the elements are initialized to zero.
5173 A null pointer is returned if insufficient memory is available to
5176 -- Function: gsl_block * gsl_block_calloc (size_t N)
5177 This function allocates memory for a block and initializes all the
5178 elements of the block to zero.
5180 -- Function: void gsl_block_free (gsl_block * B)
5181 This function frees the memory used by a block B previously
5182 allocated with `gsl_block_alloc' or `gsl_block_calloc'. The block
5183 B must be a valid block object (a null pointer is not allowed).
5186 File: gsl-ref.info, Node: Reading and writing blocks, Next: Example programs for blocks, Prev: Block allocation, Up: Blocks
5188 8.2.2 Reading and writing blocks
5189 --------------------------------
5191 The library provides functions for reading and writing blocks to a file
5192 as binary data or formatted text.
5194 -- Function: int gsl_block_fwrite (FILE * STREAM, const gsl_block * B)
5195 This function writes the elements of the block B to the stream
5196 STREAM in binary format. The return value is 0 for success and
5197 `GSL_EFAILED' if there was a problem writing to the file. Since
5198 the data is written in the native binary format it may not be
5199 portable between different architectures.
5201 -- Function: int gsl_block_fread (FILE * STREAM, gsl_block * B)
5202 This function reads into the block B from the open stream STREAM
5203 in binary format. The block B must be preallocated with the
5204 correct length since the function uses the size of B to determine
5205 how many bytes to read. The return value is 0 for success and
5206 `GSL_EFAILED' if there was a problem reading from the file. The
5207 data is assumed to have been written in the native binary format
5208 on the same architecture.
5210 -- Function: int gsl_block_fprintf (FILE * STREAM, const gsl_block *
5211 B, const char * FORMAT)
5212 This function writes the elements of the block B line-by-line to
5213 the stream STREAM using the format specifier FORMAT, which should
5214 be one of the `%g', `%e' or `%f' formats for floating point
5215 numbers and `%d' for integers. The function returns 0 for success
5216 and `GSL_EFAILED' if there was a problem writing to the file.
5218 -- Function: int gsl_block_fscanf (FILE * STREAM, gsl_block * B)
5219 This function reads formatted data from the stream STREAM into the
5220 block B. The block B must be preallocated with the correct length
5221 since the function uses the size of B to determine how many
5222 numbers to read. The function returns 0 for success and
5223 `GSL_EFAILED' if there was a problem reading from the file.
5226 File: gsl-ref.info, Node: Example programs for blocks, Prev: Reading and writing blocks, Up: Blocks
5228 8.2.3 Example programs for blocks
5229 ---------------------------------
5231 The following program shows how to allocate a block,
5234 #include <gsl/gsl_block.h>
5239 gsl_block * b = gsl_block_alloc (100);
5241 printf ("length of block = %u\n", b->size);
5242 printf ("block data address = %#x\n", b->data);
5248 Here is the output from the program,
5250 length of block = 100
5251 block data address = 0x804b0d8
5254 File: gsl-ref.info, Node: Vectors, Next: Matrices, Prev: Blocks, Up: Vectors and Matrices
5259 Vectors are defined by a `gsl_vector' structure which describes a slice
5260 of a block. Different vectors can be created which point to the same
5261 block. A vector slice is a set of equally-spaced elements of an area
5264 The `gsl_vector' structure contains five components, the "size", the
5265 "stride", a pointer to the memory where the elements are stored, DATA,
5266 a pointer to the block owned by the vector, BLOCK, if any, and an
5267 ownership flag, OWNER. The structure is very simple and looks like
5279 The SIZE is simply the number of vector elements. The range of valid
5280 indices runs from 0 to `size-1'. The STRIDE is the step-size from one
5281 element to the next in physical memory, measured in units of the
5282 appropriate datatype. The pointer DATA gives the location of the first
5283 element of the vector in memory. The pointer BLOCK stores the location
5284 of the memory block in which the vector elements are located (if any).
5285 If the vector owns this block then the OWNER field is set to one and
5286 the block will be deallocated when the vector is freed. If the vector
5287 points to a block owned by another object then the OWNER field is zero
5288 and any underlying block will not be deallocated with the vector.
5290 The functions for allocating and accessing vectors are defined in
5295 * Vector allocation::
5296 * Accessing vector elements::
5297 * Initializing vector elements::
5298 * Reading and writing vectors::
5301 * Exchanging elements::
5302 * Vector operations::
5303 * Finding maximum and minimum elements of vectors::
5304 * Vector properties::
5305 * Example programs for vectors::
5308 File: gsl-ref.info, Node: Vector allocation, Next: Accessing vector elements, Up: Vectors
5310 8.3.1 Vector allocation
5311 -----------------------
5313 The functions for allocating memory to a vector follow the style of
5314 `malloc' and `free'. In addition they also perform their own error
5315 checking. If there is insufficient memory available to allocate a
5316 vector then the functions call the GSL error handler (with an error
5317 number of `GSL_ENOMEM') in addition to returning a null pointer. Thus
5318 if you use the library error handler to abort your program then it
5319 isn't necessary to check every `alloc'.
5321 -- Function: gsl_vector * gsl_vector_alloc (size_t N)
5322 This function creates a vector of length N, returning a pointer to
5323 a newly initialized vector struct. A new block is allocated for the
5324 elements of the vector, and stored in the BLOCK component of the
5325 vector struct. The block is "owned" by the vector, and will be
5326 deallocated when the vector is deallocated.
5328 -- Function: gsl_vector * gsl_vector_calloc (size_t N)
5329 This function allocates memory for a vector of length N and
5330 initializes all the elements of the vector to zero.
5332 -- Function: void gsl_vector_free (gsl_vector * V)
5333 This function frees a previously allocated vector V. If the
5334 vector was created using `gsl_vector_alloc' then the block
5335 underlying the vector will also be deallocated. If the vector has
5336 been created from another object then the memory is still owned by
5337 that object and will not be deallocated. The vector V must be a
5338 valid vector object (a null pointer is not allowed).
5341 File: gsl-ref.info, Node: Accessing vector elements, Next: Initializing vector elements, Prev: Vector allocation, Up: Vectors
5343 8.3.2 Accessing vector elements
5344 -------------------------------
5346 Unlike FORTRAN compilers, C compilers do not usually provide support
5347 for range checking of vectors and matrices. Range checking is
5348 available in the GNU C Compiler bounds-checking extension, but it is not
5349 part of the default installation of GCC. The functions
5350 `gsl_vector_get' and `gsl_vector_set' can perform portable range
5351 checking for you and report an error if you attempt to access elements
5352 outside the allowed range.
5354 The functions for accessing the elements of a vector or matrix are
5355 defined in `gsl_vector.h' and declared `extern inline' to eliminate
5356 function-call overhead. You must compile your program with the macro
5357 `HAVE_INLINE' defined to use these functions.
5359 If necessary you can turn off range checking completely without
5360 modifying any source files by recompiling your program with the
5361 preprocessor definition `GSL_RANGE_CHECK_OFF'. Provided your compiler
5362 supports inline functions the effect of turning off range checking is
5363 to replace calls to `gsl_vector_get(v,i)' by `v->data[i*v->stride]' and
5364 calls to `gsl_vector_set(v,i,x)' by `v->data[i*v->stride]=x'. Thus
5365 there should be no performance penalty for using the range checking
5366 functions when range checking is turned off.
5368 -- Function: double gsl_vector_get (const gsl_vector * V, size_t I)
5369 This function returns the I-th element of a vector V. If I lies
5370 outside the allowed range of 0 to N-1 then the error handler is
5371 invoked and 0 is returned. An inline version of this function is
5372 used when `HAVE_INLINE' is defined.
5374 -- Function: void gsl_vector_set (gsl_vector * V, size_t I, double X)
5375 This function sets the value of the I-th element of a vector V to
5376 X. If I lies outside the allowed range of 0 to N-1 then the error
5377 handler is invoked. An inline version of this function is used
5378 when `HAVE_INLINE' is defined.
5380 -- Function: double * gsl_vector_ptr (gsl_vector * V, size_t I)
5381 -- Function: const double * gsl_vector_const_ptr (const gsl_vector *
5383 These functions return a pointer to the I-th element of a vector
5384 V. If I lies outside the allowed range of 0 to N-1 then the error
5385 handler is invoked and a null pointer is returned. Inline
5386 versions of these functions are used when `HAVE_INLINE' is defined.
5389 File: gsl-ref.info, Node: Initializing vector elements, Next: Reading and writing vectors, Prev: Accessing vector elements, Up: Vectors
5391 8.3.3 Initializing vector elements
5392 ----------------------------------
5394 -- Function: void gsl_vector_set_all (gsl_vector * V, double X)
5395 This function sets all the elements of the vector V to the value X.
5397 -- Function: void gsl_vector_set_zero (gsl_vector * V)
5398 This function sets all the elements of the vector V to zero.
5400 -- Function: int gsl_vector_set_basis (gsl_vector * V, size_t I)
5401 This function makes a basis vector by setting all the elements of
5402 the vector V to zero except for the I-th element which is set to
5406 File: gsl-ref.info, Node: Reading and writing vectors, Next: Vector views, Prev: Initializing vector elements, Up: Vectors
5408 8.3.4 Reading and writing vectors
5409 ---------------------------------
5411 The library provides functions for reading and writing vectors to a file
5412 as binary data or formatted text.
5414 -- Function: int gsl_vector_fwrite (FILE * STREAM, const gsl_vector *
5416 This function writes the elements of the vector V to the stream
5417 STREAM in binary format. The return value is 0 for success and
5418 `GSL_EFAILED' if there was a problem writing to the file. Since
5419 the data is written in the native binary format it may not be
5420 portable between different architectures.
5422 -- Function: int gsl_vector_fread (FILE * STREAM, gsl_vector * V)
5423 This function reads into the vector V from the open stream STREAM
5424 in binary format. The vector V must be preallocated with the
5425 correct length since the function uses the size of V to determine
5426 how many bytes to read. The return value is 0 for success and
5427 `GSL_EFAILED' if there was a problem reading from the file. The
5428 data is assumed to have been written in the native binary format
5429 on the same architecture.
5431 -- Function: int gsl_vector_fprintf (FILE * STREAM, const gsl_vector *
5432 V, const char * FORMAT)
5433 This function writes the elements of the vector V line-by-line to
5434 the stream STREAM using the format specifier FORMAT, which should
5435 be one of the `%g', `%e' or `%f' formats for floating point
5436 numbers and `%d' for integers. The function returns 0 for success
5437 and `GSL_EFAILED' if there was a problem writing to the file.
5439 -- Function: int gsl_vector_fscanf (FILE * STREAM, gsl_vector * V)
5440 This function reads formatted data from the stream STREAM into the
5441 vector V. The vector V must be preallocated with the correct
5442 length since the function uses the size of V to determine how many
5443 numbers to read. The function returns 0 for success and
5444 `GSL_EFAILED' if there was a problem reading from the file.
5447 File: gsl-ref.info, Node: Vector views, Next: Copying vectors, Prev: Reading and writing vectors, Up: Vectors
5452 In addition to creating vectors from slices of blocks it is also
5453 possible to slice vectors and create vector views. For example, a
5454 subvector of another vector can be described with a view, or two views
5455 can be made which provide access to the even and odd elements of a
5458 A vector view is a temporary object, stored on the stack, which can
5459 be used to operate on a subset of vector elements. Vector views can be
5460 defined for both constant and non-constant vectors, using separate types
5461 that preserve constness. A vector view has the type `gsl_vector_view'
5462 and a constant vector view has the type `gsl_vector_const_view'. In
5463 both cases the elements of the view can be accessed as a `gsl_vector'
5464 using the `vector' component of the view object. A pointer to a vector
5465 of type `gsl_vector *' or `const gsl_vector *' can be obtained by
5466 taking the address of this component with the `&' operator.
5468 When using this pointer it is important to ensure that the view
5469 itself remains in scope--the simplest way to do so is by always writing
5470 the pointer as `&'VIEW`.vector', and never storing this value in
5473 -- Function: gsl_vector_view gsl_vector_subvector (gsl_vector * V,
5474 size_t OFFSET, size_t N)
5475 -- Function: gsl_vector_const_view gsl_vector_const_subvector (const
5476 gsl_vector * V, size_t OFFSET, size_t N)
5477 These functions return a vector view of a subvector of another
5478 vector V. The start of the new vector is offset by OFFSET elements
5479 from the start of the original vector. The new vector has N
5480 elements. Mathematically, the I-th element of the new vector V'
5483 v'(i) = v->data[(offset + i)*v->stride]
5485 where the index I runs from 0 to `n-1'.
5487 The `data' pointer of the returned vector struct is set to null if
5488 the combined parameters (OFFSET,N) overrun the end of the original
5491 The new vector is only a view of the block underlying the original
5492 vector, V. The block containing the elements of V is not owned by
5493 the new vector. When the view goes out of scope the original
5494 vector V and its block will continue to exist. The original
5495 memory can only be deallocated by freeing the original vector. Of
5496 course, the original vector should not be deallocated while the
5497 view is still in use.
5499 The function `gsl_vector_const_subvector' is equivalent to
5500 `gsl_vector_subvector' but can be used for vectors which are
5503 -- Function: gsl_vector_view gsl_vector_subvector_with_stride
5504 (gsl_vector * V, size_t OFFSET, size_t STRIDE, size_t N)
5505 -- Function: gsl_vector_const_view
5506 gsl_vector_const_subvector_with_stride (const gsl_vector * V, size_t
5507 OFFSET, size_t STRIDE, size_t N)
5508 These functions return a vector view of a subvector of another
5509 vector V with an additional stride argument. The subvector is
5510 formed in the same way as for `gsl_vector_subvector' but the new
5511 vector has N elements with a step-size of STRIDE from one element
5512 to the next in the original vector. Mathematically, the I-th
5513 element of the new vector V' is given by,
5515 v'(i) = v->data[(offset + i*stride)*v->stride]
5517 where the index I runs from 0 to `n-1'.
5519 Note that subvector views give direct access to the underlying
5520 elements of the original vector. For example, the following code
5521 will zero the even elements of the vector `v' of length `n', while
5522 leaving the odd elements untouched,
5524 gsl_vector_view v_even
5525 = gsl_vector_subvector_with_stride (v, 0, 2, n/2);
5526 gsl_vector_set_zero (&v_even.vector);
5528 A vector view can be passed to any subroutine which takes a vector
5529 argument just as a directly allocated vector would be, using
5530 `&'VIEW`.vector'. For example, the following code computes the
5531 norm of the odd elements of `v' using the BLAS routine DNRM2,
5533 gsl_vector_view v_odd
5534 = gsl_vector_subvector_with_stride (v, 1, 2, n/2);
5535 double r = gsl_blas_dnrm2 (&v_odd.vector);
5537 The function `gsl_vector_const_subvector_with_stride' is equivalent
5538 to `gsl_vector_subvector_with_stride' but can be used for vectors
5539 which are declared `const'.
5541 -- Function: gsl_vector_view gsl_vector_complex_real
5542 (gsl_vector_complex * V)
5543 -- Function: gsl_vector_const_view gsl_vector_complex_const_real
5544 (const gsl_vector_complex * V)
5545 These functions return a vector view of the real parts of the
5548 The function `gsl_vector_complex_const_real' is equivalent to
5549 `gsl_vector_complex_real' but can be used for vectors which are
5552 -- Function: gsl_vector_view gsl_vector_complex_imag
5553 (gsl_vector_complex * V)
5554 -- Function: gsl_vector_const_view gsl_vector_complex_const_imag
5555 (const gsl_vector_complex * V)
5556 These functions return a vector view of the imaginary parts of the
5559 The function `gsl_vector_complex_const_imag' is equivalent to
5560 `gsl_vector_complex_imag' but can be used for vectors which are
5563 -- Function: gsl_vector_view gsl_vector_view_array (double * BASE,
5565 -- Function: gsl_vector_const_view gsl_vector_const_view_array (const
5566 double * BASE, size_t N)
5567 These functions return a vector view of an array. The start of
5568 the new vector is given by BASE and has N elements.
5569 Mathematically, the I-th element of the new vector V' is given by,
5573 where the index I runs from 0 to `n-1'.
5575 The array containing the elements of V is not owned by the new
5576 vector view. When the view goes out of scope the original array
5577 will continue to exist. The original memory can only be
5578 deallocated by freeing the original pointer BASE. Of course, the
5579 original array should not be deallocated while the view is still
5582 The function `gsl_vector_const_view_array' is equivalent to
5583 `gsl_vector_view_array' but can be used for arrays which are
5586 -- Function: gsl_vector_view gsl_vector_view_array_with_stride (double
5587 * BASE, size_t STRIDE, size_t N)
5588 -- Function: gsl_vector_const_view
5589 gsl_vector_const_view_array_with_stride (const double * BASE, size_t
5591 These functions return a vector view of an array BASE with an
5592 additional stride argument. The subvector is formed in the same
5593 way as for `gsl_vector_view_array' but the new vector has N
5594 elements with a step-size of STRIDE from one element to the next
5595 in the original array. Mathematically, the I-th element of the new
5596 vector V' is given by,
5598 v'(i) = base[i*stride]
5600 where the index I runs from 0 to `n-1'.
5602 Note that the view gives direct access to the underlying elements
5603 of the original array. A vector view can be passed to any
5604 subroutine which takes a vector argument just as a directly
5605 allocated vector would be, using `&'VIEW`.vector'.
5607 The function `gsl_vector_const_view_array_with_stride' is
5608 equivalent to `gsl_vector_view_array_with_stride' but can be used
5609 for arrays which are declared `const'.
5612 File: gsl-ref.info, Node: Copying vectors, Next: Exchanging elements, Prev: Vector views, Up: Vectors
5614 8.3.6 Copying vectors
5615 ---------------------
5617 Common operations on vectors such as addition and multiplication are
5618 available in the BLAS part of the library (*note BLAS Support::).
5619 However, it is useful to have a small number of utility functions which
5620 do not require the full BLAS code. The following functions fall into
5623 -- Function: int gsl_vector_memcpy (gsl_vector * DEST, const
5625 This function copies the elements of the vector SRC into the
5626 vector DEST. The two vectors must have the same length.
5628 -- Function: int gsl_vector_swap (gsl_vector * V, gsl_vector * W)
5629 This function exchanges the elements of the vectors V and W by
5630 copying. The two vectors must have the same length.
5633 File: gsl-ref.info, Node: Exchanging elements, Next: Vector operations, Prev: Copying vectors, Up: Vectors
5635 8.3.7 Exchanging elements
5636 -------------------------
5638 The following function can be used to exchange, or permute, the elements
5641 -- Function: int gsl_vector_swap_elements (gsl_vector * V, size_t I,
5643 This function exchanges the I-th and J-th elements of the vector V
5646 -- Function: int gsl_vector_reverse (gsl_vector * V)
5647 This function reverses the order of the elements of the vector V.
5650 File: gsl-ref.info, Node: Vector operations, Next: Finding maximum and minimum elements of vectors, Prev: Exchanging elements, Up: Vectors
5652 8.3.8 Vector operations
5653 -----------------------
5655 The following operations are only defined for real vectors.
5657 -- Function: int gsl_vector_add (gsl_vector * A, const gsl_vector * B)
5658 This function adds the elements of vector B to the elements of
5659 vector A, a'_i = a_i + b_i. The two vectors must have the same
5662 -- Function: int gsl_vector_sub (gsl_vector * A, const gsl_vector * B)
5663 This function subtracts the elements of vector B from the elements
5664 of vector A, a'_i = a_i - b_i. The two vectors must have the same
5667 -- Function: int gsl_vector_mul (gsl_vector * A, const gsl_vector * B)
5668 This function multiplies the elements of vector A by the elements
5669 of vector B, a'_i = a_i * b_i. The two vectors must have the same
5672 -- Function: int gsl_vector_div (gsl_vector * A, const gsl_vector * B)
5673 This function divides the elements of vector A by the elements of
5674 vector B, a'_i = a_i / b_i. The two vectors must have the same
5677 -- Function: int gsl_vector_scale (gsl_vector * A, const double X)
5678 This function multiplies the elements of vector A by the constant
5679 factor X, a'_i = x a_i.
5681 -- Function: int gsl_vector_add_constant (gsl_vector * A, const double
5683 This function adds the constant value X to the elements of the
5684 vector A, a'_i = a_i + x.
5687 File: gsl-ref.info, Node: Finding maximum and minimum elements of vectors, Next: Vector properties, Prev: Vector operations, Up: Vectors
5689 8.3.9 Finding maximum and minimum elements of vectors
5690 -----------------------------------------------------
5692 -- Function: double gsl_vector_max (const gsl_vector * V)
5693 This function returns the maximum value in the vector V.
5695 -- Function: double gsl_vector_min (const gsl_vector * V)
5696 This function returns the minimum value in the vector V.
5698 -- Function: void gsl_vector_minmax (const gsl_vector * V, double *
5699 MIN_OUT, double * MAX_OUT)
5700 This function returns the minimum and maximum values in the vector
5701 V, storing them in MIN_OUT and MAX_OUT.
5703 -- Function: size_t gsl_vector_max_index (const gsl_vector * V)
5704 This function returns the index of the maximum value in the vector
5705 V. When there are several equal maximum elements then the lowest
5708 -- Function: size_t gsl_vector_min_index (const gsl_vector * V)
5709 This function returns the index of the minimum value in the vector
5710 V. When there are several equal minimum elements then the lowest
5713 -- Function: void gsl_vector_minmax_index (const gsl_vector * V,
5714 size_t * IMIN, size_t * IMAX)
5715 This function returns the indices of the minimum and maximum
5716 values in the vector V, storing them in IMIN and IMAX. When there
5717 are several equal minimum or maximum elements then the lowest
5718 indices are returned.
5721 File: gsl-ref.info, Node: Vector properties, Next: Example programs for vectors, Prev: Finding maximum and minimum elements of vectors, Up: Vectors
5723 8.3.10 Vector properties
5724 ------------------------
5726 -- Function: int gsl_vector_isnull (const gsl_vector * V)
5727 -- Function: int gsl_vector_ispos (const gsl_vector * V)
5728 -- Function: int gsl_vector_isneg (const gsl_vector * V)
5729 -- Function: int gsl_vector_isnonneg (const gsl_vector * V)
5730 These functions return 1 if all the elements of the vector V are
5731 zero, strictly positive, strictly negative, or non-negative
5732 respectively, and 0 otherwise.
5735 File: gsl-ref.info, Node: Example programs for vectors, Prev: Vector properties, Up: Vectors
5737 8.3.11 Example programs for vectors
5738 -----------------------------------
5740 This program shows how to allocate, initialize and read from a vector
5741 using the functions `gsl_vector_alloc', `gsl_vector_set' and
5745 #include <gsl/gsl_vector.h>
5751 gsl_vector * v = gsl_vector_alloc (3);
5753 for (i = 0; i < 3; i++)
5755 gsl_vector_set (v, i, 1.23 + i);
5758 for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */
5760 printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
5763 gsl_vector_free (v);
5767 Here is the output from the program. The final loop attempts to read
5768 outside the range of the vector `v', and the error is trapped by the
5769 range-checking code in `gsl_vector_get'.
5775 gsl: vector_source.c:12: ERROR: index out of range
5776 Default GSL error handler invoked.
5777 Aborted (core dumped)
5779 The next program shows how to write a vector to a file.
5782 #include <gsl/gsl_vector.h>
5788 gsl_vector * v = gsl_vector_alloc (100);
5790 for (i = 0; i < 100; i++)
5792 gsl_vector_set (v, i, 1.23 + i);
5796 FILE * f = fopen ("test.dat", "w");
5797 gsl_vector_fprintf (f, v, "%.5g");
5801 gsl_vector_free (v);
5805 After running this program the file `test.dat' should contain the
5806 elements of `v', written using the format specifier `%.5g'. The vector
5807 could then be read back in using the function `gsl_vector_fscanf (f,
5811 #include <gsl/gsl_vector.h>
5817 gsl_vector * v = gsl_vector_alloc (10);
5820 FILE * f = fopen ("test.dat", "r");
5821 gsl_vector_fscanf (f, v);
5825 for (i = 0; i < 10; i++)
5827 printf ("%g\n", gsl_vector_get(v, i));
5830 gsl_vector_free (v);
5835 File: gsl-ref.info, Node: Matrices, Next: Vector and Matrix References and Further Reading, Prev: Vectors, Up: Vectors and Matrices
5840 Matrices are defined by a `gsl_matrix' structure which describes a
5841 generalized slice of a block. Like a vector it represents a set of
5842 elements in an area of memory, but uses two indices instead of one.
5844 The `gsl_matrix' structure contains six components, the two
5845 dimensions of the matrix, a physical dimension, a pointer to the memory
5846 where the elements of the matrix are stored, DATA, a pointer to the
5847 block owned by the matrix BLOCK, if any, and an ownership flag, OWNER.
5848 The physical dimension determines the memory layout and can differ from
5849 the matrix dimension to allow the use of submatrices. The `gsl_matrix'
5850 structure is very simple and looks like this,
5862 Matrices are stored in row-major order, meaning that each row of
5863 elements forms a contiguous block in memory. This is the standard
5864 "C-language ordering" of two-dimensional arrays. Note that FORTRAN
5865 stores arrays in column-major order. The number of rows is SIZE1. The
5866 range of valid row indices runs from 0 to `size1-1'. Similarly SIZE2
5867 is the number of columns. The range of valid column indices runs from
5868 0 to `size2-1'. The physical row dimension TDA, or "trailing
5869 dimension", specifies the size of a row of the matrix as laid out in
5872 For example, in the following matrix SIZE1 is 3, SIZE2 is 4, and TDA
5873 is 8. The physical memory layout of the matrix begins in the top left
5874 hand-corner and proceeds from left to right along each row in turn.
5876 00 01 02 03 XX XX XX XX
5877 10 11 12 13 XX XX XX XX
5878 20 21 22 23 XX XX XX XX
5880 Each unused memory location is represented by "`XX'". The pointer DATA
5881 gives the location of the first element of the matrix in memory. The
5882 pointer BLOCK stores the location of the memory block in which the
5883 elements of the matrix are located (if any). If the matrix owns this
5884 block then the OWNER field is set to one and the block will be
5885 deallocated when the matrix is freed. If the matrix is only a slice of
5886 a block owned by another object then the OWNER field is zero and any
5887 underlying block will not be freed.
5889 The functions for allocating and accessing matrices are defined in
5894 * Matrix allocation::
5895 * Accessing matrix elements::
5896 * Initializing matrix elements::
5897 * Reading and writing matrices::
5899 * Creating row and column views::
5900 * Copying matrices::
5901 * Copying rows and columns::
5902 * Exchanging rows and columns::
5903 * Matrix operations::
5904 * Finding maximum and minimum elements of matrices::
5905 * Matrix properties::
5906 * Example programs for matrices::
5909 File: gsl-ref.info, Node: Matrix allocation, Next: Accessing matrix elements, Up: Matrices
5911 8.4.1 Matrix allocation
5912 -----------------------
5914 The functions for allocating memory to a matrix follow the style of
5915 `malloc' and `free'. They also perform their own error checking. If
5916 there is insufficient memory available to allocate a matrix then the
5917 functions call the GSL error handler (with an error number of
5918 `GSL_ENOMEM') in addition to returning a null pointer. Thus if you use
5919 the library error handler to abort your program then it isn't necessary
5920 to check every `alloc'.
5922 -- Function: gsl_matrix * gsl_matrix_alloc (size_t N1, size_t N2)
5923 This function creates a matrix of size N1 rows by N2 columns,
5924 returning a pointer to a newly initialized matrix struct. A new
5925 block is allocated for the elements of the matrix, and stored in
5926 the BLOCK component of the matrix struct. The block is "owned" by
5927 the matrix, and will be deallocated when the matrix is deallocated.
5929 -- Function: gsl_matrix * gsl_matrix_calloc (size_t N1, size_t N2)
5930 This function allocates memory for a matrix of size N1 rows by N2
5931 columns and initializes all the elements of the matrix to zero.
5933 -- Function: void gsl_matrix_free (gsl_matrix * M)
5934 This function frees a previously allocated matrix M. If the
5935 matrix was created using `gsl_matrix_alloc' then the block
5936 underlying the matrix will also be deallocated. If the matrix has
5937 been created from another object then the memory is still owned by
5938 that object and will not be deallocated. The matrix M must be a
5939 valid matrix object (a null pointer is not allowed).
5942 File: gsl-ref.info, Node: Accessing matrix elements, Next: Initializing matrix elements, Prev: Matrix allocation, Up: Matrices
5944 8.4.2 Accessing matrix elements
5945 -------------------------------
5947 The functions for accessing the elements of a matrix use the same range
5948 checking system as vectors. You can turn off range checking by
5949 recompiling your program with the preprocessor definition
5950 `GSL_RANGE_CHECK_OFF'.
5952 The elements of the matrix are stored in "C-order", where the second
5953 index moves continuously through memory. More precisely, the element
5954 accessed by the function `gsl_matrix_get(m,i,j)' and
5955 `gsl_matrix_set(m,i,j,x)' is
5957 m->data[i * m->tda + j]
5959 where TDA is the physical row-length of the matrix.
5961 -- Function: double gsl_matrix_get (const gsl_matrix * M, size_t I,
5963 This function returns the (i,j)-th element of a matrix M. If I or
5964 J lie outside the allowed range of 0 to N1-1 and 0 to N2-1 then
5965 the error handler is invoked and 0 is returned. An inline version
5966 of this function is used when `HAVE_INLINE' is defined.
5968 -- Function: void gsl_matrix_set (gsl_matrix * M, size_t I, size_t J,
5970 This function sets the value of the (i,j)-th element of a matrix M
5971 to X. If I or J lies outside the allowed range of 0 to N1-1 and 0
5972 to N2-1 then the error handler is invoked. An inline version of
5973 this function is used when `HAVE_INLINE' is defined.
5975 -- Function: double * gsl_matrix_ptr (gsl_matrix * M, size_t I, size_t
5977 -- Function: const double * gsl_matrix_const_ptr (const gsl_matrix *
5978 M, size_t I, size_t J)
5979 These functions return a pointer to the (i,j)-th element of a
5980 matrix M. If I or J lie outside the allowed range of 0 to N1-1
5981 and 0 to N2-1 then the error handler is invoked and a null pointer
5982 is returned. Inline versions of these functions are used when
5983 `HAVE_INLINE' is defined.
5986 File: gsl-ref.info, Node: Initializing matrix elements, Next: Reading and writing matrices, Prev: Accessing matrix elements, Up: Matrices
5988 8.4.3 Initializing matrix elements
5989 ----------------------------------
5991 -- Function: void gsl_matrix_set_all (gsl_matrix * M, double X)
5992 This function sets all the elements of the matrix M to the value X.
5994 -- Function: void gsl_matrix_set_zero (gsl_matrix * M)
5995 This function sets all the elements of the matrix M to zero.
5997 -- Function: void gsl_matrix_set_identity (gsl_matrix * M)
5998 This function sets the elements of the matrix M to the
5999 corresponding elements of the identity matrix, m(i,j) =
6000 \delta(i,j), i.e. a unit diagonal with all off-diagonal elements
6001 zero. This applies to both square and rectangular matrices.
6004 File: gsl-ref.info, Node: Reading and writing matrices, Next: Matrix views, Prev: Initializing matrix elements, Up: Matrices
6006 8.4.4 Reading and writing matrices
6007 ----------------------------------
6009 The library provides functions for reading and writing matrices to a
6010 file as binary data or formatted text.
6012 -- Function: int gsl_matrix_fwrite (FILE * STREAM, const gsl_matrix *
6014 This function writes the elements of the matrix M to the stream
6015 STREAM in binary format. The return value is 0 for success and
6016 `GSL_EFAILED' if there was a problem writing to the file. Since
6017 the data is written in the native binary format it may not be
6018 portable between different architectures.
6020 -- Function: int gsl_matrix_fread (FILE * STREAM, gsl_matrix * M)
6021 This function reads into the matrix M from the open stream STREAM
6022 in binary format. The matrix M must be preallocated with the
6023 correct dimensions since the function uses the size of M to
6024 determine how many bytes to read. The return value is 0 for
6025 success and `GSL_EFAILED' if there was a problem reading from the
6026 file. The data is assumed to have been written in the native
6027 binary format on the same architecture.
6029 -- Function: int gsl_matrix_fprintf (FILE * STREAM, const gsl_matrix *
6030 M, const char * FORMAT)
6031 This function writes the elements of the matrix M line-by-line to
6032 the stream STREAM using the format specifier FORMAT, which should
6033 be one of the `%g', `%e' or `%f' formats for floating point
6034 numbers and `%d' for integers. The function returns 0 for success
6035 and `GSL_EFAILED' if there was a problem writing to the file.
6037 -- Function: int gsl_matrix_fscanf (FILE * STREAM, gsl_matrix * M)
6038 This function reads formatted data from the stream STREAM into the
6039 matrix M. The matrix M must be preallocated with the correct
6040 dimensions since the function uses the size of M to determine how
6041 many numbers to read. The function returns 0 for success and
6042 `GSL_EFAILED' if there was a problem reading from the file.
6045 File: gsl-ref.info, Node: Matrix views, Next: Creating row and column views, Prev: Reading and writing matrices, Up: Matrices
6050 A matrix view is a temporary object, stored on the stack, which can be
6051 used to operate on a subset of matrix elements. Matrix views can be
6052 defined for both constant and non-constant matrices using separate types
6053 that preserve constness. A matrix view has the type `gsl_matrix_view'
6054 and a constant matrix view has the type `gsl_matrix_const_view'. In
6055 both cases the elements of the view can by accessed using the `matrix'
6056 component of the view object. A pointer `gsl_matrix *' or `const
6057 gsl_matrix *' can be obtained by taking the address of the `matrix'
6058 component with the `&' operator. In addition to matrix views it is
6059 also possible to create vector views of a matrix, such as row or column
6062 -- Function: gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * M,
6063 size_t K1, size_t K2, size_t N1, size_t N2)
6064 -- Function: gsl_matrix_const_view gsl_matrix_const_submatrix (const
6065 gsl_matrix * M, size_t K1, size_t K2, size_t N1, size_t N2)
6066 These functions return a matrix view of a submatrix of the matrix
6067 M. The upper-left element of the submatrix is the element (K1,K2)
6068 of the original matrix. The submatrix has N1 rows and N2 columns.
6069 The physical number of columns in memory given by TDA is
6070 unchanged. Mathematically, the (i,j)-th element of the new matrix
6073 m'(i,j) = m->data[(k1*m->tda + k2) + i*m->tda + j]
6075 where the index I runs from 0 to `n1-1' and the index J runs from
6078 The `data' pointer of the returned matrix struct is set to null if
6079 the combined parameters (I,J,N1,N2,TDA) overrun the ends of the
6082 The new matrix view is only a view of the block underlying the
6083 existing matrix, M. The block containing the elements of M is not
6084 owned by the new matrix view. When the view goes out of scope the
6085 original matrix M and its block will continue to exist. The
6086 original memory can only be deallocated by freeing the original
6087 matrix. Of course, the original matrix should not be deallocated
6088 while the view is still in use.
6090 The function `gsl_matrix_const_submatrix' is equivalent to
6091 `gsl_matrix_submatrix' but can be used for matrices which are
6094 -- Function: gsl_matrix_view gsl_matrix_view_array (double * BASE,
6095 size_t N1, size_t N2)
6096 -- Function: gsl_matrix_const_view gsl_matrix_const_view_array (const
6097 double * BASE, size_t N1, size_t N2)
6098 These functions return a matrix view of the array BASE. The
6099 matrix has N1 rows and N2 columns. The physical number of columns
6100 in memory is also given by N2. Mathematically, the (i,j)-th
6101 element of the new matrix is given by,
6103 m'(i,j) = base[i*n2 + j]
6105 where the index I runs from 0 to `n1-1' and the index J runs from
6108 The new matrix is only a view of the array BASE. When the view
6109 goes out of scope the original array BASE will continue to exist.
6110 The original memory can only be deallocated by freeing the original
6111 array. Of course, the original array should not be deallocated
6112 while the view is still in use.
6114 The function `gsl_matrix_const_view_array' is equivalent to
6115 `gsl_matrix_view_array' but can be used for matrices which are
6118 -- Function: gsl_matrix_view gsl_matrix_view_array_with_tda (double *
6119 BASE, size_t N1, size_t N2, size_t TDA)
6120 -- Function: gsl_matrix_const_view
6121 gsl_matrix_const_view_array_with_tda (const double * BASE,
6122 size_t N1, size_t N2, size_t TDA)
6123 These functions return a matrix view of the array BASE with a
6124 physical number of columns TDA which may differ from the
6125 corresponding dimension of the matrix. The matrix has N1 rows and
6126 N2 columns, and the physical number of columns in memory is given
6127 by TDA. Mathematically, the (i,j)-th element of the new matrix is
6130 m'(i,j) = base[i*tda + j]
6132 where the index I runs from 0 to `n1-1' and the index J runs from
6135 The new matrix is only a view of the array BASE. When the view
6136 goes out of scope the original array BASE will continue to exist.
6137 The original memory can only be deallocated by freeing the original
6138 array. Of course, the original array should not be deallocated
6139 while the view is still in use.
6141 The function `gsl_matrix_const_view_array_with_tda' is equivalent
6142 to `gsl_matrix_view_array_with_tda' but can be used for matrices
6143 which are declared `const'.
6145 -- Function: gsl_matrix_view gsl_matrix_view_vector (gsl_vector * V,
6146 size_t N1, size_t N2)
6147 -- Function: gsl_matrix_const_view gsl_matrix_const_view_vector (const
6148 gsl_vector * V, size_t N1, size_t N2)
6149 These functions return a matrix view of the vector V. The matrix
6150 has N1 rows and N2 columns. The vector must have unit stride. The
6151 physical number of columns in memory is also given by N2.
6152 Mathematically, the (i,j)-th element of the new matrix is given by,
6154 m'(i,j) = v->data[i*n2 + j]
6156 where the index I runs from 0 to `n1-1' and the index J runs from
6159 The new matrix is only a view of the vector V. When the view goes
6160 out of scope the original vector V will continue to exist. The
6161 original memory can only be deallocated by freeing the original
6162 vector. Of course, the original vector should not be deallocated
6163 while the view is still in use.
6165 The function `gsl_matrix_const_view_vector' is equivalent to
6166 `gsl_matrix_view_vector' but can be used for matrices which are
6169 -- Function: gsl_matrix_view gsl_matrix_view_vector_with_tda
6170 (gsl_vector * V, size_t N1, size_t N2, size_t TDA)
6171 -- Function: gsl_matrix_const_view
6172 gsl_matrix_const_view_vector_with_tda (const gsl_vector * V, size_t N1,
6173 size_t N2, size_t TDA)
6174 These functions return a matrix view of the vector V with a
6175 physical number of columns TDA which may differ from the
6176 corresponding matrix dimension. The vector must have unit stride.
6177 The matrix has N1 rows and N2 columns, and the physical number of
6178 columns in memory is given by TDA. Mathematically, the (i,j)-th
6179 element of the new matrix is given by,
6181 m'(i,j) = v->data[i*tda + j]
6183 where the index I runs from 0 to `n1-1' and the index J runs from
6186 The new matrix is only a view of the vector V. When the view goes
6187 out of scope the original vector V will continue to exist. The
6188 original memory can only be deallocated by freeing the original
6189 vector. Of course, the original vector should not be deallocated
6190 while the view is still in use.
6192 The function `gsl_matrix_const_view_vector_with_tda' is equivalent
6193 to `gsl_matrix_view_vector_with_tda' but can be used for matrices
6194 which are declared `const'.
6197 File: gsl-ref.info, Node: Creating row and column views, Next: Copying matrices, Prev: Matrix views, Up: Matrices
6199 8.4.6 Creating row and column views
6200 -----------------------------------
6202 In general there are two ways to access an object, by reference or by
6203 copying. The functions described in this section create vector views
6204 which allow access to a row or column of a matrix by reference.
6205 Modifying elements of the view is equivalent to modifying the matrix,
6206 since both the vector view and the matrix point to the same memory
6209 -- Function: gsl_vector_view gsl_matrix_row (gsl_matrix * M, size_t I)
6210 -- Function: gsl_vector_const_view gsl_matrix_const_row (const
6211 gsl_matrix * M, size_t I)
6212 These functions return a vector view of the I-th row of the matrix
6213 M. The `data' pointer of the new vector is set to null if I is
6216 The function `gsl_vector_const_row' is equivalent to
6217 `gsl_matrix_row' but can be used for matrices which are declared
6220 -- Function: gsl_vector_view gsl_matrix_column (gsl_matrix * M, size_t
6222 -- Function: gsl_vector_const_view gsl_matrix_const_column (const
6223 gsl_matrix * M, size_t J)
6224 These functions return a vector view of the J-th column of the
6225 matrix M. The `data' pointer of the new vector is set to null if
6228 The function `gsl_vector_const_column' is equivalent to
6229 `gsl_matrix_column' but can be used for matrices which are declared
6232 -- Function: gsl_vector_view gsl_matrix_subrow (gsl_matrix * M, size_t
6233 I, size_t OFFSET, size_t N)
6234 -- Function: gsl_vector_const_view gsl_matrix_const_subrow (const
6235 gsl_matrix * M, size_t I, size_t OFFSET, size_t N)
6236 These functions return a vector view of the I-th row of the matrix
6237 M beginning at OFFSET elements past the first column and
6238 containing N elements. The `data' pointer of the new vector is set
6239 to null if I, OFFSET, or N are out of range.
6241 The function `gsl_vector_const_subrow' is equivalent to
6242 `gsl_matrix_subrow' but can be used for matrices which are declared
6245 -- Function: gsl_vector_view gsl_matrix_subcolumn (gsl_matrix * M,
6246 size_t J, size_t OFFSET, size_t N)
6247 -- Function: gsl_vector_const_view gsl_matrix_const_subcolumn (const
6248 gsl_matrix * M, size_t J, size_t OFFSET, size_t N)
6249 These functions return a vector view of the J-th column of the
6250 matrix M beginning at OFFSET elements past the first row and
6251 containing N elements. The `data' pointer of the new vector is set
6252 to null if J, OFFSET, or N are out of range.
6254 The function `gsl_vector_const_subcolumn' is equivalent to
6255 `gsl_matrix_subcolumn' but can be used for matrices which are
6258 -- Function: gsl_vector_view gsl_matrix_diagonal (gsl_matrix * M)
6259 -- Function: gsl_vector_const_view gsl_matrix_const_diagonal (const
6261 These functions returns a vector view of the diagonal of the matrix
6262 M. The matrix M is not required to be square. For a rectangular
6263 matrix the length of the diagonal is the same as the smaller
6264 dimension of the matrix.
6266 The function `gsl_matrix_const_diagonal' is equivalent to
6267 `gsl_matrix_diagonal' but can be used for matrices which are
6270 -- Function: gsl_vector_view gsl_matrix_subdiagonal (gsl_matrix * M,
6272 -- Function: gsl_vector_const_view gsl_matrix_const_subdiagonal (const
6273 gsl_matrix * M, size_t K)
6274 These functions return a vector view of the K-th subdiagonal of
6275 the matrix M. The matrix M is not required to be square. The
6276 diagonal of the matrix corresponds to k = 0.
6278 The function `gsl_matrix_const_subdiagonal' is equivalent to
6279 `gsl_matrix_subdiagonal' but can be used for matrices which are
6282 -- Function: gsl_vector_view gsl_matrix_superdiagonal (gsl_matrix * M,
6284 -- Function: gsl_vector_const_view gsl_matrix_const_superdiagonal
6285 (const gsl_matrix * M, size_t K)
6286 These functions return a vector view of the K-th superdiagonal of
6287 the matrix M. The matrix M is not required to be square. The
6288 diagonal of the matrix corresponds to k = 0.
6290 The function `gsl_matrix_const_superdiagonal' is equivalent to
6291 `gsl_matrix_superdiagonal' but can be used for matrices which are
6295 File: gsl-ref.info, Node: Copying matrices, Next: Copying rows and columns, Prev: Creating row and column views, Up: Matrices
6297 8.4.7 Copying matrices
6298 ----------------------
6300 -- Function: int gsl_matrix_memcpy (gsl_matrix * DEST, const
6302 This function copies the elements of the matrix SRC into the
6303 matrix DEST. The two matrices must have the same size.
6305 -- Function: int gsl_matrix_swap (gsl_matrix * M1, gsl_matrix * M2)
6306 This function exchanges the elements of the matrices M1 and M2 by
6307 copying. The two matrices must have the same size.
6310 File: gsl-ref.info, Node: Copying rows and columns, Next: Exchanging rows and columns, Prev: Copying matrices, Up: Matrices
6312 8.4.8 Copying rows and columns
6313 ------------------------------
6315 The functions described in this section copy a row or column of a matrix
6316 into a vector. This allows the elements of the vector and the matrix to
6317 be modified independently. Note that if the matrix and the vector point
6318 to overlapping regions of memory then the result will be undefined. The
6319 same effect can be achieved with more generality using
6320 `gsl_vector_memcpy' with vector views of rows and columns.
6322 -- Function: int gsl_matrix_get_row (gsl_vector * V, const gsl_matrix
6324 This function copies the elements of the I-th row of the matrix M
6325 into the vector V. The length of the vector must be the same as
6326 the length of the row.
6328 -- Function: int gsl_matrix_get_col (gsl_vector * V, const gsl_matrix
6330 This function copies the elements of the J-th column of the matrix
6331 M into the vector V. The length of the vector must be the same as
6332 the length of the column.
6334 -- Function: int gsl_matrix_set_row (gsl_matrix * M, size_t I, const
6336 This function copies the elements of the vector V into the I-th
6337 row of the matrix M. The length of the vector must be the same as
6338 the length of the row.
6340 -- Function: int gsl_matrix_set_col (gsl_matrix * M, size_t J, const
6342 This function copies the elements of the vector V into the J-th
6343 column of the matrix M. The length of the vector must be the same
6344 as the length of the column.
6347 File: gsl-ref.info, Node: Exchanging rows and columns, Next: Matrix operations, Prev: Copying rows and columns, Up: Matrices
6349 8.4.9 Exchanging rows and columns
6350 ---------------------------------
6352 The following functions can be used to exchange the rows and columns of
6355 -- Function: int gsl_matrix_swap_rows (gsl_matrix * M, size_t I,
6357 This function exchanges the I-th and J-th rows of the matrix M
6360 -- Function: int gsl_matrix_swap_columns (gsl_matrix * M, size_t I,
6362 This function exchanges the I-th and J-th columns of the matrix M
6365 -- Function: int gsl_matrix_swap_rowcol (gsl_matrix * M, size_t I,
6367 This function exchanges the I-th row and J-th column of the matrix
6368 M in-place. The matrix must be square for this operation to be
6371 -- Function: int gsl_matrix_transpose_memcpy (gsl_matrix * DEST, const
6373 This function makes the matrix DEST the transpose of the matrix
6374 SRC by copying the elements of SRC into DEST. This function works
6375 for all matrices provided that the dimensions of the matrix DEST
6376 match the transposed dimensions of the matrix SRC.
6378 -- Function: int gsl_matrix_transpose (gsl_matrix * M)
6379 This function replaces the matrix M by its transpose by copying
6380 the elements of the matrix in-place. The matrix must be square
6381 for this operation to be possible.
6384 File: gsl-ref.info, Node: Matrix operations, Next: Finding maximum and minimum elements of matrices, Prev: Exchanging rows and columns, Up: Matrices
6386 8.4.10 Matrix operations
6387 ------------------------
6389 The following operations are defined for real and complex matrices.
6391 -- Function: int gsl_matrix_add (gsl_matrix * A, const gsl_matrix * B)
6392 This function adds the elements of matrix B to the elements of
6393 matrix A, a'(i,j) = a(i,j) + b(i,j). The two matrices must have the
6396 -- Function: int gsl_matrix_sub (gsl_matrix * A, const gsl_matrix * B)
6397 This function subtracts the elements of matrix B from the elements
6398 of matrix A, a'(i,j) = a(i,j) - b(i,j). The two matrices must have
6399 the same dimensions.
6401 -- Function: int gsl_matrix_mul_elements (gsl_matrix * A, const
6403 This function multiplies the elements of matrix A by the elements
6404 of matrix B, a'(i,j) = a(i,j) * b(i,j). The two matrices must have
6405 the same dimensions.
6407 -- Function: int gsl_matrix_div_elements (gsl_matrix * A, const
6409 This function divides the elements of matrix A by the elements of
6410 matrix B, a'(i,j) = a(i,j) / b(i,j). The two matrices must have the
6413 -- Function: int gsl_matrix_scale (gsl_matrix * A, const double X)
6414 This function multiplies the elements of matrix A by the constant
6415 factor X, a'(i,j) = x a(i,j).
6417 -- Function: int gsl_matrix_add_constant (gsl_matrix * A, const double
6419 This function adds the constant value X to the elements of the
6420 matrix A, a'(i,j) = a(i,j) + x.
6423 File: gsl-ref.info, Node: Finding maximum and minimum elements of matrices, Next: Matrix properties, Prev: Matrix operations, Up: Matrices
6425 8.4.11 Finding maximum and minimum elements of matrices
6426 -------------------------------------------------------
6428 The following operations are only defined for real matrices.
6430 -- Function: double gsl_matrix_max (const gsl_matrix * M)
6431 This function returns the maximum value in the matrix M.
6433 -- Function: double gsl_matrix_min (const gsl_matrix * M)
6434 This function returns the minimum value in the matrix M.
6436 -- Function: void gsl_matrix_minmax (const gsl_matrix * M, double *
6437 MIN_OUT, double * MAX_OUT)
6438 This function returns the minimum and maximum values in the matrix
6439 M, storing them in MIN_OUT and MAX_OUT.
6441 -- Function: void gsl_matrix_max_index (const gsl_matrix * M, size_t *
6442 IMAX, size_t * JMAX)
6443 This function returns the indices of the maximum value in the
6444 matrix M, storing them in IMAX and JMAX. When there are several
6445 equal maximum elements then the first element found is returned,
6446 searching in row-major order.
6448 -- Function: void gsl_matrix_min_index (const gsl_matrix * M, size_t *
6449 IMIN, size_t * JMIN)
6450 This function returns the indices of the minimum value in the
6451 matrix M, storing them in IMIN and JMIN. When there are several
6452 equal minimum elements then the first element found is returned,
6453 searching in row-major order.
6455 -- Function: void gsl_matrix_minmax_index (const gsl_matrix * M,
6456 size_t * IMIN, size_t * JMIN, size_t * IMAX, size_t * JMAX)
6457 This function returns the indices of the minimum and maximum
6458 values in the matrix M, storing them in (IMIN,JMIN) and
6459 (IMAX,JMAX). When there are several equal minimum or maximum
6460 elements then the first elements found are returned, searching in
6464 File: gsl-ref.info, Node: Matrix properties, Next: Example programs for matrices, Prev: Finding maximum and minimum elements of matrices, Up: Matrices
6466 8.4.12 Matrix properties
6467 ------------------------
6469 -- Function: int gsl_matrix_isnull (const gsl_matrix * M)
6470 -- Function: int gsl_matrix_ispos (const gsl_matrix * M)
6471 -- Function: int gsl_matrix_isneg (const gsl_matrix * M)
6472 -- Function: int gsl_matrix_isnonneg (const gsl_matrix * M)
6473 These functions return 1 if all the elements of the matrix M are
6474 zero, strictly positive, strictly negative, or non-negative
6475 respectively, and 0 otherwise. To test whether a matrix is
6476 positive-definite, use the Cholesky decomposition (*note Cholesky
6480 File: gsl-ref.info, Node: Example programs for matrices, Prev: Matrix properties, Up: Matrices
6482 8.4.13 Example programs for matrices
6483 ------------------------------------
6485 The program below shows how to allocate, initialize and read from a
6486 matrix using the functions `gsl_matrix_alloc', `gsl_matrix_set' and
6490 #include <gsl/gsl_matrix.h>
6496 gsl_matrix * m = gsl_matrix_alloc (10, 3);
6498 for (i = 0; i < 10; i++)
6499 for (j = 0; j < 3; j++)
6500 gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
6502 for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */
6503 for (j = 0; j < 3; j++)
6504 printf ("m(%d,%d) = %g\n", i, j,
6505 gsl_matrix_get (m, i, j));
6507 gsl_matrix_free (m);
6512 Here is the output from the program. The final loop attempts to read
6513 outside the range of the matrix `m', and the error is trapped by the
6514 range-checking code in `gsl_matrix_get'.
6525 gsl: matrix_source.c:13: ERROR: first index out of range
6526 Default GSL error handler invoked.
6527 Aborted (core dumped)
6529 The next program shows how to write a matrix to a file.
6532 #include <gsl/gsl_matrix.h>
6538 gsl_matrix * m = gsl_matrix_alloc (100, 100);
6539 gsl_matrix * a = gsl_matrix_alloc (100, 100);
6541 for (i = 0; i < 100; i++)
6542 for (j = 0; j < 100; j++)
6543 gsl_matrix_set (m, i, j, 0.23 + i + j);
6546 FILE * f = fopen ("test.dat", "wb");
6547 gsl_matrix_fwrite (f, m);
6552 FILE * f = fopen ("test.dat", "rb");
6553 gsl_matrix_fread (f, a);
6557 for (i = 0; i < 100; i++)
6558 for (j = 0; j < 100; j++)
6560 double mij = gsl_matrix_get (m, i, j);
6561 double aij = gsl_matrix_get (a, i, j);
6562 if (mij != aij) k++;
6565 gsl_matrix_free (m);
6566 gsl_matrix_free (a);
6568 printf ("differences = %d (should be zero)\n", k);
6572 After running this program the file `test.dat' should contain the
6573 elements of `m', written in binary format. The matrix which is read
6574 back in using the function `gsl_matrix_fread' should be exactly equal
6575 to the original matrix.
6577 The following program demonstrates the use of vector views. The
6578 program computes the column norms of a matrix.
6582 #include <gsl/gsl_matrix.h>
6583 #include <gsl/gsl_blas.h>
6590 gsl_matrix *m = gsl_matrix_alloc (10, 10);
6592 for (i = 0; i < 10; i++)
6593 for (j = 0; j < 10; j++)
6594 gsl_matrix_set (m, i, j, sin (i) + cos (j));
6596 for (j = 0; j < 10; j++)
6598 gsl_vector_view column = gsl_matrix_column (m, j);
6601 d = gsl_blas_dnrm2 (&column.vector);
6603 printf ("matrix column %d, norm = %g\n", j, d);
6606 gsl_matrix_free (m);
6611 Here is the output of the program,
6614 matrix column 0, norm = 4.31461
6615 matrix column 1, norm = 3.1205
6616 matrix column 2, norm = 2.19316
6617 matrix column 3, norm = 3.26114
6618 matrix column 4, norm = 2.53416
6619 matrix column 5, norm = 2.57281
6620 matrix column 6, norm = 4.20469
6621 matrix column 7, norm = 3.65202
6622 matrix column 8, norm = 2.08524
6623 matrix column 9, norm = 3.07313
6625 The results can be confirmed using GNU OCTAVE,
6628 GNU Octave, version 2.0.16.92
6629 octave> m = sin(0:9)' * ones(1,10)
6630 + ones(10,1) * cos(0:9);
6631 octave> sqrt(sum(m.^2))
6633 4.3146 3.1205 2.1932 3.2611 2.5342 2.5728
6634 4.2047 3.6520 2.0852 3.0731
6637 File: gsl-ref.info, Node: Vector and Matrix References and Further Reading, Prev: Matrices, Up: Vectors and Matrices
6639 8.5 References and Further Reading
6640 ==================================
6642 The block, vector and matrix objects in GSL follow the `valarray' model
6643 of C++. A description of this model can be found in the following
6646 B. Stroustrup, `The C++ Programming Language' (3rd Ed), Section
6647 22.4 Vector Arithmetic. Addison-Wesley 1997, ISBN 0-201-88954-4.
6650 File: gsl-ref.info, Node: Permutations, Next: Combinations, Prev: Vectors and Matrices, Up: Top
6655 This chapter describes functions for creating and manipulating
6656 permutations. A permutation p is represented by an array of n integers
6657 in the range 0 to n-1, where each value p_i occurs once and only once.
6658 The application of a permutation p to a vector v yields a new vector v'
6659 where v'_i = v_{p_i}. For example, the array (0,1,3,2) represents a
6660 permutation which exchanges the last two elements of a four element
6661 vector. The corresponding identity permutation is (0,1,2,3).
6663 Note that the permutations produced by the linear algebra routines
6664 correspond to the exchange of matrix columns, and so should be
6665 considered as applying to row-vectors in the form v' = v P rather than
6666 column-vectors, when permuting the elements of a vector.
6668 The functions described in this chapter are defined in the header
6669 file `gsl_permutation.h'.
6673 * The Permutation struct::
6674 * Permutation allocation::
6675 * Accessing permutation elements::
6676 * Permutation properties::
6677 * Permutation functions::
6678 * Applying Permutations::
6679 * Reading and writing permutations::
6680 * Permutations in cyclic form::
6681 * Permutation Examples::
6682 * Permutation References and Further Reading::
6685 File: gsl-ref.info, Node: The Permutation struct, Next: Permutation allocation, Up: Permutations
6687 9.1 The Permutation struct
6688 ==========================
6690 A permutation is defined by a structure containing two components, the
6691 size of the permutation and a pointer to the permutation array. The
6692 elements of the permutation array are all of type `size_t'. The
6693 `gsl_permutation' structure looks like this,
6703 File: gsl-ref.info, Node: Permutation allocation, Next: Accessing permutation elements, Prev: The Permutation struct, Up: Permutations
6705 9.2 Permutation allocation
6706 ==========================
6708 -- Function: gsl_permutation * gsl_permutation_alloc (size_t N)
6709 This function allocates memory for a new permutation of size N.
6710 The permutation is not initialized and its elements are undefined.
6711 Use the function `gsl_permutation_calloc' if you want to create a
6712 permutation which is initialized to the identity. A null pointer is
6713 returned if insufficient memory is available to create the
6716 -- Function: gsl_permutation * gsl_permutation_calloc (size_t N)
6717 This function allocates memory for a new permutation of size N and
6718 initializes it to the identity. A null pointer is returned if
6719 insufficient memory is available to create the permutation.
6721 -- Function: void gsl_permutation_init (gsl_permutation * P)
6722 This function initializes the permutation P to the identity, i.e.
6725 -- Function: void gsl_permutation_free (gsl_permutation * P)
6726 This function frees all the memory used by the permutation P.
6728 -- Function: int gsl_permutation_memcpy (gsl_permutation * DEST, const
6729 gsl_permutation * SRC)
6730 This function copies the elements of the permutation SRC into the
6731 permutation DEST. The two permutations must have the same size.
6734 File: gsl-ref.info, Node: Accessing permutation elements, Next: Permutation properties, Prev: Permutation allocation, Up: Permutations
6736 9.3 Accessing permutation elements
6737 ==================================
6739 The following functions can be used to access and manipulate
6742 -- Function: size_t gsl_permutation_get (const gsl_permutation * P,
6744 This function returns the value of the I-th element of the
6745 permutation P. If I lies outside the allowed range of 0 to N-1
6746 then the error handler is invoked and 0 is returned. An inline
6747 version of this function is used when `HAVE_INLINE' is defined.
6749 -- Function: int gsl_permutation_swap (gsl_permutation * P, const
6750 size_t I, const size_t J)
6751 This function exchanges the I-th and J-th elements of the
6755 File: gsl-ref.info, Node: Permutation properties, Next: Permutation functions, Prev: Accessing permutation elements, Up: Permutations
6757 9.4 Permutation properties
6758 ==========================
6760 -- Function: size_t gsl_permutation_size (const gsl_permutation * P)
6761 This function returns the size of the permutation P.
6763 -- Function: size_t * gsl_permutation_data (const gsl_permutation * P)
6764 This function returns a pointer to the array of elements in the
6767 -- Function: int gsl_permutation_valid (const gsl_permutation * P)
6768 This function checks that the permutation P is valid. The N
6769 elements should contain each of the numbers 0 to N-1 once and only
6773 File: gsl-ref.info, Node: Permutation functions, Next: Applying Permutations, Prev: Permutation properties, Up: Permutations
6775 9.5 Permutation functions
6776 =========================
6778 -- Function: void gsl_permutation_reverse (gsl_permutation * P)
6779 This function reverses the elements of the permutation P.
6781 -- Function: int gsl_permutation_inverse (gsl_permutation * INV, const
6782 gsl_permutation * P)
6783 This function computes the inverse of the permutation P, storing
6786 -- Function: int gsl_permutation_next (gsl_permutation * P)
6787 This function advances the permutation P to the next permutation
6788 in lexicographic order and returns `GSL_SUCCESS'. If no further
6789 permutations are available it returns `GSL_FAILURE' and leaves P
6790 unmodified. Starting with the identity permutation and repeatedly
6791 applying this function will iterate through all possible
6792 permutations of a given order.
6794 -- Function: int gsl_permutation_prev (gsl_permutation * P)
6795 This function steps backwards from the permutation P to the
6796 previous permutation in lexicographic order, returning
6797 `GSL_SUCCESS'. If no previous permutation is available it returns
6798 `GSL_FAILURE' and leaves P unmodified.
6801 File: gsl-ref.info, Node: Applying Permutations, Next: Reading and writing permutations, Prev: Permutation functions, Up: Permutations
6803 9.6 Applying Permutations
6804 =========================
6806 -- Function: int gsl_permute (const size_t * P, double * DATA, size_t
6808 This function applies the permutation P to the array DATA of size
6809 N with stride STRIDE.
6811 -- Function: int gsl_permute_inverse (const size_t * P, double * DATA,
6812 size_t STRIDE, size_t N)
6813 This function applies the inverse of the permutation P to the
6814 array DATA of size N with stride STRIDE.
6816 -- Function: int gsl_permute_vector (const gsl_permutation * P,
6818 This function applies the permutation P to the elements of the
6819 vector V, considered as a row-vector acted on by a permutation
6820 matrix from the right, v' = v P. The j-th column of the
6821 permutation matrix P is given by the p_j-th column of the identity
6822 matrix. The permutation P and the vector V must have the same
6825 -- Function: int gsl_permute_vector_inverse (const gsl_permutation *
6827 This function applies the inverse of the permutation P to the
6828 elements of the vector V, considered as a row-vector acted on by
6829 an inverse permutation matrix from the right, v' = v P^T. Note
6830 that for permutation matrices the inverse is the same as the
6831 transpose. The j-th column of the permutation matrix P is given by
6832 the p_j-th column of the identity matrix. The permutation P and
6833 the vector V must have the same length.
6835 -- Function: int gsl_permutation_mul (gsl_permutation * P, const
6836 gsl_permutation * PA, const gsl_permutation * PB)
6837 This function combines the two permutations PA and PB into a
6838 single permutation P, where p = pa . pb. The permutation P is
6839 equivalent to applying pb first and then PA.
6842 File: gsl-ref.info, Node: Reading and writing permutations, Next: Permutations in cyclic form, Prev: Applying Permutations, Up: Permutations
6844 9.7 Reading and writing permutations
6845 ====================================
6847 The library provides functions for reading and writing permutations to a
6848 file as binary data or formatted text.
6850 -- Function: int gsl_permutation_fwrite (FILE * STREAM, const
6851 gsl_permutation * P)
6852 This function writes the elements of the permutation P to the
6853 stream STREAM in binary format. The function returns
6854 `GSL_EFAILED' if there was a problem writing to the file. Since
6855 the data is written in the native binary format it may not be
6856 portable between different architectures.
6858 -- Function: int gsl_permutation_fread (FILE * STREAM, gsl_permutation
6860 This function reads into the permutation P from the open stream
6861 STREAM in binary format. The permutation P must be preallocated
6862 with the correct length since the function uses the size of P to
6863 determine how many bytes to read. The function returns
6864 `GSL_EFAILED' if there was a problem reading from the file. The
6865 data is assumed to have been written in the native binary format
6866 on the same architecture.
6868 -- Function: int gsl_permutation_fprintf (FILE * STREAM, const
6869 gsl_permutation * P, const char * FORMAT)
6870 This function writes the elements of the permutation P
6871 line-by-line to the stream STREAM using the format specifier
6872 FORMAT, which should be suitable for a type of SIZE_T. In ISO C99
6873 the type modifier `z' represents `size_t', so `"%zu\n"' is a
6874 suitable format.(1) The function returns `GSL_EFAILED' if there
6875 was a problem writing to the file.
6877 -- Function: int gsl_permutation_fscanf (FILE * STREAM,
6878 gsl_permutation * P)
6879 This function reads formatted data from the stream STREAM into the
6880 permutation P. The permutation P must be preallocated with the
6881 correct length since the function uses the size of P to determine
6882 how many numbers to read. The function returns `GSL_EFAILED' if
6883 there was a problem reading from the file.
6885 ---------- Footnotes ----------
6887 (1) In versions of the GNU C library prior to the ISO C99 standard,
6888 the type modifier `Z' was used instead.
6891 File: gsl-ref.info, Node: Permutations in cyclic form, Next: Permutation Examples, Prev: Reading and writing permutations, Up: Permutations
6893 9.8 Permutations in cyclic form
6894 ===============================
6896 A permutation can be represented in both "linear" and "cyclic"
6897 notations. The functions described in this section convert between the
6898 two forms. The linear notation is an index mapping, and has already
6899 been described above. The cyclic notation expresses a permutation as a
6900 series of circular rearrangements of groups of elements, or "cycles".
6902 For example, under the cycle (1 2 3), 1 is replaced by 2, 2 is
6903 replaced by 3 and 3 is replaced by 1 in a circular fashion. Cycles of
6904 different sets of elements can be combined independently, for example
6905 (1 2 3) (4 5) combines the cycle (1 2 3) with the cycle (4 5), which is
6906 an exchange of elements 4 and 5. A cycle of length one represents an
6907 element which is unchanged by the permutation and is referred to as a
6910 It can be shown that every permutation can be decomposed into
6911 combinations of cycles. The decomposition is not unique, but can always
6912 be rearranged into a standard "canonical form" by a reordering of
6913 elements. The library uses the canonical form defined in Knuth's `Art
6914 of Computer Programming' (Vol 1, 3rd Ed, 1997) Section 1.3.3, p.178.
6916 The procedure for obtaining the canonical form given by Knuth is,
6918 1. Write all singleton cycles explicitly
6920 2. Within each cycle, put the smallest number first
6922 3. Order the cycles in decreasing order of the first number in the
6925 For example, the linear representation (2 4 3 0 1) is represented as (1
6926 4) (0 2 3) in canonical form. The permutation corresponds to an
6927 exchange of elements 1 and 4, and rotation of elements 0, 2 and 3.
6929 The important property of the canonical form is that it can be
6930 reconstructed from the contents of each cycle without the brackets. In
6931 addition, by removing the brackets it can be considered as a linear
6932 representation of a different permutation. In the example given above
6933 the permutation (2 4 3 0 1) would become (1 4 0 2 3). This mapping has
6934 many applications in the theory of permutations.
6936 -- Function: int gsl_permutation_linear_to_canonical (gsl_permutation
6937 * Q, const gsl_permutation * P)
6938 This function computes the canonical form of the permutation P and
6939 stores it in the output argument Q.
6941 -- Function: int gsl_permutation_canonical_to_linear (gsl_permutation
6942 * P, const gsl_permutation * Q)
6943 This function converts a permutation Q in canonical form back into
6944 linear form storing it in the output argument P.
6946 -- Function: size_t gsl_permutation_inversions (const gsl_permutation
6948 This function counts the number of inversions in the permutation
6949 P. An inversion is any pair of elements that are not in order.
6950 For example, the permutation 2031 has three inversions,
6951 corresponding to the pairs (2,0) (2,1) and (3,1). The identity
6952 permutation has no inversions.
6954 -- Function: size_t gsl_permutation_linear_cycles (const
6955 gsl_permutation * P)
6956 This function counts the number of cycles in the permutation P,
6957 given in linear form.
6959 -- Function: size_t gsl_permutation_canonical_cycles (const
6960 gsl_permutation * Q)
6961 This function counts the number of cycles in the permutation Q,
6962 given in canonical form.
6965 File: gsl-ref.info, Node: Permutation Examples, Next: Permutation References and Further Reading, Prev: Permutations in cyclic form, Up: Permutations
6970 The example program below creates a random permutation (by shuffling the
6971 elements of the identity) and finds its inverse.
6974 #include <gsl/gsl_rng.h>
6975 #include <gsl/gsl_randist.h>
6976 #include <gsl/gsl_permutation.h>
6981 const size_t N = 10;
6982 const gsl_rng_type * T;
6985 gsl_permutation * p = gsl_permutation_alloc (N);
6986 gsl_permutation * q = gsl_permutation_alloc (N);
6988 gsl_rng_env_setup();
6989 T = gsl_rng_default;
6990 r = gsl_rng_alloc (T);
6992 printf ("initial permutation:");
6993 gsl_permutation_init (p);
6994 gsl_permutation_fprintf (stdout, p, " %u");
6997 printf (" random permutation:");
6998 gsl_ran_shuffle (r, p->data, N, sizeof(size_t));
6999 gsl_permutation_fprintf (stdout, p, " %u");
7002 printf ("inverse permutation:");
7003 gsl_permutation_inverse (q, p);
7004 gsl_permutation_fprintf (stdout, q, " %u");
7007 gsl_permutation_free (p);
7008 gsl_permutation_free (q);
7014 Here is the output from the program,
7017 initial permutation: 0 1 2 3 4 5 6 7 8 9
7018 random permutation: 1 3 5 2 7 6 0 4 9 8
7019 inverse permutation: 6 0 3 1 7 2 5 4 9 8
7021 The random permutation `p[i]' and its inverse `q[i]' are related
7022 through the identity `p[q[i]] = i', which can be verified from the
7025 The next example program steps forwards through all possible third
7026 order permutations, starting from the identity,
7029 #include <gsl/gsl_permutation.h>
7034 gsl_permutation * p = gsl_permutation_alloc (3);
7036 gsl_permutation_init (p);
7040 gsl_permutation_fprintf (stdout, p, " %u");
7043 while (gsl_permutation_next(p) == GSL_SUCCESS);
7045 gsl_permutation_free (p);
7050 Here is the output from the program,
7060 The permutations are generated in lexicographic order. To reverse the
7061 sequence, begin with the final permutation (which is the reverse of the
7062 identity) and replace `gsl_permutation_next' with
7063 `gsl_permutation_prev'.
7066 File: gsl-ref.info, Node: Permutation References and Further Reading, Prev: Permutation Examples, Up: Permutations
7068 9.10 References and Further Reading
7069 ===================================
7071 The subject of permutations is covered extensively in Knuth's `Sorting
7074 Donald E. Knuth, `The Art of Computer Programming: Sorting and
7075 Searching' (Vol 3, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850.
7077 For the definition of the "canonical form" see,
7079 Donald E. Knuth, `The Art of Computer Programming: Fundamental
7080 Algorithms' (Vol 1, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850.
7081 Section 1.3.3, `An Unusual Correspondence', p.178-179.
7084 File: gsl-ref.info, Node: Combinations, Next: Sorting, Prev: Permutations, Up: Top
7089 This chapter describes functions for creating and manipulating
7090 combinations. A combination c is represented by an array of k integers
7091 in the range 0 to n-1, where each value c_i occurs at most once. The
7092 combination c corresponds to indices of k elements chosen from an n
7093 element vector. Combinations are useful for iterating over all
7094 k-element subsets of a set.
7096 The functions described in this chapter are defined in the header
7097 file `gsl_combination.h'.
7101 * The Combination struct::
7102 * Combination allocation::
7103 * Accessing combination elements::
7104 * Combination properties::
7105 * Combination functions::
7106 * Reading and writing combinations::
7107 * Combination Examples::
7108 * Combination References and Further Reading::
7111 File: gsl-ref.info, Node: The Combination struct, Next: Combination allocation, Up: Combinations
7113 10.1 The Combination struct
7114 ===========================
7116 A combination is defined by a structure containing three components, the
7117 values of n and k, and a pointer to the combination array. The
7118 elements of the combination array are all of type `size_t', and are
7119 stored in increasing order. The `gsl_combination' structure looks like
7131 File: gsl-ref.info, Node: Combination allocation, Next: Accessing combination elements, Prev: The Combination struct, Up: Combinations
7133 10.2 Combination allocation
7134 ===========================
7136 -- Function: gsl_combination * gsl_combination_alloc (size_t N, size_t
7138 This function allocates memory for a new combination with
7139 parameters N, K. The combination is not initialized and its
7140 elements are undefined. Use the function `gsl_combination_calloc'
7141 if you want to create a combination which is initialized to the
7142 lexicographically first combination. A null pointer is returned if
7143 insufficient memory is available to create the combination.
7145 -- Function: gsl_combination * gsl_combination_calloc (size_t N,
7147 This function allocates memory for a new combination with
7148 parameters N, K and initializes it to the lexicographically first
7149 combination. A null pointer is returned if insufficient memory is
7150 available to create the combination.
7152 -- Function: void gsl_combination_init_first (gsl_combination * C)
7153 This function initializes the combination C to the
7154 lexicographically first combination, i.e. (0,1,2,...,k-1).
7156 -- Function: void gsl_combination_init_last (gsl_combination * C)
7157 This function initializes the combination C to the
7158 lexicographically last combination, i.e. (n-k,n-k+1,...,n-1).
7160 -- Function: void gsl_combination_free (gsl_combination * C)
7161 This function frees all the memory used by the combination C.
7163 -- Function: int gsl_combination_memcpy (gsl_combination * DEST, const
7164 gsl_combination * SRC)
7165 This function copies the elements of the combination SRC into the
7166 combination DEST. The two combinations must have the same size.
7169 File: gsl-ref.info, Node: Accessing combination elements, Next: Combination properties, Prev: Combination allocation, Up: Combinations
7171 10.3 Accessing combination elements
7172 ===================================
7174 The following function can be used to access the elements of a
7177 -- Function: size_t gsl_combination_get (const gsl_combination * C,
7179 This function returns the value of the I-th element of the
7180 combination C. If I lies outside the allowed range of 0 to K-1
7181 then the error handler is invoked and 0 is returned. An inline
7182 version of this function is used when `HAVE_INLINE' is defined.
7185 File: gsl-ref.info, Node: Combination properties, Next: Combination functions, Prev: Accessing combination elements, Up: Combinations
7187 10.4 Combination properties
7188 ===========================
7190 -- Function: size_t gsl_combination_n (const gsl_combination * C)
7191 This function returns the range (n) of the combination C.
7193 -- Function: size_t gsl_combination_k (const gsl_combination * C)
7194 This function returns the number of elements (k) in the
7197 -- Function: size_t * gsl_combination_data (const gsl_combination * C)
7198 This function returns a pointer to the array of elements in the
7201 -- Function: int gsl_combination_valid (gsl_combination * C)
7202 This function checks that the combination C is valid. The K
7203 elements should lie in the range 0 to N-1, with each value
7204 occurring once at most and in increasing order.
7207 File: gsl-ref.info, Node: Combination functions, Next: Reading and writing combinations, Prev: Combination properties, Up: Combinations
7209 10.5 Combination functions
7210 ==========================
7212 -- Function: int gsl_combination_next (gsl_combination * C)
7213 This function advances the combination C to the next combination
7214 in lexicographic order and returns `GSL_SUCCESS'. If no further
7215 combinations are available it returns `GSL_FAILURE' and leaves C
7216 unmodified. Starting with the first combination and repeatedly
7217 applying this function will iterate through all possible
7218 combinations of a given order.
7220 -- Function: int gsl_combination_prev (gsl_combination * C)
7221 This function steps backwards from the combination C to the
7222 previous combination in lexicographic order, returning
7223 `GSL_SUCCESS'. If no previous combination is available it returns
7224 `GSL_FAILURE' and leaves C unmodified.
7227 File: gsl-ref.info, Node: Reading and writing combinations, Next: Combination Examples, Prev: Combination functions, Up: Combinations
7229 10.6 Reading and writing combinations
7230 =====================================
7232 The library provides functions for reading and writing combinations to a
7233 file as binary data or formatted text.
7235 -- Function: int gsl_combination_fwrite (FILE * STREAM, const
7236 gsl_combination * C)
7237 This function writes the elements of the combination C to the
7238 stream STREAM in binary format. The function returns
7239 `GSL_EFAILED' if there was a problem writing to the file. Since
7240 the data is written in the native binary format it may not be
7241 portable between different architectures.
7243 -- Function: int gsl_combination_fread (FILE * STREAM, gsl_combination
7245 This function reads elements from the open stream STREAM into the
7246 combination C in binary format. The combination C must be
7247 preallocated with correct values of n and k since the function
7248 uses the size of C to determine how many bytes to read. The
7249 function returns `GSL_EFAILED' if there was a problem reading from
7250 the file. The data is assumed to have been written in the native
7251 binary format on the same architecture.
7253 -- Function: int gsl_combination_fprintf (FILE * STREAM, const
7254 gsl_combination * C, const char * FORMAT)
7255 This function writes the elements of the combination C
7256 line-by-line to the stream STREAM using the format specifier
7257 FORMAT, which should be suitable for a type of SIZE_T. In ISO C99
7258 the type modifier `z' represents `size_t', so `"%zu\n"' is a
7259 suitable format.(1) The function returns `GSL_EFAILED' if there
7260 was a problem writing to the file.
7262 -- Function: int gsl_combination_fscanf (FILE * STREAM,
7263 gsl_combination * C)
7264 This function reads formatted data from the stream STREAM into the
7265 combination C. The combination C must be preallocated with
7266 correct values of n and k since the function uses the size of C to
7267 determine how many numbers to read. The function returns
7268 `GSL_EFAILED' if there was a problem reading from the file.
7270 ---------- Footnotes ----------
7272 (1) In versions of the GNU C library prior to the ISO C99 standard,
7273 the type modifier `Z' was used instead.
7276 File: gsl-ref.info, Node: Combination Examples, Next: Combination References and Further Reading, Prev: Reading and writing combinations, Up: Combinations
7281 The example program below prints all subsets of the set {0,1,2,3}
7282 ordered by size. Subsets of the same size are ordered
7286 #include <gsl/gsl_combination.h>
7291 gsl_combination * c;
7294 printf ("All subsets of {0,1,2,3} by size:\n") ;
7295 for (i = 0; i <= 4; i++)
7297 c = gsl_combination_calloc (4, i);
7301 gsl_combination_fprintf (stdout, c, " %u");
7304 while (gsl_combination_next (c) == GSL_SUCCESS);
7305 gsl_combination_free (c);
7311 Here is the output from the program,
7314 All subsets of {0,1,2,3} by size:
7332 All 16 subsets are generated, and the subsets of each size are sorted
7336 File: gsl-ref.info, Node: Combination References and Further Reading, Prev: Combination Examples, Up: Combinations
7338 10.8 References and Further Reading
7339 ===================================
7341 Further information on combinations can be found in,
7343 Donald L. Kreher, Douglas R. Stinson, `Combinatorial Algorithms:
7344 Generation, Enumeration and Search', 1998, CRC Press LLC, ISBN
7349 File: gsl-ref.info, Node: Sorting, Next: BLAS Support, Prev: Combinations, Up: Top
7354 This chapter describes functions for sorting data, both directly and
7355 indirectly (using an index). All the functions use the "heapsort"
7356 algorithm. Heapsort is an O(N \log N) algorithm which operates
7357 in-place and does not require any additional storage. It also provides
7358 consistent performance, the running time for its worst-case (ordered
7359 data) being not significantly longer than the average and best cases.
7360 Note that the heapsort algorithm does not preserve the relative ordering
7361 of equal elements--it is an "unstable" sort. However the resulting
7362 order of equal elements will be consistent across different platforms
7363 when using these functions.
7369 * Selecting the k smallest or largest elements::
7370 * Computing the rank::
7371 * Sorting Examples::
7372 * Sorting References and Further Reading::
7375 File: gsl-ref.info, Node: Sorting objects, Next: Sorting vectors, Up: Sorting
7377 11.1 Sorting objects
7378 ====================
7380 The following function provides a simple alternative to the standard
7381 library function `qsort'. It is intended for systems lacking `qsort',
7382 not as a replacement for it. The function `qsort' should be used
7383 whenever possible, as it will be faster and can provide stable ordering
7384 of equal elements. Documentation for `qsort' is available in the `GNU
7385 C Library Reference Manual'.
7387 The functions described in this section are defined in the header
7388 file `gsl_heapsort.h'.
7390 -- Function: void gsl_heapsort (void * ARRAY, size_t COUNT, size_t
7391 SIZE, gsl_comparison_fn_t COMPARE)
7392 This function sorts the COUNT elements of the array ARRAY, each of
7393 size SIZE, into ascending order using the comparison function
7394 COMPARE. The type of the comparison function is defined by,
7396 int (*gsl_comparison_fn_t) (const void * a,
7399 A comparison function should return a negative integer if the first
7400 argument is less than the second argument, `0' if the two arguments
7401 are equal and a positive integer if the first argument is greater
7402 than the second argument.
7404 For example, the following function can be used to sort doubles
7405 into ascending numerical order.
7408 compare_doubles (const double * a,
7419 The appropriate function call to perform the sort is,
7421 gsl_heapsort (array, count, sizeof(double),
7424 Note that unlike `qsort' the heapsort algorithm cannot be made into
7425 a stable sort by pointer arithmetic. The trick of comparing
7426 pointers for equal elements in the comparison function does not
7427 work for the heapsort algorithm. The heapsort algorithm performs
7428 an internal rearrangement of the data which destroys its initial
7431 -- Function: int gsl_heapsort_index (size_t * P, const void * ARRAY,
7432 size_t COUNT, size_t SIZE, gsl_comparison_fn_t COMPARE)
7433 This function indirectly sorts the COUNT elements of the array
7434 ARRAY, each of size SIZE, into ascending order using the
7435 comparison function COMPARE. The resulting permutation is stored
7436 in P, an array of length N. The elements of P give the index of
7437 the array element which would have been stored in that position if
7438 the array had been sorted in place. The first element of P gives
7439 the index of the least element in ARRAY, and the last element of P
7440 gives the index of the greatest element in ARRAY. The array
7441 itself is not changed.