diff options
Diffstat (limited to 'gsl-1.9/doc/gsl-ref.info-2')
| -rw-r--r-- | gsl-1.9/doc/gsl-ref.info-2 | 2862 |
1 files changed, 1435 insertions, 1427 deletions
diff --git a/gsl-1.9/doc/gsl-ref.info-2 b/gsl-1.9/doc/gsl-ref.info-2 index 3fd15e9..f27f665 100644 --- a/gsl-1.9/doc/gsl-ref.info-2 +++ b/gsl-1.9/doc/gsl-ref.info-2 @@ -1,12 +1,7 @@ -This is gsl-ref.info, produced by makeinfo version 4.8 from +This is gsl-ref.info, produced by makeinfo version 5.1 from gsl-ref.texi. -INFO-DIR-SECTION Scientific software -START-INFO-DIR-ENTRY -* gsl-ref: (gsl-ref). GNU Scientific Library - Reference -END-INFO-DIR-ENTRY - - Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, +Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 The GSL Team. Permission is granted to copy, distribute and/or modify this document @@ -20,6 +15,10 @@ License". (a) The Back-Cover Text is: "You have freedom to copy and modify this GNU Manual, like GNU software." +INFO-DIR-SECTION Scientific software +START-INFO-DIR-ENTRY +* gsl-ref: (gsl-ref). GNU Scientific Library - Reference +END-INFO-DIR-ENTRY File: gsl-ref.info, Node: Selecting the k smallest or largest elements, Next: Computing the rank, Prev: Sorting vectors, Up: Sorting @@ -27,47 +26,47 @@ File: gsl-ref.info, Node: Selecting the k smallest or largest elements, Next: 11.3 Selecting the k smallest or largest elements ================================================= -The functions described in this section select the k smallest or -largest elements of a data set of size N. The routines use an O(kN) -direct insertion algorithm which is suited to subsets that are small -compared with the total size of the dataset. For example, the routines -are useful for selecting the 10 largest values from one million data -points, but not for selecting the largest 100,000 values. If the -subset is a significant part of the total dataset it may be faster to -sort all the elements of the dataset directly with an O(N \log N) -algorithm and obtain the smallest or largest values that way. +The functions described in this section select the k smallest or largest +elements of a data set of size N. The routines use an O(kN) direct +insertion algorithm which is suited to subsets that are small compared +with the total size of the dataset. For example, the routines are +useful for selecting the 10 largest values from one million data points, +but not for selecting the largest 100,000 values. If the subset is a +significant part of the total dataset it may be faster to sort all the +elements of the dataset directly with an O(N \log N) algorithm and +obtain the smallest or largest values that way. -- Function: int gsl_sort_smallest (double * DEST, size_t K, const double * SRC, size_t STRIDE, size_t N) This function copies the K smallest elements of the array SRC, of size N and stride STRIDE, in ascending numerical order into the - array DEST. The size K of the subset must be less than or equal - to N. The data SRC is not modified by this operation. + array DEST. The size K of the subset must be less than or equal to + N. The data SRC is not modified by this operation. -- Function: int gsl_sort_largest (double * DEST, size_t K, const double * SRC, size_t STRIDE, size_t N) This function copies the K largest elements of the array SRC, of size N and stride STRIDE, in descending numerical order into the - array DEST. K must be less than or equal to N. The data SRC is not - modified by this operation. + array DEST. K must be less than or equal to N. The data SRC is + not modified by this operation. -- Function: int gsl_sort_vector_smallest (double * DEST, size_t K, const gsl_vector * V) -- Function: int gsl_sort_vector_largest (double * DEST, size_t K, const gsl_vector * V) These functions copy the K smallest or largest elements of the - vector V into the array DEST. K must be less than or equal to the + vector V into the array DEST. K must be less than or equal to the length of the vector V. - The following functions find the indices of the k smallest or -largest elements of a dataset, + The following functions find the indices of the k smallest or largest +elements of a dataset, -- Function: int gsl_sort_smallest_index (size_t * P, size_t K, const double * SRC, size_t STRIDE, size_t N) This function stores the indices of the K smallest elements of the array SRC, of size N and stride STRIDE, in the array P. The indices are chosen so that the corresponding data is in ascending - numerical order. K must be less than or equal to N. The data SRC + numerical order. K must be less than or equal to N. The data SRC is not modified by this operation. -- Function: int gsl_sort_largest_index (size_t * P, size_t K, const @@ -75,7 +74,7 @@ largest elements of a dataset, This function stores the indices of the K largest elements of the array SRC, of size N and stride STRIDE, in the array P. The indices are chosen so that the corresponding data is in descending - numerical order. K must be less than or equal to N. The data SRC + numerical order. K must be less than or equal to N. The data SRC is not modified by this operation. -- Function: int gsl_sort_vector_smallest_index (size_t * P, size_t K, @@ -83,7 +82,7 @@ largest elements of a dataset, -- Function: int gsl_sort_vector_largest_index (size_t * P, size_t K, const gsl_vector * V) These functions store the indices of the K smallest or largest - elements of the vector V in the array P. K must be less than or + elements of the vector V in the array P. K must be less than or equal to the length of the vector V. @@ -103,7 +102,7 @@ following algorithm, } This can be computed directly from the function -`gsl_permutation_inverse(rank,p)'. +'gsl_permutation_inverse(rank,p)'. The following function will print the rank of each element of the vector V, @@ -147,7 +146,7 @@ elements of the vector V in ascending order, printf ("order = %d, value = %g\n", i, vpi); } -The next example uses the function `gsl_sort_smallest' to select the 5 +The next example uses the function 'gsl_sort_smallest' to select the 5 smallest numbers from 100000 uniform random variates stored in an array, #include <gsl/gsl_rng.h> @@ -204,15 +203,15 @@ File: gsl-ref.info, Node: Sorting References and Further Reading, Prev: Sortin 11.6 References and Further Reading =================================== -The subject of sorting is covered extensively in Knuth's `Sorting and +The subject of sorting is covered extensively in Knuth's 'Sorting and Searching', - Donald E. Knuth, `The Art of Computer Programming: Sorting and + Donald E. Knuth, 'The Art of Computer Programming: Sorting and Searching' (Vol 3, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896850. The Heapsort algorithm is described in the following book, - Robert Sedgewick, `Algorithms in C', Addison-Wesley, ISBN + Robert Sedgewick, 'Algorithms in C', Addison-Wesley, ISBN 0201514257. @@ -230,34 +229,31 @@ the C-language BLAS standard, referred to here as "CBLAS", and a higher-level interface for operations on GSL vectors and matrices. Users who are interested in simple operations on GSL vector and matrix objects should use the high-level layer, which is declared in the file -`gsl_blas.h'. This should satisfy the needs of most users. Note that +'gsl_blas.h'. This should satisfy the needs of most users. Note that GSL matrices are implemented using dense-storage so the interface only includes the corresponding dense-storage BLAS functions. The full BLAS functionality for band-format and packed-format matrices is available through the low-level CBLAS interface. - The interface for the `gsl_cblas' layer is specified in the file -`gsl_cblas.h'. This interface corresponds to the BLAS Technical -Forum's draft standard for the C interface to legacy BLAS -implementations. Users who have access to other conforming CBLAS -implementations can use these in place of the version provided by the -library. Note that users who have only a Fortran BLAS library can use -a CBLAS conformant wrapper to convert it into a CBLAS library. A -reference CBLAS wrapper for legacy Fortran implementations exists as -part of the draft CBLAS standard and can be obtained from Netlib. The -complete set of CBLAS functions is listed in an appendix (*note GSL -CBLAS Library::). + The interface for the 'gsl_cblas' layer is specified in the file +'gsl_cblas.h'. This interface corresponds to the BLAS Technical Forum's +draft standard for the C interface to legacy BLAS implementations. +Users who have access to other conforming CBLAS implementations can use +these in place of the version provided by the library. Note that users +who have only a Fortran BLAS library can use a CBLAS conformant wrapper +to convert it into a CBLAS library. A reference CBLAS wrapper for +legacy Fortran implementations exists as part of the draft CBLAS +standard and can be obtained from Netlib. The complete set of CBLAS +functions is listed in an appendix (*note GSL CBLAS Library::). There are three levels of BLAS operations, Level 1 - Vector operations, e.g. y = \alpha x + y - + Vector operations, e.g. y = \alpha x + y Level 2 - Matrix-vector operations, e.g. y = \alpha A x + \beta y - + Matrix-vector operations, e.g. y = \alpha A x + \beta y Level 3 - Matrix-matrix operations, e.g. C = \alpha A B + C + Matrix-matrix operations, e.g. C = \alpha A B + C Each routine has a name which specifies the operation, the type of matrices involved and their precisions. Some of the most common @@ -265,19 +261,14 @@ operations and their names are given below, DOT scalar product, x^T y - AXPY vector sum, \alpha x + y - MV matrix-vector product, A x - SV matrix-vector solve, inv(A) x - MM matrix-matrix product, A B - SM matrix-matrix solve, inv(A) B @@ -285,34 +276,24 @@ The types of matrices are, GE general - GB general band - SY symmetric - SB symmetric band - SP symmetric packed - HE hermitian - HB hermitian band - HP hermitian packed - TR triangular - TB triangular band - TP triangular packed @@ -320,13 +301,10 @@ Each operation is defined for four precisions, S single real - D double real - C single complex - Z double complex @@ -336,9 +314,9 @@ matrix-matrix multiply". * Menu: -* GSL BLAS Interface:: -* BLAS Examples:: -* BLAS References and Further Reading:: +* GSL BLAS Interface:: +* BLAS Examples:: +* BLAS References and Further Reading:: File: gsl-ref.info, Node: GSL BLAS Interface, Next: BLAS Examples, Up: BLAS Support @@ -349,13 +327,13 @@ File: gsl-ref.info, Node: GSL BLAS Interface, Next: BLAS Examples, Up: BLAS S GSL provides dense vector and matrix objects, based on the relevant built-in types. The library provides an interface to the BLAS operations which apply to these objects. The interface to this -functionality is given in the file `gsl_blas.h'. +functionality is given in the file 'gsl_blas.h'. * Menu: -* Level 1 GSL BLAS Interface:: -* Level 2 GSL BLAS Interface:: -* Level 3 GSL BLAS Interface:: +* Level 1 GSL BLAS Interface:: +* Level 2 GSL BLAS Interface:: +* Level 3 GSL BLAS Interface:: File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS Interface, Up: GSL BLAS Interface @@ -363,10 +341,10 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I 12.1.1 Level 1 -------------- - -- Function: int gsl_blas_sdsdot (float ALPHA, const gsl_vector_float - * X, const gsl_vector_float * Y, float * RESULT) - This function computes the sum \alpha + x^T y for the vectors X - and Y, returning the result in RESULT. + -- Function: int gsl_blas_sdsdot (float ALPHA, const gsl_vector_float * + X, const gsl_vector_float * Y, float * RESULT) + This function computes the sum \alpha + x^T y for the vectors X and + Y, returning the result in RESULT. -- Function: int gsl_blas_sdot (const gsl_vector_float * X, const gsl_vector_float * Y, float * RESULT) @@ -404,8 +382,8 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I -- Function: float gsl_blas_sasum (const gsl_vector_float * X) -- Function: double gsl_blas_dasum (const gsl_vector * X) - These functions compute the absolute sum \sum |x_i| of the - elements of the vector X. + These functions compute the absolute sum \sum |x_i| of the elements + of the vector X. -- Function: float gsl_blas_scasum (const gsl_vector_complex_float * X) -- Function: double gsl_blas_dzasum (const gsl_vector_complex * X) @@ -420,14 +398,14 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I -- Function: CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * X) These functions return the index of the largest element of the - vector X. The largest element is determined by its absolute + vector X. The largest element is determined by its absolute magnitude for real vectors and by the sum of the magnitudes of the real and imaginary parts |\Re(x_i)| + |\Im(x_i)| for complex vectors. If the largest value occurs several times then the index of the first occurrence is returned. - -- Function: int gsl_blas_sswap (gsl_vector_float * X, - gsl_vector_float * Y) + -- Function: int gsl_blas_sswap (gsl_vector_float * X, gsl_vector_float + * Y) -- Function: int gsl_blas_dswap (gsl_vector * X, gsl_vector * Y) -- Function: int gsl_blas_cswap (gsl_vector_complex_float * X, gsl_vector_complex_float * Y) @@ -469,8 +447,8 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I These functions rescale the vector X by the multiplicative factor ALPHA. - -- Function: int gsl_blas_srotg (float A[], float B[], float C[], - float S[]) + -- Function: int gsl_blas_srotg (float A[], float B[], float C[], float + S[]) -- Function: int gsl_blas_drotg (double A[], double B[], double C[], double S[]) These functions compute a Givens rotation (c,s) which zeroes the @@ -485,8 +463,8 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I * Y, float C, float S) -- Function: int gsl_blas_drot (gsl_vector * X, gsl_vector * Y, const double C, const double S) - These functions apply a Givens rotation (x', y') = (c x + s y, -s - x + c y) to the vectors X, Y. + These functions apply a Givens rotation (x', y') = (c x + s y, -s x + + c y) to the vectors X, Y. -- Function: int gsl_blas_srotmg (float D1[], float D2[], float B1[], float B2, float P[]) @@ -496,8 +474,8 @@ File: gsl-ref.info, Node: Level 1 GSL BLAS Interface, Next: Level 2 GSL BLAS I modified Givens transformation is defined in the original Level-1 BLAS specification, given in the references. - -- Function: int gsl_blas_srotm (gsl_vector_float * X, - gsl_vector_float * Y, const float P[]) + -- Function: int gsl_blas_srotm (gsl_vector_float * X, gsl_vector_float + * Y, const float P[]) -- Function: int gsl_blas_drotm (gsl_vector * X, gsl_vector * Y, const double P[]) These functions apply a modified Givens transformation. @@ -508,9 +486,9 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I 12.1.2 Level 2 -------------- - -- Function: int gsl_blas_sgemv (CBLAS_TRANSPOSE_t TRANSA, float - ALPHA, const gsl_matrix_float * A, const gsl_vector_float * - X, float BETA, gsl_vector_float * Y) + -- Function: int gsl_blas_sgemv (CBLAS_TRANSPOSE_t TRANSA, float ALPHA, + const gsl_matrix_float * A, const gsl_vector_float * X, float + BETA, gsl_vector_float * Y) -- Function: int gsl_blas_dgemv (CBLAS_TRANSPOSE_t TRANSA, double ALPHA, const gsl_matrix * A, const gsl_vector * X, double BETA, gsl_vector * Y) @@ -524,7 +502,7 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I gsl_vector_complex * Y) These functions compute the matrix-vector product and sum y = \alpha op(A) x + \beta y, where op(A) = A, A^T, A^H for TRANSA = - `CblasNoTrans', `CblasTrans', `CblasConjTrans'. + 'CblasNoTrans', 'CblasTrans', 'CblasConjTrans'. -- Function: int gsl_blas_strmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_float * A, @@ -533,19 +511,19 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix * A, gsl_vector * X) -- Function: int gsl_blas_ctrmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t - TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float * - A, gsl_vector_complex_float * X) + TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float * A, + gsl_vector_complex_float * X) -- Function: int gsl_blas_ztrmv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex * A, gsl_vector_complex * X) These functions compute the matrix-vector product x = op(A) x for the triangular matrix A, where op(A) = A, A^T, A^H for TRANSA = - `CblasNoTrans', `CblasTrans', `CblasConjTrans'. When UPLO is - `CblasUpper' then the upper triangle of A is used, and when UPLO - is `CblasLower' then the lower triangle of A is used. If DIAG is - `CblasNonUnit' then the diagonal of the matrix is used, but if - DIAG is `CblasUnit' then the diagonal elements of the matrix A are - taken as unity and are not referenced. + 'CblasNoTrans', 'CblasTrans', 'CblasConjTrans'. When UPLO is + 'CblasUpper' then the upper triangle of A is used, and when UPLO is + 'CblasLower' then the lower triangle of A is used. If DIAG is + 'CblasNonUnit' then the diagonal of the matrix is used, but if DIAG + is 'CblasUnit' then the diagonal elements of the matrix A are taken + as unity and are not referenced. -- Function: int gsl_blas_strsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_float * A, @@ -554,31 +532,31 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix * A, gsl_vector * X) -- Function: int gsl_blas_ctrsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t - TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float * - A, gsl_vector_complex_float * X) + TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex_float * A, + gsl_vector_complex_float * X) -- Function: int gsl_blas_ztrsv (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_matrix_complex * A, gsl_vector_complex * X) These functions compute inv(op(A)) x for X, where op(A) = A, A^T, - A^H for TRANSA = `CblasNoTrans', `CblasTrans', `CblasConjTrans'. - When UPLO is `CblasUpper' then the upper triangle of A is used, - and when UPLO is `CblasLower' then the lower triangle of A is - used. If DIAG is `CblasNonUnit' then the diagonal of the matrix - is used, but if DIAG is `CblasUnit' then the diagonal elements of - the matrix A are taken as unity and are not referenced. + A^H for TRANSA = 'CblasNoTrans', 'CblasTrans', 'CblasConjTrans'. + When UPLO is 'CblasUpper' then the upper triangle of A is used, and + when UPLO is 'CblasLower' then the lower triangle of A is used. If + DIAG is 'CblasNonUnit' then the diagonal of the matrix is used, but + if DIAG is 'CblasUnit' then the diagonal elements of the matrix A + are taken as unity and are not referenced. -- Function: int gsl_blas_ssymv (CBLAS_UPLO_t UPLO, float ALPHA, const gsl_matrix_float * A, const gsl_vector_float * X, float BETA, gsl_vector_float * Y) - -- Function: int gsl_blas_dsymv (CBLAS_UPLO_t UPLO, double ALPHA, - const gsl_matrix * A, const gsl_vector * X, double BETA, - gsl_vector * Y) + -- Function: int gsl_blas_dsymv (CBLAS_UPLO_t UPLO, double ALPHA, const + gsl_matrix * A, const gsl_vector * X, double BETA, gsl_vector + * Y) These functions compute the matrix-vector product and sum y = \alpha A x + \beta y for the symmetric matrix A. Since the matrix - A is symmetric only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of A are used, and when UPLO is `CblasLower' then the - lower triangle and diagonal of A are used. + A is symmetric only its upper half or lower half need to be stored. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. -- Function: int gsl_blas_chemv (CBLAS_UPLO_t UPLO, const gsl_complex_float ALPHA, const gsl_matrix_complex_float * A, @@ -589,20 +567,19 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I * X, const gsl_complex BETA, gsl_vector_complex * Y) These functions compute the matrix-vector product and sum y = \alpha A x + \beta y for the hermitian matrix A. Since the matrix - A is hermitian only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of A are used, and when UPLO is `CblasLower' then the - lower triangle and diagonal of A are used. The imaginary elements - of the diagonal are automatically assumed to be zero and are not - referenced. + A is hermitian only its upper half or lower half need to be stored. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. The imaginary elements of the diagonal are + automatically assumed to be zero and are not referenced. -- Function: int gsl_blas_sger (float ALPHA, const gsl_vector_float * X, const gsl_vector_float * Y, gsl_matrix_float * A) -- Function: int gsl_blas_dger (double ALPHA, const gsl_vector * X, const gsl_vector * Y, gsl_matrix * A) -- Function: int gsl_blas_cgeru (const gsl_complex_float ALPHA, const - gsl_vector_complex_float * X, const gsl_vector_complex_float - * Y, gsl_matrix_complex_float * A) + gsl_vector_complex_float * X, const gsl_vector_complex_float * + Y, gsl_matrix_complex_float * A) -- Function: int gsl_blas_zgeru (const gsl_complex ALPHA, const gsl_vector_complex * X, const gsl_vector_complex * Y, gsl_matrix_complex * A) @@ -610,8 +587,8 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I the matrix A. -- Function: int gsl_blas_cgerc (const gsl_complex_float ALPHA, const - gsl_vector_complex_float * X, const gsl_vector_complex_float - * Y, gsl_matrix_complex_float * A) + gsl_vector_complex_float * X, const gsl_vector_complex_float * + Y, gsl_matrix_complex_float * A) -- Function: int gsl_blas_zgerc (const gsl_complex ALPHA, const gsl_vector_complex * X, const gsl_vector_complex * Y, gsl_matrix_complex * A) @@ -623,50 +600,50 @@ File: gsl-ref.info, Node: Level 2 GSL BLAS Interface, Next: Level 3 GSL BLAS I -- Function: int gsl_blas_dsyr (CBLAS_UPLO_t UPLO, double ALPHA, const gsl_vector * X, gsl_matrix * A) These functions compute the symmetric rank-1 update A = \alpha x - x^T + A of the symmetric matrix A. Since the matrix A is - symmetric only its upper half or lower half need to be stored. - When UPLO is `CblasUpper' then the upper triangle and diagonal of - A are used, and when UPLO is `CblasLower' then the lower triangle - and diagonal of A are used. + x^T + A of the symmetric matrix A. Since the matrix A is symmetric + only its upper half or lower half need to be stored. When UPLO is + 'CblasUpper' then the upper triangle and diagonal of A are used, + and when UPLO is 'CblasLower' then the lower triangle and diagonal + of A are used. -- Function: int gsl_blas_cher (CBLAS_UPLO_t UPLO, float ALPHA, const gsl_vector_complex_float * X, gsl_matrix_complex_float * A) -- Function: int gsl_blas_zher (CBLAS_UPLO_t UPLO, double ALPHA, const gsl_vector_complex * X, gsl_matrix_complex * A) These functions compute the hermitian rank-1 update A = \alpha x - x^H + A of the hermitian matrix A. Since the matrix A is - hermitian only its upper half or lower half need to be stored. - When UPLO is `CblasUpper' then the upper triangle and diagonal of - A are used, and when UPLO is `CblasLower' then the lower triangle - and diagonal of A are used. The imaginary elements of the - diagonal are automatically set to zero. + x^H + A of the hermitian matrix A. Since the matrix A is hermitian + only its upper half or lower half need to be stored. When UPLO is + 'CblasUpper' then the upper triangle and diagonal of A are used, + and when UPLO is 'CblasLower' then the lower triangle and diagonal + of A are used. The imaginary elements of the diagonal are + automatically set to zero. -- Function: int gsl_blas_ssyr2 (CBLAS_UPLO_t UPLO, float ALPHA, const gsl_vector_float * X, const gsl_vector_float * Y, gsl_matrix_float * A) - -- Function: int gsl_blas_dsyr2 (CBLAS_UPLO_t UPLO, double ALPHA, - const gsl_vector * X, const gsl_vector * Y, gsl_matrix * A) + -- Function: int gsl_blas_dsyr2 (CBLAS_UPLO_t UPLO, double ALPHA, const + gsl_vector * X, const gsl_vector * Y, gsl_matrix * A) These functions compute the symmetric rank-2 update A = \alpha x - y^T + \alpha y x^T + A of the symmetric matrix A. Since the - matrix A is symmetric only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of A are used, and when UPLO is `CblasLower' then the - lower triangle and diagonal of A are used. + y^T + \alpha y x^T + A of the symmetric matrix A. Since the matrix + A is symmetric only its upper half or lower half need to be stored. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. -- Function: int gsl_blas_cher2 (CBLAS_UPLO_t UPLO, const gsl_complex_float ALPHA, const gsl_vector_complex_float * X, - const gsl_vector_complex_float * Y, gsl_matrix_complex_float - * A) + const gsl_vector_complex_float * Y, gsl_matrix_complex_float * + A) -- Function: int gsl_blas_zher2 (CBLAS_UPLO_t UPLO, const gsl_complex ALPHA, const gsl_vector_complex * X, const gsl_vector_complex * Y, gsl_matrix_complex * A) These functions compute the hermitian rank-2 update A = \alpha x - y^H + \alpha^* y x^H A of the hermitian matrix A. Since the - matrix A is hermitian only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of A are used, and when UPLO is `CblasLower' then the - lower triangle and diagonal of A are used. The imaginary elements - of the diagonal are automatically set to zero. + y^H + \alpha^* y x^H A of the hermitian matrix A. Since the matrix + A is hermitian only its upper half or lower half need to be stored. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. The imaginary elements of the diagonal are + automatically set to zero. File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS Interface, Up: GSL BLAS Interface @@ -682,17 +659,16 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I CBLAS_TRANSPOSE_t TRANSB, double ALPHA, const gsl_matrix * A, const gsl_matrix * B, double BETA, gsl_matrix * C) -- Function: int gsl_blas_cgemm (CBLAS_TRANSPOSE_t TRANSA, - CBLAS_TRANSPOSE_t TRANSB, const gsl_complex_float ALPHA, - const gsl_matrix_complex_float * A, const - gsl_matrix_complex_float * B, const gsl_complex_float BETA, - gsl_matrix_complex_float * C) + CBLAS_TRANSPOSE_t TRANSB, const gsl_complex_float ALPHA, const + gsl_matrix_complex_float * A, const gsl_matrix_complex_float * + B, const gsl_complex_float BETA, gsl_matrix_complex_float * C) -- Function: int gsl_blas_zgemm (CBLAS_TRANSPOSE_t TRANSA, CBLAS_TRANSPOSE_t TRANSB, const gsl_complex ALPHA, const gsl_matrix_complex * A, const gsl_matrix_complex * B, const gsl_complex BETA, gsl_matrix_complex * C) These functions compute the matrix-matrix product and sum C = - \alpha op(A) op(B) + \beta C where op(A) = A, A^T, A^H for TRANSA - = `CblasNoTrans', `CblasTrans', `CblasConjTrans' and similarly for + \alpha op(A) op(B) + \beta C where op(A) = A, A^T, A^H for TRANSA = + 'CblasNoTrans', 'CblasTrans', 'CblasConjTrans' and similarly for the parameter TRANSB. -- Function: int gsl_blas_ssymm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, @@ -710,11 +686,11 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I gsl_matrix_complex * B, const gsl_complex BETA, gsl_matrix_complex * C) These functions compute the matrix-matrix product and sum C = - \alpha A B + \beta C for SIDE is `CblasLeft' and C = \alpha B A + - \beta C for SIDE is `CblasRight', where the matrix A is symmetric. - When UPLO is `CblasUpper' then the upper triangle and diagonal of - A are used, and when UPLO is `CblasLower' then the lower triangle - and diagonal of A are used. + \alpha A B + \beta C for SIDE is 'CblasLeft' and C = \alpha B A + + \beta C for SIDE is 'CblasRight', where the matrix A is symmetric. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. -- Function: int gsl_blas_chemm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, const gsl_complex_float ALPHA, const gsl_matrix_complex_float @@ -725,12 +701,12 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I gsl_matrix_complex * B, const gsl_complex BETA, gsl_matrix_complex * C) These functions compute the matrix-matrix product and sum C = - \alpha A B + \beta C for SIDE is `CblasLeft' and C = \alpha B A + - \beta C for SIDE is `CblasRight', where the matrix A is hermitian. - When UPLO is `CblasUpper' then the upper triangle and diagonal of - A are used, and when UPLO is `CblasLower' then the lower triangle - and diagonal of A are used. The imaginary elements of the - diagonal are automatically set to zero. + \alpha A B + \beta C for SIDE is 'CblasLeft' and C = \alpha B A + + \beta C for SIDE is 'CblasRight', where the matrix A is hermitian. + When UPLO is 'CblasUpper' then the upper triangle and diagonal of A + are used, and when UPLO is 'CblasLower' then the lower triangle and + diagonal of A are used. The imaginary elements of the diagonal are + automatically set to zero. -- Function: int gsl_blas_strmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, float ALPHA, @@ -743,18 +719,17 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I gsl_complex_float ALPHA, const gsl_matrix_complex_float * A, gsl_matrix_complex_float * B) -- Function: int gsl_blas_ztrmm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, - CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const - gsl_complex ALPHA, const gsl_matrix_complex * A, - gsl_matrix_complex * B) + CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_complex + ALPHA, const gsl_matrix_complex * A, gsl_matrix_complex * B) These functions compute the matrix-matrix product B = \alpha op(A) - B for SIDE is `CblasLeft' and B = \alpha B op(A) for SIDE is - `CblasRight'. The matrix A is triangular and op(A) = A, A^T, A^H - for TRANSA = `CblasNoTrans', `CblasTrans', `CblasConjTrans'. When - UPLO is `CblasUpper' then the upper triangle of A is used, and - when UPLO is `CblasLower' then the lower triangle of A is used. - If DIAG is `CblasNonUnit' then the diagonal of A is used, but if - DIAG is `CblasUnit' then the diagonal elements of the matrix A are - taken as unity and are not referenced. + B for SIDE is 'CblasLeft' and B = \alpha B op(A) for SIDE is + 'CblasRight'. The matrix A is triangular and op(A) = A, A^T, A^H + for TRANSA = 'CblasNoTrans', 'CblasTrans', 'CblasConjTrans'. When + UPLO is 'CblasUpper' then the upper triangle of A is used, and when + UPLO is 'CblasLower' then the lower triangle of A is used. If DIAG + is 'CblasNonUnit' then the diagonal of A is used, but if DIAG is + 'CblasUnit' then the diagonal elements of the matrix A are taken as + unity and are not referenced. -- Function: int gsl_blas_strsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, float ALPHA, @@ -767,19 +742,17 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I gsl_complex_float ALPHA, const gsl_matrix_complex_float * A, gsl_matrix_complex_float * B) -- Function: int gsl_blas_ztrsm (CBLAS_SIDE_t SIDE, CBLAS_UPLO_t UPLO, - CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const - gsl_complex ALPHA, const gsl_matrix_complex * A, - gsl_matrix_complex * B) + CBLAS_TRANSPOSE_t TRANSA, CBLAS_DIAG_t DIAG, const gsl_complex + ALPHA, const gsl_matrix_complex * A, gsl_matrix_complex * B) These functions compute the inverse-matrix matrix product B = - \alpha op(inv(A))B for SIDE is `CblasLeft' and B = \alpha B - op(inv(A)) for SIDE is `CblasRight'. The matrix A is triangular - and op(A) = A, A^T, A^H for TRANSA = `CblasNoTrans', `CblasTrans', - `CblasConjTrans'. When UPLO is `CblasUpper' then the upper - triangle of A is used, and when UPLO is `CblasLower' then the - lower triangle of A is used. If DIAG is `CblasNonUnit' then the - diagonal of A is used, but if DIAG is `CblasUnit' then the - diagonal elements of the matrix A are taken as unity and are not - referenced. + \alpha op(inv(A))B for SIDE is 'CblasLeft' and B = \alpha B + op(inv(A)) for SIDE is 'CblasRight'. The matrix A is triangular + and op(A) = A, A^T, A^H for TRANSA = 'CblasNoTrans', 'CblasTrans', + 'CblasConjTrans'. When UPLO is 'CblasUpper' then the upper + triangle of A is used, and when UPLO is 'CblasLower' then the lower + triangle of A is used. If DIAG is 'CblasNonUnit' then the diagonal + of A is used, but if DIAG is 'CblasUnit' then the diagonal elements + of the matrix A are taken as unity and are not referenced. -- Function: int gsl_blas_ssyrk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANS, float ALPHA, const gsl_matrix_float * A, float BETA, @@ -795,11 +768,11 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A, const gsl_complex BETA, gsl_matrix_complex * C) These functions compute a rank-k update of the symmetric matrix C, - C = \alpha A A^T + \beta C when TRANS is `CblasNoTrans' and C = - \alpha A^T A + \beta C when TRANS is `CblasTrans'. Since the + C = \alpha A A^T + \beta C when TRANS is 'CblasNoTrans' and C = + \alpha A^T A + \beta C when TRANS is 'CblasTrans'. Since the matrix C is symmetric only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of C are used, and when UPLO is `CblasLower' then the + stored. When UPLO is 'CblasUpper' then the upper triangle and + diagonal of C are used, and when UPLO is 'CblasLower' then the lower triangle and diagonal of C are used. -- Function: int gsl_blas_cherk (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t @@ -809,11 +782,11 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I TRANS, double ALPHA, const gsl_matrix_complex * A, double BETA, gsl_matrix_complex * C) These functions compute a rank-k update of the hermitian matrix C, - C = \alpha A A^H + \beta C when TRANS is `CblasNoTrans' and C = - \alpha A^H A + \beta C when TRANS is `CblasTrans'. Since the + C = \alpha A A^H + \beta C when TRANS is 'CblasNoTrans' and C = + \alpha A^H A + \beta C when TRANS is 'CblasTrans'. Since the matrix C is hermitian only its upper half or lower half need to be - stored. When UPLO is `CblasUpper' then the upper triangle and - diagonal of C are used, and when UPLO is `CblasLower' then the + stored. When UPLO is 'CblasUpper' then the upper triangle and + diagonal of C are used, and when UPLO is 'CblasLower' then the lower triangle and diagonal of C are used. The imaginary elements of the diagonal are automatically set to zero. @@ -825,37 +798,36 @@ File: gsl-ref.info, Node: Level 3 GSL BLAS Interface, Prev: Level 2 GSL BLAS I B, double BETA, gsl_matrix * C) -- Function: int gsl_blas_csyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANS, const gsl_complex_float ALPHA, const - gsl_matrix_complex_float * A, const gsl_matrix_complex_float - * B, const gsl_complex_float BETA, gsl_matrix_complex_float * - C) + gsl_matrix_complex_float * A, const gsl_matrix_complex_float * + B, const gsl_complex_float BETA, gsl_matrix_complex_float * C) -- Function: int gsl_blas_zsyr2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A, const gsl_matrix_complex * B, const gsl_complex BETA, gsl_matrix_complex * C) These functions compute a rank-2k update of the symmetric matrix C, C = \alpha A B^T + \alpha B A^T + \beta C when TRANS is - `CblasNoTrans' and C = \alpha A^T B + \alpha B^T A + \beta C when - TRANS is `CblasTrans'. Since the matrix C is symmetric only its + 'CblasNoTrans' and C = \alpha A^T B + \alpha B^T A + \beta C when + TRANS is 'CblasTrans'. Since the matrix C is symmetric only its upper half or lower half need to be stored. When UPLO is - `CblasUpper' then the upper triangle and diagonal of C are used, - and when UPLO is `CblasLower' then the lower triangle and diagonal + 'CblasUpper' then the upper triangle and diagonal of C are used, + and when UPLO is 'CblasLower' then the lower triangle and diagonal of C are used. -- Function: int gsl_blas_cher2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANS, const gsl_complex_float ALPHA, const - gsl_matrix_complex_float * A, const gsl_matrix_complex_float - * B, float BETA, gsl_matrix_complex_float * C) + gsl_matrix_complex_float * A, const gsl_matrix_complex_float * + B, float BETA, gsl_matrix_complex_float * C) -- Function: int gsl_blas_zher2k (CBLAS_UPLO_t UPLO, CBLAS_TRANSPOSE_t TRANS, const gsl_complex ALPHA, const gsl_matrix_complex * A, const gsl_matrix_complex * B, double BETA, gsl_matrix_complex * C) These functions compute a rank-2k update of the hermitian matrix C, C = \alpha A B^H + \alpha^* B A^H + \beta C when TRANS is - `CblasNoTrans' and C = \alpha A^H B + \alpha^* B^H A + \beta C when - TRANS is `CblasConjTrans'. Since the matrix C is hermitian only + 'CblasNoTrans' and C = \alpha A^H B + \alpha^* B^H A + \beta C when + TRANS is 'CblasConjTrans'. Since the matrix C is hermitian only its upper half or lower half need to be stored. When UPLO is - `CblasUpper' then the upper triangle and diagonal of C are used, - and when UPLO is `CblasLower' then the lower triangle and diagonal + 'CblasUpper' then the upper triangle and diagonal of C are used, + and when UPLO is 'CblasLower' then the lower triangle and diagonal of C are used. The imaginary elements of the diagonal are automatically set to zero. @@ -920,33 +892,31 @@ File: gsl-ref.info, Node: BLAS References and Further Reading, Prev: BLAS Exam =================================== Information on the BLAS standards, including both the legacy and draft -interface standards, is available online from the BLAS Homepage and -BLAS Technical Forum web-site. +interface standards, is available online from the BLAS Homepage and BLAS +Technical Forum web-site. - `BLAS Homepage' - `http://www.netlib.org/blas/' - - `BLAS Technical Forum' - `http://www.netlib.org/cgi-bin/checkout/blast/blast.pl' + 'BLAS Homepage' + <http://www.netlib.org/blas/> + 'BLAS Technical Forum' + <http://www.netlib.org/cgi-bin/checkout/blast/blast.pl> The following papers contain the specifications for Level 1, Level 2 and Level 3 BLAS. C. Lawson, R. Hanson, D. Kincaid, F. Krogh, "Basic Linear Algebra - Subprograms for Fortran Usage", `ACM Transactions on Mathematical + Subprograms for Fortran Usage", 'ACM Transactions on Mathematical Software', Vol. 5 (1979), Pages 308-325. J.J. Dongarra, J. DuCroz, S. Hammarling, R. Hanson, "An Extended - Set of Fortran Basic Linear Algebra Subprograms", `ACM - Transactions on Mathematical Software', Vol. 14, No. 1 (1988), - Pages 1-32. + Set of Fortran Basic Linear Algebra Subprograms", 'ACM Transactions + on Mathematical Software', Vol. 14, No. 1 (1988), Pages 1-32. - J.J. Dongarra, I. Duff, J. DuCroz, S. Hammarling, "A Set of Level - 3 Basic Linear Algebra Subprograms", `ACM Transactions on + J.J. Dongarra, I. Duff, J. DuCroz, S. Hammarling, "A Set of Level 3 + Basic Linear Algebra Subprograms", 'ACM Transactions on Mathematical Software', Vol. 16 (1990), Pages 1-28. Postscript versions of the latter two papers are available from -`http://www.netlib.org/blas/'. A CBLAS wrapper for Fortran BLAS +<http://www.netlib.org/blas/>. A CBLAS wrapper for Fortran BLAS libraries is available from the same location. @@ -956,34 +926,34 @@ File: gsl-ref.info, Node: Linear Algebra, Next: Eigensystems, Prev: BLAS Supp ***************** This chapter describes functions for solving linear systems. The -library provides linear algebra operations which operate directly on -the `gsl_vector' and `gsl_matrix' objects. These routines use the -standard algorithms from Golub & Van Loan's `Matrix Computations'. +library provides linear algebra operations which operate directly on the +'gsl_vector' and 'gsl_matrix' objects. These routines use the standard +algorithms from Golub & Van Loan's 'Matrix Computations'. When dealing with very large systems the routines found in LAPACK should be considered. These support specialized data representations and other optimizations. The functions described in this chapter are declared in the header -file `gsl_linalg.h'. +file 'gsl_linalg.h'. * Menu: -* LU Decomposition:: -* QR Decomposition:: -* QR Decomposition with Column Pivoting:: -* Singular Value Decomposition:: -* Cholesky Decomposition:: -* Tridiagonal Decomposition of Real Symmetric Matrices:: -* Tridiagonal Decomposition of Hermitian Matrices:: +* LU Decomposition:: +* QR Decomposition:: +* QR Decomposition with Column Pivoting:: +* Singular Value Decomposition:: +* Cholesky Decomposition:: +* Tridiagonal Decomposition of Real Symmetric Matrices:: +* Tridiagonal Decomposition of Hermitian Matrices:: * Hessenberg Decomposition of Real Matrices:: -* Bidiagonalization:: -* Householder Transformations:: -* Householder solver for linear systems:: -* Tridiagonal Systems:: +* Bidiagonalization:: +* Householder Transformations:: +* Householder solver for linear systems:: +* Tridiagonal Systems:: * Balancing:: -* Linear Algebra Examples:: -* Linear Algebra References and Further Reading:: +* Linear Algebra Examples:: +* Linear Algebra References and Further Reading:: File: gsl-ref.info, Node: LU Decomposition, Next: QR Decomposition, Up: Linear Algebra @@ -996,68 +966,67 @@ triangular matrices, P A = L U -where P is a permutation matrix, L is unit lower triangular matrix and -U is upper triangular matrix. For square matrices this decomposition -can be used to convert the linear system A x = b into a pair of -triangular systems (L y = P b, U x = y), which can be solved by forward -and back-substitution. Note that the LU decomposition is valid for -singular matrices. +where P is a permutation matrix, L is unit lower triangular matrix and U +is upper triangular matrix. For square matrices this decomposition can +be used to convert the linear system A x = b into a pair of triangular +systems (L y = P b, U x = y), which can be solved by forward and +back-substitution. Note that the LU decomposition is valid for singular +matrices. -- Function: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation * P, int * SIGNUM) -- Function: int gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * P, int * SIGNUM) These functions factorize the square matrix A into the LU - decomposition PA = LU. On output the diagonal and upper - triangular part of the input matrix A contain the matrix U. The - lower triangular part of the input matrix (excluding the diagonal) + decomposition PA = LU. On output the diagonal and upper triangular + part of the input matrix A contain the matrix U. The lower + triangular part of the input matrix (excluding the diagonal) contains L. The diagonal elements of L are unity, and are not stored. - The permutation matrix P is encoded in the permutation P. The j-th + The permutation matrix P is encoded in the permutation P. The j-th column of the matrix P is given by the k-th column of the identity matrix, where k = p_j the j-th element of the permutation vector. - The sign of the permutation is given by SIGNUM. It has the value + The sign of the permutation is given by SIGNUM. It has the value (-1)^n, where n is the number of interchanges in the permutation. The algorithm used in the decomposition is Gaussian Elimination - with partial pivoting (Golub & Van Loan, `Matrix Computations', + with partial pivoting (Golub & Van Loan, 'Matrix Computations', Algorithm 3.4.1). -- Function: int gsl_linalg_LU_solve (const gsl_matrix * LU, const gsl_permutation * P, const gsl_vector * B, gsl_vector * X) -- Function: int gsl_linalg_complex_LU_solve (const gsl_matrix_complex - * LU, const gsl_permutation * P, const gsl_vector_complex * - B, gsl_vector_complex * X) + * LU, const gsl_permutation * P, const gsl_vector_complex * B, + gsl_vector_complex * X) These functions solve the square system A x = b using the LU - decomposition of A into (LU, P) given by `gsl_linalg_LU_decomp' or - `gsl_linalg_complex_LU_decomp'. + decomposition of A into (LU, P) given by 'gsl_linalg_LU_decomp' or + 'gsl_linalg_complex_LU_decomp'. -- Function: int gsl_linalg_LU_svx (const gsl_matrix * LU, const gsl_permutation * P, gsl_vector * X) -- Function: int gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * P, gsl_vector_complex * X) These functions solve the square system A x = b in-place using the - LU decomposition of A into (LU,P). On input X should contain the + LU decomposition of A into (LU,P). On input X should contain the right-hand side b, which is replaced by the solution on output. -- Function: int gsl_linalg_LU_refine (const gsl_matrix * A, const - gsl_matrix * LU, const gsl_permutation * P, const gsl_vector - * B, gsl_vector * X, gsl_vector * RESIDUAL) - -- Function: int gsl_linalg_complex_LU_refine (const - gsl_matrix_complex * A, const gsl_matrix_complex * LU, const - gsl_permutation * P, const gsl_vector_complex * B, - gsl_vector_complex * X, gsl_vector_complex * RESIDUAL) + gsl_matrix * LU, const gsl_permutation * P, const gsl_vector * + B, gsl_vector * X, gsl_vector * RESIDUAL) + -- Function: int gsl_linalg_complex_LU_refine (const gsl_matrix_complex + * A, const gsl_matrix_complex * LU, const gsl_permutation * P, + const gsl_vector_complex * B, gsl_vector_complex * X, + gsl_vector_complex * RESIDUAL) These functions apply an iterative improvement to X, the solution - of A x = b, using the LU decomposition of A into (LU,P). The + of A x = b, using the LU decomposition of A into (LU,P). The initial residual r = A x - b is also computed and stored in RESIDUAL. -- Function: int gsl_linalg_LU_invert (const gsl_matrix * LU, const gsl_permutation * P, gsl_matrix * INVERSE) - -- Function: int gsl_linalg_complex_LU_invert (const - gsl_matrix_complex * LU, const gsl_permutation * P, - gsl_matrix_complex * INVERSE) + -- Function: int gsl_linalg_complex_LU_invert (const gsl_matrix_complex + * LU, const gsl_permutation * P, gsl_matrix_complex * INVERSE) These functions compute the inverse of a matrix A from its LU decomposition (LU,P), storing the result in the matrix INVERSE. The inverse is computed by solving the system A x = b for each @@ -1071,7 +1040,7 @@ singular matrices. -- Function: gsl_complex gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int SIGNUM) These functions compute the determinant of a matrix A from its LU - decomposition, LU. The determinant is computed as the product of + decomposition, LU. The determinant is computed as the product of the diagonal elements of U and the sign of the row permutation SIGNUM. @@ -1079,16 +1048,15 @@ singular matrices. -- Function: double gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU) These functions compute the logarithm of the absolute value of the - determinant of a matrix A, \ln|\det(A)|, from its LU - decomposition, LU. This function may be useful if the direct - computation of the determinant would overflow or underflow. + determinant of a matrix A, \ln|\det(A)|, from its LU decomposition, + LU. This function may be useful if the direct computation of the + determinant would overflow or underflow. -- Function: int gsl_linalg_LU_sgndet (gsl_matrix * LU, int SIGNUM) -- Function: gsl_complex gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int SIGNUM) - These functions compute the sign or phase factor of the - determinant of a matrix A, \det(A)/|\det(A)|, from its LU - decomposition, LU. + These functions compute the sign or phase factor of the determinant + of a matrix A, \det(A)/|\det(A)|, from its LU decomposition, LU. File: gsl-ref.info, Node: QR Decomposition, Next: QR Decomposition with Column Pivoting, Prev: LU Decomposition, Up: Linear Algebra @@ -1097,47 +1065,46 @@ File: gsl-ref.info, Node: QR Decomposition, Next: QR Decomposition with Column ===================== A general rectangular M-by-N matrix A has a QR decomposition into the -product of an orthogonal M-by-M square matrix Q (where Q^T Q = I) and -an M-by-N right-triangular matrix R, +product of an orthogonal M-by-M square matrix Q (where Q^T Q = I) and an +M-by-N right-triangular matrix R, A = Q R -This decomposition can be used to convert the linear system A x = b -into the triangular system R x = Q^T b, which can be solved by -back-substitution. Another use of the QR decomposition is to compute an -orthonormal basis for a set of vectors. The first N columns of Q form +This decomposition can be used to convert the linear system A x = b into +the triangular system R x = Q^T b, which can be solved by +back-substitution. Another use of the QR decomposition is to compute an +orthonormal basis for a set of vectors. The first N columns of Q form an orthonormal basis for the range of A, ran(A), when A has full column rank. -- Function: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * TAU) This function factorizes the M-by-N matrix A into the QR - decomposition A = Q R. On output the diagonal and upper - triangular part of the input matrix contain the matrix R. The - vector TAU and the columns of the lower triangular part of the - matrix A contain the Householder coefficients and Householder - vectors which encode the orthogonal matrix Q. The vector TAU must - be of length k=\min(M,N). The matrix Q is related to these - components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i - v_i^T and v_i is the Householder vector v_i = - (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage - scheme as used by LAPACK. + decomposition A = Q R. On output the diagonal and upper triangular + part of the input matrix contain the matrix R. The vector TAU and + the columns of the lower triangular part of the matrix A contain + the Householder coefficients and Householder vectors which encode + the orthogonal matrix Q. The vector TAU must be of length + k=\min(M,N). The matrix Q is related to these components by, Q = + Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the + Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). + This is the same storage scheme as used by LAPACK. The algorithm used to perform the decomposition is Householder QR - (Golub & Van Loan, `Matrix Computations', Algorithm 5.2.1). + (Golub & Van Loan, 'Matrix Computations', Algorithm 5.2.1). -- Function: int gsl_linalg_QR_solve (const gsl_matrix * QR, const gsl_vector * TAU, const gsl_vector * B, gsl_vector * X) This function solves the square system A x = b using the QR - decomposition of A into (QR, TAU) given by `gsl_linalg_QR_decomp'. + decomposition of A into (QR, TAU) given by 'gsl_linalg_QR_decomp'. The least-squares solution for rectangular systems can be found - using `gsl_linalg_QR_lssolve'. + using 'gsl_linalg_QR_lssolve'. -- Function: int gsl_linalg_QR_svx (const gsl_matrix * QR, const gsl_vector * TAU, gsl_vector * X) This function solves the square system A x = b in-place using the QR decomposition of A into (QR,TAU) given by - `gsl_linalg_QR_decomp'. On input X should contain the right-hand + 'gsl_linalg_QR_decomp'. On input X should contain the right-hand side b, which is replaced by the solution on output. -- Function: int gsl_linalg_QR_lssolve (const gsl_matrix * QR, const @@ -1147,7 +1114,7 @@ rank. overdetermined system A x = b where the matrix A has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual, ||Ax - b||.The routine uses the QR decomposition - of A into (QR, TAU) given by `gsl_linalg_QR_decomp'. The solution + of A into (QR, TAU) given by 'gsl_linalg_QR_decomp'. The solution is returned in X. The residual is computed as a by-product and stored in RESIDUAL. @@ -1155,8 +1122,8 @@ rank. gsl_vector * TAU, gsl_vector * V) This function applies the matrix Q^T encoded in the decomposition (QR,TAU) to the vector V, storing the result Q^T v in V. The - matrix multiplication is carried out directly using the encoding - of the Householder vectors without needing to form the full matrix + matrix multiplication is carried out directly using the encoding of + the Householder vectors without needing to form the full matrix Q^T. -- Function: int gsl_linalg_QR_Qvec (const gsl_matrix * QR, const @@ -1168,34 +1135,34 @@ rank. -- Function: int gsl_linalg_QR_Rsolve (const gsl_matrix * QR, const gsl_vector * B, gsl_vector * X) - This function solves the triangular system R x = b for X. It may - be useful if the product b' = Q^T b has already been computed - using `gsl_linalg_QR_QTvec'. + This function solves the triangular system R x = b for X. It may + be useful if the product b' = Q^T b has already been computed using + 'gsl_linalg_QR_QTvec'. -- Function: int gsl_linalg_QR_Rsvx (const gsl_matrix * QR, gsl_vector * X) This function solves the triangular system R x = b for X in-place. On input X should contain the right-hand side b and is replaced by - the solution on output. This function may be useful if the product - b' = Q^T b has already been computed using `gsl_linalg_QR_QTvec'. + the solution on output. This function may be useful if the product + b' = Q^T b has already been computed using 'gsl_linalg_QR_QTvec'. -- Function: int gsl_linalg_QR_unpack (const gsl_matrix * QR, const gsl_vector * TAU, gsl_matrix * Q, gsl_matrix * R) This function unpacks the encoded QR decomposition (QR,TAU) into the matrices Q and R, where Q is M-by-M and R is M-by-N. - -- Function: int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * - R, const gsl_vector * B, gsl_vector * X) - This function solves the system R x = Q^T b for X. It can be used - when the QR decomposition of a matrix is available in unpacked - form as (Q, R). + -- Function: int gsl_linalg_QR_QRsolve (gsl_matrix * Q, gsl_matrix * R, + const gsl_vector * B, gsl_vector * X) + This function solves the system R x = Q^T b for X. It can be used + when the QR decomposition of a matrix is available in unpacked form + as (Q, R). -- Function: int gsl_linalg_QR_update (gsl_matrix * Q, gsl_matrix * R, gsl_vector * W, const gsl_vector * V) This function performs a rank-1 update w v^T of the QR - decomposition (Q, R). The update is given by Q'R' = Q R + w v^T + decomposition (Q, R). The update is given by Q'R' = Q R + w v^T where the output matrices Q' and R' are also orthogonal and right - triangular. Note that W is destroyed by the update. + triangular. Note that W is destroyed by the update. -- Function: int gsl_linalg_R_solve (const gsl_matrix * R, const gsl_vector * B, gsl_vector * X) @@ -1204,7 +1171,7 @@ rank. -- Function: int gsl_linalg_R_svx (const gsl_matrix * R, gsl_vector * X) - This function solves the triangular system R x = b in-place. On + This function solves the triangular system R x = b in-place. On input X should contain the right-hand side b, which is replaced by the solution on output. @@ -1220,31 +1187,31 @@ introducing a column permutation P, A P = Q R The first r columns of Q form an orthonormal basis for the range of A -for a matrix with column rank r. This decomposition can also be used -to convert the linear system A x = b into the triangular system R y = -Q^T b, x = P y, which can be solved by back-substitution and -permutation. We denote the QR decomposition with column pivoting by -QRP^T since A = Q R P^T. +for a matrix with column rank r. This decomposition can also be used to +convert the linear system A x = b into the triangular system R y = Q^T +b, x = P y, which can be solved by back-substitution and permutation. +We denote the QR decomposition with column pivoting by QRP^T since A = Q +R P^T. -- Function: int gsl_linalg_QRPT_decomp (gsl_matrix * A, gsl_vector * TAU, gsl_permutation * P, int * SIGNUM, gsl_vector * NORM) This function factorizes the M-by-N matrix A into the QRP^T decomposition A = Q R P^T. On output the diagonal and upper - triangular part of the input matrix contain the matrix R. The + triangular part of the input matrix contain the matrix R. The permutation matrix P is stored in the permutation P. The sign of - the permutation is given by SIGNUM. It has the value (-1)^n, where - n is the number of interchanges in the permutation. The vector TAU + the permutation is given by SIGNUM. It has the value (-1)^n, where + n is the number of interchanges in the permutation. The vector TAU and the columns of the lower triangular part of the matrix A contain the Householder coefficients and vectors which encode the - orthogonal matrix Q. The vector TAU must be of length - k=\min(M,N). The matrix Q is related to these components by, Q = - Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the - Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). - This is the same storage scheme as used by LAPACK. The vector - NORM is a workspace of length N used for column pivoting. + orthogonal matrix Q. The vector TAU must be of length k=\min(M,N). + The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1 + where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder vector + v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same + storage scheme as used by LAPACK. The vector NORM is a workspace + of length N used for column pivoting. The algorithm used to perform the decomposition is Householder QR - with column pivoting (Golub & Van Loan, `Matrix Computations', + with column pivoting (Golub & Van Loan, 'Matrix Computations', Algorithm 5.4.1). -- Function: int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A, @@ -1259,29 +1226,29 @@ QRP^T since A = Q R P^T. * B, gsl_vector * X) This function solves the square system A x = b using the QRP^T decomposition of A into (QR, TAU, P) given by - `gsl_linalg_QRPT_decomp'. + 'gsl_linalg_QRPT_decomp'. -- Function: int gsl_linalg_QRPT_svx (const gsl_matrix * QR, const gsl_vector * TAU, const gsl_permutation * P, gsl_vector * X) This function solves the square system A x = b in-place using the - QRP^T decomposition of A into (QR,TAU,P). On input X should - contain the right-hand side b, which is replaced by the solution - on output. + QRP^T decomposition of A into (QR,TAU,P). On input X should + contain the right-hand side b, which is replaced by the solution on + output. -- Function: int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, const gsl_matrix * R, const gsl_permutation * P, const gsl_vector * B, gsl_vector * X) - This function solves the square system R P^T x = Q^T b for X. It + This function solves the square system R P^T x = Q^T b for X. It can be used when the QR decomposition of a matrix is available in unpacked form as (Q, R). -- Function: int gsl_linalg_QRPT_update (gsl_matrix * Q, gsl_matrix * - R, const gsl_permutation * P, gsl_vector * U, const - gsl_vector * V) + R, const gsl_permutation * P, gsl_vector * U, const gsl_vector + * V) This function performs a rank-1 update w v^T of the QRP^T - decomposition (Q, R, P). The update is given by Q'R' = Q R + w v^T + decomposition (Q, R, P). The update is given by Q'R' = Q R + w v^T where the output matrices Q' and R' are also orthogonal and right - triangular. Note that W is destroyed by the update. The + triangular. Note that W is destroyed by the update. The permutation P is not changed. -- Function: int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR, const @@ -1291,9 +1258,9 @@ QRP^T since A = Q R P^T. -- Function: int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR, const gsl_permutation * P, gsl_vector * X) - This function solves the triangular system R P^T x = b in-place - for the N-by-N matrix R contained in QR. On input X should contain - the right-hand side b, which is replaced by the solution on output. + This function solves the triangular system R P^T x = b in-place for + the N-by-N matrix R contained in QR. On input X should contain the + right-hand side b, which is replaced by the solution on output. File: gsl-ref.info, Node: Singular Value Decomposition, Next: Cholesky Decomposition, Prev: QR Decomposition with Column Pivoting, Up: Linear Algebra @@ -1301,10 +1268,10 @@ File: gsl-ref.info, Node: Singular Value Decomposition, Next: Cholesky Decompo 13.4 Singular Value Decomposition ================================= -A general rectangular M-by-N matrix A has a singular value -decomposition (SVD) into the product of an M-by-N orthogonal matrix U, -an N-by-N diagonal matrix of singular values S and the transpose of an -N-by-N orthogonal square matrix V, +A general rectangular M-by-N matrix A has a singular value decomposition +(SVD) into the product of an M-by-N orthogonal matrix U, an N-by-N +diagonal matrix of singular values S and the transpose of an N-by-N +orthogonal square matrix V, A = U S V^T @@ -1314,8 +1281,8 @@ generally chosen to form a non-increasing sequence \sigma_1 >= \sigma_2 The singular value decomposition of a matrix has many practical uses. The condition number of the matrix is given by the ratio of the largest -singular value to the smallest singular value. The presence of a zero -singular value indicates that the matrix is singular. The number of +singular value to the smallest singular value. The presence of a zero +singular value indicates that the matrix is singular. The number of non-zero singular values indicates the rank of the matrix. In practice singular value decomposition of a rank-deficient matrix will not produce exact zeroes for singular values, due to finite numerical precision. @@ -1329,18 +1296,18 @@ singular values. -- Function: int gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * WORK) This function factorizes the M-by-N matrix A into the singular - value decomposition A = U S V^T for M >= N. On output the matrix - A is replaced by U. The diagonal elements of the singular value - matrix S are stored in the vector S. The singular values are + value decomposition A = U S V^T for M >= N. On output the matrix A + is replaced by U. The diagonal elements of the singular value + matrix S are stored in the vector S. The singular values are non-negative and form a non-increasing sequence from S_1 to S_N. - The matrix V contains the elements of V in untransposed form. To + The matrix V contains the elements of V in untransposed form. To form the product U S V^T it is necessary to take the transpose of V. A workspace of length N is required in WORK. This routine uses the Golub-Reinsch SVD algorithm. - -- Function: int gsl_linalg_SV_decomp_mod (gsl_matrix * A, gsl_matrix - * X, gsl_matrix * V, gsl_vector * S, gsl_vector * WORK) + -- Function: int gsl_linalg_SV_decomp_mod (gsl_matrix * A, gsl_matrix * + X, gsl_matrix * V, gsl_vector * S, gsl_vector * WORK) This function computes the SVD using the modified Golub-Reinsch algorithm, which is faster for M>>N. It requires the vector WORK of length N and the N-by-N matrix X as additional working space. @@ -1355,7 +1322,7 @@ singular values. -- Function: int gsl_linalg_SV_solve (gsl_matrix * U, gsl_matrix * V, gsl_vector * S, const gsl_vector * B, gsl_vector * X) This function solves the system A x = b using the singular value - decomposition (U, S, V) of A given by `gsl_linalg_SV_decomp'. + decomposition (U, S, V) of A given by 'gsl_linalg_SV_decomp'. Only non-zero singular values are used in computing the solution. The parts of the solution corresponding to singular values of zero @@ -1378,33 +1345,34 @@ transpose L^T, A = L L^T -This is sometimes referred to as taking the square-root of a matrix. The -Cholesky decomposition can only be carried out when all the eigenvalues -of the matrix are positive. This decomposition can be used to convert -the linear system A x = b into a pair of triangular systems (L y = b, -L^T x = y), which can be solved by forward and back-substitution. +This is sometimes referred to as taking the square-root of a matrix. +The Cholesky decomposition can only be carried out when all the +eigenvalues of the matrix are positive. This decomposition can be used +to convert the linear system A x = b into a pair of triangular systems +(L y = b, L^T x = y), which can be solved by forward and +back-substitution. -- Function: int gsl_linalg_cholesky_decomp (gsl_matrix * A) This function factorizes the positive-definite symmetric square - matrix A into the Cholesky decomposition A = L L^T. On input only + matrix A into the Cholesky decomposition A = L L^T. On input only the diagonal and lower-triangular part of the matrix A are needed. On output the diagonal and lower triangular part of the input - matrix A contain the matrix L. The upper triangular part of the + matrix A contain the matrix L. The upper triangular part of the input matrix contains L^T, the diagonal terms being identical for both L and L^T. If the matrix is not positive-definite then the - decomposition will fail, returning the error code `GSL_EDOM'. + decomposition will fail, returning the error code 'GSL_EDOM'. -- Function: int gsl_linalg_cholesky_solve (const gsl_matrix * CHOLESKY, const gsl_vector * B, gsl_vector * X) This function solves the system A x = b using the Cholesky decomposition of A into the matrix CHOLESKY given by - `gsl_linalg_cholesky_decomp'. + 'gsl_linalg_cholesky_decomp'. -- Function: int gsl_linalg_cholesky_svx (const gsl_matrix * CHOLESKY, gsl_vector * X) - This function solves the system A x = b in-place using the - Cholesky decomposition of A into the matrix CHOLESKY given by - `gsl_linalg_cholesky_decomp'. On input X should contain the + This function solves the system A x = b in-place using the Cholesky + decomposition of A into the matrix CHOLESKY given by + 'gsl_linalg_cholesky_decomp'. On input X should contain the right-hand side b, which is replaced by the solution on output. @@ -1420,8 +1388,8 @@ into the form, where Q is an orthogonal matrix and T is a symmetric tridiagonal matrix. - -- Function: int gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector - * TAU) + -- Function: int gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * + TAU) This function factorizes the symmetric square matrix A into the symmetric tridiagonal decomposition Q T Q^T. On output the diagonal and subdiagonal part of the input matrix A contain the @@ -1435,7 +1403,7 @@ where Q is an orthogonal matrix and T is a symmetric tridiagonal matrix. gsl_vector * TAU, gsl_matrix * Q, gsl_vector * DIAG, gsl_vector * SUBDIAG) This function unpacks the encoded symmetric tridiagonal - decomposition (A, TAU) obtained from `gsl_linalg_symmtd_decomp' + decomposition (A, TAU) obtained from 'gsl_linalg_symmtd_decomp' into the orthogonal matrix Q, the vector of diagonal elements DIAG and the vector of subdiagonal elements SUBDIAG. @@ -1443,7 +1411,7 @@ where Q is an orthogonal matrix and T is a symmetric tridiagonal matrix. gsl_vector * DIAG, gsl_vector * SUBDIAG) This function unpacks the diagonal and subdiagonal of the encoded symmetric tridiagonal decomposition (A, TAU) obtained from - `gsl_linalg_symmtd_decomp' into the vectors DIAG and SUBDIAG. + 'gsl_linalg_symmtd_decomp' into the vectors DIAG and SUBDIAG. File: gsl-ref.info, Node: Tridiagonal Decomposition of Hermitian Matrices, Next: Hessenberg Decomposition of Real Matrices, Prev: Tridiagonal Decomposition of Real Symmetric Matrices, Up: Linear Algebra @@ -1462,12 +1430,12 @@ matrix. -- Function: int gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * TAU) This function factorizes the hermitian matrix A into the symmetric - tridiagonal decomposition U T U^T. On output the real parts of - the diagonal and subdiagonal part of the input matrix A contain - the tridiagonal matrix T. The remaining lower triangular part of - the input matrix contains the Householder vectors which, together - with the Householder coefficients TAU, encode the orthogonal matrix - Q. This storage scheme is the same as used by LAPACK. The upper + tridiagonal decomposition U T U^T. On output the real parts of the + diagonal and subdiagonal part of the input matrix A contain the + tridiagonal matrix T. The remaining lower triangular part of the + input matrix contains the Householder vectors which, together with + the Householder coefficients TAU, encode the orthogonal matrix Q. + This storage scheme is the same as used by LAPACK. The upper triangular part of A and imaginary parts of the diagonal are not referenced. @@ -1475,15 +1443,15 @@ matrix. A, const gsl_vector_complex * TAU, gsl_matrix_complex * Q, gsl_vector * DIAG, gsl_vector * SUBDIAG) This function unpacks the encoded tridiagonal decomposition (A, - TAU) obtained from `gsl_linalg_hermtd_decomp' into the unitary + TAU) obtained from 'gsl_linalg_hermtd_decomp' into the unitary matrix U, the real vector of diagonal elements DIAG and the real vector of subdiagonal elements SUBDIAG. - -- Function: int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex - * A, gsl_vector * DIAG, gsl_vector * SUBDIAG) + -- Function: int gsl_linalg_hermtd_unpack_T (const gsl_matrix_complex * + A, gsl_vector * DIAG, gsl_vector * SUBDIAG) This function unpacks the diagonal and subdiagonal of the encoded tridiagonal decomposition (A, TAU) obtained from the - `gsl_linalg_hermtd_decomp' into the real vectors DIAG and SUBDIAG. + 'gsl_linalg_hermtd_decomp' into the real vectors DIAG and SUBDIAG. File: gsl-ref.info, Node: Hessenberg Decomposition of Real Matrices, Next: Bidiagonalization, Prev: Tridiagonal Decomposition of Hermitian Matrices, Up: Linear Algebra @@ -1497,40 +1465,40 @@ transformations into the form A = U H U^T where U is orthogonal and H is an upper Hessenberg matrix, meaning -that it has zeros below the first subdiagonal. The Hessenberg reduction +that it has zeros below the first subdiagonal. The Hessenberg reduction is the first step in the Schur decomposition for the nonsymmetric eigenvalue problem, but has applications in other areas as well. -- Function: int gsl_linalg_hessenberg (gsl_matrix * A, gsl_vector * TAU) - This function computes the Hessenberg decomposition of the matrix - A by applying the similarity transformation H = U^T A U. On - output, H is stored in the upper portion of A. The information - required to construct the matrix U is stored in the lower - triangular portion of A. U is a product of N - 2 Householder - matrices. The Householder vectors are stored in the lower portion - of A (below the subdiagonal) and the Householder coefficients are - stored in the vector TAU. TAU must be of length N. + This function computes the Hessenberg decomposition of the matrix A + by applying the similarity transformation H = U^T A U. On output, + H is stored in the upper portion of A. The information required to + construct the matrix U is stored in the lower triangular portion of + A. U is a product of N - 2 Householder matrices. The Householder + vectors are stored in the lower portion of A (below the + subdiagonal) and the Householder coefficients are stored in the + vector TAU. TAU must be of length N. -- Function: int gsl_linalg_hessenberg_unpack (gsl_matrix * H, gsl_vector * TAU, gsl_matrix * U) This function constructs the orthogonal matrix U from the - information stored in the Hessenberg matrix H along with the - vector TAU. H and TAU are outputs from `gsl_linalg_hessenberg'. + information stored in the Hessenberg matrix H along with the vector + TAU. H and TAU are outputs from 'gsl_linalg_hessenberg'. -- Function: int gsl_linalg_hessenberg_unpack_accum (gsl_matrix * H, gsl_vector * TAU, gsl_matrix * V) - This function is similar to `gsl_linalg_hessenberg_unpack', except + This function is similar to 'gsl_linalg_hessenberg_unpack', except it accumulates the matrix U into V, so that V' = VU. The matrix V must be initialized prior to calling this function. Setting V to the identity matrix provides the same result as - `gsl_linalg_hessenberg_unpack'. If H is order N, then V must have + 'gsl_linalg_hessenberg_unpack'. If H is order N, then V must have N columns but may have any number of rows. -- Function: void gsl_linalg_hessenberg_set_zero (gsl_matrix * H) This function sets the lower triangular portion of H, below the - subdiagonal, to zero. It is useful for clearing out the - Householder vectors after calling `gsl_linalg_hessenberg'. + subdiagonal, to zero. It is useful for clearing out the + Householder vectors after calling 'gsl_linalg_hessenberg'. File: gsl-ref.info, Node: Bidiagonalization, Next: Householder Transformations, Prev: Hessenberg Decomposition of Real Matrices, Up: Linear Algebra @@ -1547,37 +1515,37 @@ where U and V are orthogonal matrices and B is a N-by-N bidiagonal matrix with non-zero entries only on the diagonal and superdiagonal. The size of U is M-by-N and the size of V is N-by-N. - -- Function: int gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector - * TAU_U, gsl_vector * TAU_V) - This function factorizes the M-by-N matrix A into bidiagonal form - U B V^T. The diagonal and superdiagonal of the matrix B are - stored in the diagonal and superdiagonal of A. The orthogonal - matrices U and V are stored as compressed Householder vectors in - the remaining elements of A. The Householder coefficients are - stored in the vectors TAU_U and TAU_V. The length of TAU_U must - equal the number of elements in the diagonal of A and the length - of TAU_V should be one element shorter. + -- Function: int gsl_linalg_bidiag_decomp (gsl_matrix * A, gsl_vector * + TAU_U, gsl_vector * TAU_V) + This function factorizes the M-by-N matrix A into bidiagonal form U + B V^T. The diagonal and superdiagonal of the matrix B are stored + in the diagonal and superdiagonal of A. The orthogonal matrices U + and V are stored as compressed Householder vectors in the remaining + elements of A. The Householder coefficients are stored in the + vectors TAU_U and TAU_V. The length of TAU_U must equal the number + of elements in the diagonal of A and the length of TAU_V should be + one element shorter. -- Function: int gsl_linalg_bidiag_unpack (const gsl_matrix * A, const gsl_vector * TAU_U, gsl_matrix * U, const gsl_vector * TAU_V, gsl_matrix * V, gsl_vector * DIAG, gsl_vector * SUPERDIAG) This function unpacks the bidiagonal decomposition of A given by - `gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate + 'gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate orthogonal matrices U, V and the diagonal vector DIAG and - superdiagonal SUPERDIAG. Note that U is stored as a compact - M-by-N orthogonal matrix satisfying U^T U = I for efficiency. + superdiagonal SUPERDIAG. Note that U is stored as a compact M-by-N + orthogonal matrix satisfying U^T U = I for efficiency. -- Function: int gsl_linalg_bidiag_unpack2 (gsl_matrix * A, gsl_vector * TAU_U, gsl_vector * TAU_V, gsl_matrix * V) This function unpacks the bidiagonal decomposition of A given by - `gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate + 'gsl_linalg_bidiag_decomp', (A, TAU_U, TAU_V) into the separate orthogonal matrices U, V and the diagonal vector DIAG and superdiagonal SUPERDIAG. The matrix U is stored in-place in A. -- Function: int gsl_linalg_bidiag_unpack_B (const gsl_matrix * A, gsl_vector * DIAG, gsl_vector * SUPERDIAG) This function unpacks the diagonal and superdiagonal of the - bidiagonal decomposition of A given by `gsl_linalg_bidiag_decomp', + bidiagonal decomposition of A given by 'gsl_linalg_bidiag_decomp', into the diagonal vector DIAG and superdiagonal vector SUPERDIAG. @@ -1593,9 +1561,9 @@ Householder matrix P takes the form, P = I - \tau v v^T where v is a vector (called the "Householder vector") and \tau = 2/(v^T -v). The functions described in this section use the rank-1 structure -of the Householder matrix to create and apply Householder -transformations efficiently. +v). The functions described in this section use the rank-1 structure of +the Householder matrix to create and apply Householder transformations +efficiently. -- Function: double gsl_linalg_householder_transform (gsl_vector * V) This function prepares a Householder transformation P = I - \tau v @@ -1612,8 +1580,8 @@ transformations efficiently. -- Function: int gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A) This function applies the Householder matrix P defined by the - scalar TAU and the vector V to the right-hand side of the matrix - A. On output the result A P is stored in A. + scalar TAU and the vector V to the right-hand side of the matrix A. + On output the result A P is stored in A. -- Function: int gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w) @@ -1630,8 +1598,8 @@ File: gsl-ref.info, Node: Householder solver for linear systems, Next: Tridiag -- Function: int gsl_linalg_HH_solve (gsl_matrix * A, const gsl_vector * B, gsl_vector * X) This function solves the system A x = b directly using Householder - transformations. On output the solution is stored in X and B is - not modified. The matrix A is destroyed by the Householder + transformations. On output the solution is stored in X and B is + not modified. The matrix A is destroyed by the Householder transformations. -- Function: int gsl_linalg_HH_svx (gsl_matrix * A, gsl_vector * X) @@ -1647,10 +1615,10 @@ File: gsl-ref.info, Node: Tridiagonal Systems, Next: Balancing, Prev: Househo ========================= -- Function: int gsl_linalg_solve_tridiag (const gsl_vector * DIAG, - const gsl_vector * E, const gsl_vector * F, const gsl_vector - * B, gsl_vector * X) + const gsl_vector * E, const gsl_vector * F, const gsl_vector * + B, gsl_vector * X) This function solves the general N-by-N system A x = b where A is - tridiagonal (N >= 2). The super-diagonal and sub-diagonal vectors + tridiagonal (N >= 2). The super-diagonal and sub-diagonal vectors E and F must be one element shorter than the diagonal vector DIAG. The form of A for the 4-by-4 case is shown below, @@ -1660,8 +1628,8 @@ File: gsl-ref.info, Node: Tridiagonal Systems, Next: Balancing, Prev: Househo ( 0 0 f_2 d_3 ) -- Function: int gsl_linalg_solve_symm_tridiag (const gsl_vector * - DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector - * X) + DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector * + X) This function solves the general N-by-N system A x = b where A is symmetric tridiagonal (N >= 2). The off-diagonal vector E must be one element shorter than the diagonal vector DIAG. The form of A @@ -1675,9 +1643,9 @@ File: gsl-ref.info, Node: Tridiagonal Systems, Next: Balancing, Prev: Househo which can cause division by zero if the matrix is not positive definite. - -- Function: int gsl_linalg_solve_cyc_tridiag (const gsl_vector * - DIAG, const gsl_vector * E, const gsl_vector * F, const - gsl_vector * B, gsl_vector * X) + -- Function: int gsl_linalg_solve_cyc_tridiag (const gsl_vector * DIAG, + const gsl_vector * E, const gsl_vector * F, const gsl_vector * + B, gsl_vector * X) This function solves the general N-by-N system A x = b where A is cyclic tridiagonal (N >= 3). The cyclic super-diagonal and sub-diagonal vectors E and F must have the same number of elements @@ -1690,8 +1658,8 @@ File: gsl-ref.info, Node: Tridiagonal Systems, Next: Balancing, Prev: Househo ( e_3 0 f_2 d_3 ) -- Function: int gsl_linalg_solve_symm_cyc_tridiag (const gsl_vector * - DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector - * X) + DIAG, const gsl_vector * E, const gsl_vector * B, gsl_vector * + X) This function solves the general N-by-N system A x = b where A is symmetric cyclic tridiagonal (N >= 3). The cyclic off-diagonal vector E must have the same number of elements as the diagonal @@ -1709,15 +1677,15 @@ File: gsl-ref.info, Node: Balancing, Next: Linear Algebra Examples, Prev: Tri =============== The process of balancing a matrix applies similarity transformations to -make the rows and columns have comparable norms. This is useful, for +make the rows and columns have comparable norms. This is useful, for example, to reduce roundoff errors in the solution of eigenvalue -problems. Balancing a matrix A consists of replacing A with a similar +problems. Balancing a matrix A consists of replacing A with a similar matrix A' = D^(-1) A D - where D is a diagonal matrix whose entries are powers of the -floating point radix. + where D is a diagonal matrix whose entries are powers of the floating +point radix. -- Function: int gsl_linalg_balance_matrix (gsl_matrix * A, gsl_vector * D) @@ -1731,7 +1699,7 @@ File: gsl-ref.info, Node: Linear Algebra Examples, Next: Linear Algebra Refere 13.14 Examples ============== -The following program solves the linear system A x = b. The system to +The following program solves the linear system A x = b. The system to be solved is, [ 0.18 0.60 0.57 0.96 ] [x0] [1.0] @@ -1814,15 +1782,15 @@ File: gsl-ref.info, Node: Linear Algebra References and Further Reading, Prev: Further information on the algorithms described in this section can be found in the following book, - G. H. Golub, C. F. Van Loan, `Matrix Computations' (3rd Ed, 1996), + G. H. Golub, C. F. Van Loan, 'Matrix Computations' (3rd Ed, 1996), Johns Hopkins University Press, ISBN 0-8018-5414-8. The LAPACK library is described in the following manual, - `LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM, + 'LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM, ISBN 0-89871-447-8. - `http://www.netlib.org/lapack' + <http://www.netlib.org/lapack> The LAPACK source code can be found at the website above, along with an online copy of the users guide. @@ -1831,20 +1799,20 @@ The Modified Golub-Reinsch algorithm is described in the following paper, T.F. Chan, "An Improved Algorithm for Computing the Singular Value - Decomposition", `ACM Transactions on Mathematical Software', 8 + Decomposition", 'ACM Transactions on Mathematical Software', 8 (1982), pp 72-83. The Jacobi algorithm for singular value decomposition is described in the following papers, J.C. Nash, "A one-sided transformation method for the singular - value decomposition and algebraic eigenproblem", `Computer + value decomposition and algebraic eigenproblem", 'Computer Journal', Volume 18, Number 1 (1973), p 74-76 James Demmel, Kresimir Veselic, "Jacobi's Method is more accurate - than QR", `Lapack Working Note 15' (LAWN-15), October 1989. - Available from netlib, `http://www.netlib.org/lapack/' in the - `lawns' or `lawnspdf' directories. + than QR", 'Lapack Working Note 15' (LAWN-15), October 1989. + Available from netlib, <http://www.netlib.org/lapack/> in the + 'lawns' or 'lawnspdf' directories. File: gsl-ref.info, Node: Eigensystems, Next: Fast Fourier Transforms, Prev: Linear Algebra, Up: Top @@ -1854,28 +1822,28 @@ File: gsl-ref.info, Node: Eigensystems, Next: Fast Fourier Transforms, Prev: This chapter describes functions for computing eigenvalues and eigenvectors of matrices. There are routines for real symmetric, real -nonsymmetric, and complex hermitian matrices. Eigenvalues can be -computed with or without eigenvectors. The hermitian matrix algorithms -used are symmetric bidiagonalization followed by QR reduction. The +nonsymmetric, and complex hermitian matrices. Eigenvalues can be +computed with or without eigenvectors. The hermitian matrix algorithms +used are symmetric bidiagonalization followed by QR reduction. The nonsymmetric algorithm is the Francis QR double-shift. These routines are intended for "small" systems where simple -algorithms are acceptable. Anyone interested in finding eigenvalues -and eigenvectors of large matrices will want to use the sophisticated -routines found in LAPACK. The Fortran version of LAPACK is recommended +algorithms are acceptable. Anyone interested in finding eigenvalues and +eigenvectors of large matrices will want to use the sophisticated +routines found in LAPACK. The Fortran version of LAPACK is recommended as the standard package for large-scale linear algebra. The functions described in this chapter are declared in the header -file `gsl_eigen.h'. +file 'gsl_eigen.h'. * Menu: -* Real Symmetric Matrices:: -* Complex Hermitian Matrices:: -* Real Nonsymmetric Matrices:: -* Sorting Eigenvalues and Eigenvectors:: -* Eigenvalue and Eigenvector Examples:: -* Eigenvalue and Eigenvector References:: +* Real Symmetric Matrices:: +* Complex Hermitian Matrices:: +* Real Nonsymmetric Matrices:: +* Sorting Eigenvalues and Eigenvectors:: +* Eigenvalue and Eigenvector Examples:: +* Eigenvalue and Eigenvector References:: File: gsl-ref.info, Node: Real Symmetric Matrices, Next: Complex Hermitian Matrices, Up: Eigensystems @@ -1897,9 +1865,9 @@ File: gsl-ref.info, Node: Real Symmetric Matrices, Next: Complex Hermitian Mat This function computes the eigenvalues of the real symmetric matrix A. Additional workspace of the appropriate size must be provided in W. The diagonal and lower triangular part of A are destroyed - during the computation, but the strict upper triangular part is - not referenced. The eigenvalues are stored in the vector EVAL and - are unordered. + during the computation, but the strict upper triangular part is not + referenced. The eigenvalues are stored in the vector EVAL and are + unordered. -- Function: gsl_eigen_symmv_workspace * gsl_eigen_symmv_alloc (const size_t N) @@ -1914,8 +1882,8 @@ File: gsl-ref.info, Node: Real Symmetric Matrices, Next: Complex Hermitian Mat gsl_matrix * EVEC, gsl_eigen_symmv_workspace * W) This function computes the eigenvalues and eigenvectors of the real symmetric matrix A. Additional workspace of the appropriate size - must be provided in W. The diagonal and lower triangular part of - A are destroyed during the computation, but the strict upper + must be provided in W. The diagonal and lower triangular part of A + are destroyed during the computation, but the strict upper triangular part is not referenced. The eigenvalues are stored in the vector EVAL and are unordered. The corresponding eigenvectors are stored in the columns of the matrix EVEC. For example, the @@ -1945,14 +1913,14 @@ File: gsl-ref.info, Node: Complex Hermitian Matrices, Next: Real Nonsymmetric provided in W. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The imaginary parts of the diagonal are - assumed to be zero and are not referenced. The eigenvalues are + assumed to be zero and are not referenced. The eigenvalues are stored in the vector EVAL and are unordered. -- Function: gsl_eigen_hermv_workspace * gsl_eigen_hermv_alloc (const size_t N) This function allocates a workspace for computing eigenvalues and - eigenvectors of N-by-N complex hermitian matrices. The size of - the workspace is O(5n). + eigenvectors of N-by-N complex hermitian matrices. The size of the + workspace is O(5n). -- Function: void gsl_eigen_hermv_free (gsl_eigen_hermv_workspace * W) This function frees the memory associated with the workspace W. @@ -1964,7 +1932,7 @@ File: gsl-ref.info, Node: Complex Hermitian Matrices, Next: Real Nonsymmetric complex hermitian matrix A. Additional workspace of the appropriate size must be provided in W. The diagonal and lower triangular part of A are destroyed during the computation, but the - strict upper triangular part is not referenced. The imaginary + strict upper triangular part is not referenced. The imaginary parts of the diagonal are assumed to be zero and are not referenced. The eigenvalues are stored in the vector EVAL and are unordered. The corresponding complex eigenvectors are stored in @@ -1979,106 +1947,106 @@ File: gsl-ref.info, Node: Real Nonsymmetric Matrices, Next: Sorting Eigenvalue 14.3 Real Nonsymmetric Matrices =============================== -The solution of the real nonsymmetric eigensystem problem for a matrix -A involves computing the Schur decomposition +The solution of the real nonsymmetric eigensystem problem for a matrix A +involves computing the Schur decomposition A = Z T Z^T where Z is an orthogonal matrix of Schur vectors and T, the Schur form, is quasi upper triangular with diagonal 1-by-1 blocks which are real eigenvalues of A, and diagonal 2-by-2 blocks whose eigenvalues are -complex conjugate eigenvalues of A. The algorithm used is the double +complex conjugate eigenvalues of A. The algorithm used is the double shift Francis method. -- Function: gsl_eigen_nonsymm_workspace * gsl_eigen_nonsymm_alloc (const size_t N) This function allocates a workspace for computing eigenvalues of - N-by-N real nonsymmetric matrices. The size of the workspace is + N-by-N real nonsymmetric matrices. The size of the workspace is O(2n). - -- Function: void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace - * W) + -- Function: void gsl_eigen_nonsymm_free (gsl_eigen_nonsymm_workspace * + W) This function frees the memory associated with the workspace W. -- Function: void gsl_eigen_nonsymm_params (const int COMPUTE_T, const int BALANCE, gsl_eigen_nonsymm_workspace * W) This function sets some parameters which determine how the eigenvalue problem is solved in subsequent calls to - `gsl_eigen_nonsymm'. + 'gsl_eigen_nonsymm'. - If COMPUTE_T is set to 1, the full Schur form T will be computed - by `gsl_eigen_nonsymm'. If it is set to 0, T will not be computed - (this is the default setting). Computing the full Schur form T + If COMPUTE_T is set to 1, the full Schur form T will be computed by + 'gsl_eigen_nonsymm'. If it is set to 0, T will not be computed + (this is the default setting). Computing the full Schur form T requires approximately 1.5-2 times the number of flops. If BALANCE is set to 1, a balancing transformation is applied to - the matrix prior to computing eigenvalues. This transformation is + the matrix prior to computing eigenvalues. This transformation is designed to make the rows and columns of the matrix have comparable norms, and can result in more accurate eigenvalues for matrices - whose entries vary widely in magnitude. See *Note Balancing:: for - more information. Note that the balancing transformation does not + whose entries vary widely in magnitude. See *note Balancing:: for + more information. Note that the balancing transformation does not preserve the orthogonality of the Schur vectors, so if you wish to - compute the Schur vectors with `gsl_eigen_nonsymm_Z' you will + compute the Schur vectors with 'gsl_eigen_nonsymm_Z' you will obtain the Schur vectors of the balanced matrix instead of the - original matrix. The relationship will be + original matrix. The relationship will be T = Q^t D^(-1) A D Q where Q is the matrix of Schur vectors for the balanced matrix, and - D is the balancing transformation. Then `gsl_eigen_nonsymm_Z' will + D is the balancing transformation. Then 'gsl_eigen_nonsymm_Z' will compute a matrix Z which satisfies T = Z^(-1) A Z - with Z = D Q. Note that Z will not be orthogonal. For this reason, - balancing is not performed by default. + with Z = D Q. Note that Z will not be orthogonal. For this + reason, balancing is not performed by default. -- Function: int gsl_eigen_nonsymm (gsl_matrix * A, gsl_vector_complex * EVAL, gsl_eigen_nonsymm_workspace * W) This function computes the eigenvalues of the real nonsymmetric - matrix A and stores them in the vector EVAL. If T is desired, it + matrix A and stores them in the vector EVAL. If T is desired, it is stored in the upper portion of A on output. Otherwise, on output, the diagonal of A will contain the 1-by-1 real eigenvalues and 2-by-2 complex conjugate eigenvalue systems, and the rest of A - is destroyed. In rare cases, this function will fail to find all - eigenvalues. If this happens, an error code is returned and the - number of converged eigenvalues is stored in `w->n_evals'. The + is destroyed. In rare cases, this function will fail to find all + eigenvalues. If this happens, an error code is returned and the + number of converged eigenvalues is stored in 'w->n_evals'. The converged eigenvalues are stored in the beginning of EVAL. -- Function: int gsl_eigen_nonsymm_Z (gsl_matrix * A, gsl_vector_complex * EVAL, gsl_matrix * Z, gsl_eigen_nonsymm_workspace * W) - This function is identical to `gsl_eigen_nonsymm' except it also + This function is identical to 'gsl_eigen_nonsymm' except it also computes the Schur vectors and stores them into Z. -- Function: gsl_eigen_nonsymmv_workspace * gsl_eigen_nonsymmv_alloc (const size_t N) This function allocates a workspace for computing eigenvalues and - eigenvectors of N-by-N real nonsymmetric matrices. The size of the + eigenvectors of N-by-N real nonsymmetric matrices. The size of the workspace is O(5n). - -- Function: void gsl_eigen_nonsymmv_free - (gsl_eigen_nonsymmv_workspace * W) + -- Function: void gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace + * W) This function frees the memory associated with the workspace W. - -- Function: int gsl_eigen_nonsymmv (gsl_matrix * A, - gsl_vector_complex * EVAL, gsl_matrix_complex * EVEC, + -- Function: int gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex + * EVAL, gsl_matrix_complex * EVEC, gsl_eigen_nonsymmv_workspace * W) This function computes eigenvalues and right eigenvectors of the - N-by-N real nonsymmetric matrix A. It first calls - `gsl_eigen_nonsymm' to compute the eigenvalues, Schur form T, and - Schur vectors. Then it finds eigenvectors of T and backtransforms - them using the Schur vectors. The Schur vectors are destroyed in - the process, but can be saved by using `gsl_eigen_nonsymmv_Z'. The - computed eigenvectors are normalized to have Euclidean norm 1. On - output, the upper portion of A contains the Schur form T. If - `gsl_eigen_nonsymm' fails, no eigenvectors are computed, and an + N-by-N real nonsymmetric matrix A. It first calls + 'gsl_eigen_nonsymm' to compute the eigenvalues, Schur form T, and + Schur vectors. Then it finds eigenvectors of T and backtransforms + them using the Schur vectors. The Schur vectors are destroyed in + the process, but can be saved by using 'gsl_eigen_nonsymmv_Z'. The + computed eigenvectors are normalized to have Euclidean norm 1. On + output, the upper portion of A contains the Schur form T. If + 'gsl_eigen_nonsymm' fails, no eigenvectors are computed, and an error code is returned. -- Function: int gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * EVAL, gsl_matrix_complex * EVEC, gsl_matrix * Z, gsl_eigen_nonsymmv_workspace * W) - This function is identical to `gsl_eigen_nonsymmv' except it also + This function is identical to 'gsl_eigen_nonsymmv' except it also saves the Schur vectors into Z. @@ -2094,19 +2062,15 @@ File: gsl-ref.info, Node: Sorting Eigenvalues and Eigenvectors, Next: Eigenval columns of the matrix EVEC into ascending or descending order according to the value of the parameter SORT_TYPE, - `GSL_EIGEN_SORT_VAL_ASC' + 'GSL_EIGEN_SORT_VAL_ASC' ascending order in numerical value - - `GSL_EIGEN_SORT_VAL_DESC' + 'GSL_EIGEN_SORT_VAL_DESC' descending order in numerical value - - `GSL_EIGEN_SORT_ABS_ASC' + 'GSL_EIGEN_SORT_ABS_ASC' ascending order in magnitude - - `GSL_EIGEN_SORT_ABS_DESC' + 'GSL_EIGEN_SORT_ABS_DESC' descending order in magnitude - -- Function: int gsl_eigen_hermv_sort (gsl_vector * EVAL, gsl_matrix_complex * EVEC, gsl_eigen_sort_t SORT_TYPE) This function simultaneously sorts the eigenvalues stored in the @@ -2218,7 +2182,7 @@ sign of an eigenvector is arbitrary. The following program illustrates the use of the nonsymmetric eigensolver, by computing the eigenvalues and eigenvectors of the -Vandermonde matrix V(x;i,j) = x_i^n - j with x = (-1,-2,3,4). +Vandermonde matrix V(x;i,j) = x_i^{n - j} with x = (-1,-2,3,4). #include <stdio.h> #include <gsl/gsl_math.h> @@ -2319,8 +2283,9 @@ This can be compared with the corresponding output from GNU OCTAVE, 0.35660 + 0.00000i 0.91937 + 0.00000i 0.08118 + 0.00000i + Note that the eigenvectors corresponding to the eigenvalue 5.54555 + -3.08545i are slightly different. This is because they differ by the +3.08545i are slightly different. This is because they differ by the multiplicative constant 0.9999984 + 0.0017674i which has magnitude 1. @@ -2332,15 +2297,15 @@ File: gsl-ref.info, Node: Eigenvalue and Eigenvector References, Prev: Eigenva Further information on the algorithms described in this section can be found in the following book, - G. H. Golub, C. F. Van Loan, `Matrix Computations' (3rd Ed, 1996), + G. H. Golub, C. F. Van Loan, 'Matrix Computations' (3rd Ed, 1996), Johns Hopkins University Press, ISBN 0-8018-5414-8. The LAPACK library is described in, - `LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM, + 'LAPACK Users' Guide' (Third Edition, 1999), Published by SIAM, ISBN 0-89871-447-8. - `http://www.netlib.org/lapack' + <http://www.netlib.org/lapack> The LAPACK source code can be found at the website above along with an online copy of the users guide. @@ -2358,21 +2323,21 @@ efficiency there are separate versions of the routines for real data and for complex data. The mixed-radix routines are a reimplementation of the FFTPACK library of Paul Swarztrauber. Fortran code for FFTPACK is available on Netlib (FFTPACK also includes some routines for sine and -cosine transforms but these are currently not available in GSL). For +cosine transforms but these are currently not available in GSL). For details and derivations of the underlying algorithms consult the -document `GSL FFT Algorithms' (*note FFT References and Further +document 'GSL FFT Algorithms' (*note FFT References and Further Reading::) * Menu: -* Mathematical Definitions:: -* Overview of complex data FFTs:: -* Radix-2 FFT routines for complex data:: -* Mixed-radix FFT routines for complex data:: -* Overview of real data FFTs:: -* Radix-2 FFT routines for real data:: -* Mixed-radix FFT routines for real data:: -* FFT References and Further Reading:: +* Mathematical Definitions:: +* Overview of complex data FFTs:: +* Radix-2 FFT routines for complex data:: +* Mixed-radix FFT routines for complex data:: +* Overview of real data FFTs:: +* Radix-2 FFT routines for real data:: +* Mixed-radix FFT routines for real data:: +* FFT References and Further Reading:: File: gsl-ref.info, Node: Mathematical Definitions, Next: Overview of complex data FFTs, Up: Fast Fourier Transforms @@ -2388,14 +2353,13 @@ discrete fourier transform (DFT), The DFT usually arises as an approximation to the continuous fourier transform when functions are sampled at discrete intervals in space or time. The naive evaluation of the discrete fourier transform is a -matrix-vector multiplication W\vec{z}. A general matrix-vector +matrix-vector multiplication W\vec{z}. A general matrix-vector multiplication takes O(N^2) operations for N data-points. Fast fourier transform algorithms use a divide-and-conquer strategy to factorize the -matrix W into smaller sub-matrices, corresponding to the integer -factors of the length N. If N can be factorized into a product of -integers f_1 f_2 ... f_n then the DFT can be computed in O(N \sum f_i) -operations. For a radix-2 FFT this gives an operation count of O(N -\log_2 N). +matrix W into smaller sub-matrices, corresponding to the integer factors +of the length N. If N can be factorized into a product of integers f_1 +f_2 ... f_n then the DFT can be computed in O(N \sum f_i) operations. +For a radix-2 FFT this gives an operation count of O(N \log_2 N). All the FFT functions offer three types of transform: forwards, inverse and backwards, based on the same mathematical definitions. The @@ -2408,20 +2372,20 @@ and the definition of the "inverse fourier transform", x = IFFT(z), is, z_j = {1 \over N} \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N). The factor of 1/N makes this a true inverse. For example, a call to -`gsl_fft_complex_forward' followed by a call to -`gsl_fft_complex_inverse' should return the original data (within +'gsl_fft_complex_forward' followed by a call to +'gsl_fft_complex_inverse' should return the original data (within numerical errors). In general there are two possible choices for the sign of the -exponential in the transform/ inverse-transform pair. GSL follows the +exponential in the transform/ inverse-transform pair. GSL follows the same convention as FFTPACK, using a negative exponential for the forward transform. The advantage of this convention is that the inverse -transform recreates the original function with simple fourier -synthesis. Numerical Recipes uses the opposite convention, a positive -exponential in the forward transform. +transform recreates the original function with simple fourier synthesis. +Numerical Recipes uses the opposite convention, a positive exponential +in the forward transform. - The "backwards FFT" is simply our terminology for an unscaled -version of the inverse FFT, + The "backwards FFT" is simply our terminology for an unscaled version +of the inverse FFT, z^{backwards}_j = \sum_{k=0}^{N-1} x_k \exp(2\pi i j k / N). @@ -2444,7 +2408,7 @@ length 6, double x[3*2]; gsl_complex_packed_array data = x; -can be used to hold an array of three complex numbers, `z[3]', in the +can be used to hold an array of three complex numbers, 'z[3]', in the following way, data[0] = Re(z[0]) @@ -2455,16 +2419,16 @@ following way, data[5] = Im(z[2]) The array indices for the data have the same ordering as those in the -definition of the DFT--i.e. there are no index transformations or +definition of the DFT--i.e. there are no index transformations or permutations of the data. A "stride" parameter allows the user to perform transforms on the -elements `z[stride*i]' instead of `z[i]'. A stride greater than 1 can -be used to take an in-place FFT of the column of a matrix. A stride of +elements 'z[stride*i]' instead of 'z[i]'. A stride greater than 1 can +be used to take an in-place FFT of the column of a matrix. A stride of 1 accesses the array without any additional spacing between elements. - To perform an FFT on a vector argument, such as `gsl_vector_complex -* v', use the following definitions (or their equivalents) when calling + To perform an FFT on a vector argument, such as 'gsl_vector_complex * +v', use the following definitions (or their equivalents) when calling the functions described in this chapter: gsl_complex_packed_array data = v->data; @@ -2476,9 +2440,9 @@ appearing in the DFT does not correspond directly to a physical frequency. If the time-step of the DFT is \Delta then the frequency-domain includes both positive and negative frequencies, ranging from -1/(2\Delta) through 0 to +1/(2\Delta). The positive -frequencies are stored from the beginning of the array up to the -middle, and the negative frequencies are stored backwards from the end -of the array. +frequencies are stored from the beginning of the array up to the middle, +and the negative frequencies are stored backwards from the end of the +array. Here is a table which shows the layout of the array DATA, and the correspondence between the time-domain data z, and the frequency-domain @@ -2511,50 +2475,57 @@ File: gsl-ref.info, Node: Radix-2 FFT routines for complex data, Next: Mixed-r The radix-2 algorithms described in this section are simple and compact, although not necessarily the most efficient. They use the Cooley-Tukey algorithm to compute in-place complex FFTs for lengths which are a power -of 2--no additional storage is required. The corresponding -self-sorting mixed-radix routines offer better performance at the -expense of requiring additional working space. +of 2--no additional storage is required. The corresponding self-sorting +mixed-radix routines offer better performance at the expense of +requiring additional working space. All the functions described in this section are declared in the -header file `gsl_fft_complex.h'. +header file 'gsl_fft_complex.h'. -- Function: int gsl_fft_complex_radix2_forward (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + -- Function: int gsl_fft_complex_radix2_transform (gsl_complex_packed_array DATA, size_t STRIDE, size_t N, gsl_fft_direction SIGN) + -- Function: int gsl_fft_complex_radix2_backward (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + -- Function: int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + These functions compute forward, backward and inverse FFTs of - length N with stride STRIDE, on the packed complex array DATA - using an in-place radix-2 decimation-in-time algorithm. The - length of the transform N is restricted to powers of two. For the - `transform' version of the function the SIGN argument can be - either `forward' (-1) or `backward' (+1). - - The functions return a value of `GSL_SUCCESS' if no errors were - detected, or `GSL_EDOM' if the length of the data N is not a power + length N with stride STRIDE, on the packed complex array DATA using + an in-place radix-2 decimation-in-time algorithm. The length of + the transform N is restricted to powers of two. For the + 'transform' version of the function the SIGN argument can be either + 'forward' (-1) or 'backward' (+1). + + The functions return a value of 'GSL_SUCCESS' if no errors were + detected, or 'GSL_EDOM' if the length of the data N is not a power of two. -- Function: int gsl_fft_complex_radix2_dif_forward (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + -- Function: int gsl_fft_complex_radix2_dif_transform (gsl_complex_packed_array DATA, size_t STRIDE, size_t N, gsl_fft_direction SIGN) + -- Function: int gsl_fft_complex_radix2_dif_backward (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + -- Function: int gsl_fft_complex_radix2_dif_inverse (gsl_complex_packed_array DATA, size_t STRIDE, size_t N) + These are decimation-in-frequency versions of the radix-2 FFT functions. - - Here is an example program which computes the FFT of a short pulse -in a sample of length 128. To make the resulting fourier transform -real the pulse is defined for equal positive and negative times (-10 -... 10), where the negative times wrap around the end of the array. + Here is an example program which computes the FFT of a short pulse in +a sample of length 128. To make the resulting fourier transform real +the pulse is defined for equal positive and negative times (-10 ... 10), +where the negative times wrap around the end of the array. #include <stdio.h> #include <math.h> @@ -2601,16 +2572,16 @@ real the pulse is defined for equal positive and negative times (-10 } Note that we have assumed that the program is using the default error -handler (which calls `abort' for any errors). If you are not using a +handler (which calls 'abort' for any errors). If you are not using a safe error handler you would need to check the return status of -`gsl_fft_complex_radix2_forward'. +'gsl_fft_complex_radix2_forward'. The transformed data is rescaled by 1/\sqrt N so that it fits on the same plot as the input. Only the real part is shown, by the choice of -the input data the imaginary part is zero. Allowing for the -wrap-around of negative times at t=128, and working in units of k/N, -the DFT approximates the continuum fourier transform, giving a -modulated sine function. +the input data the imaginary part is zero. Allowing for the wrap-around +of negative times at t=128, and working in units of k/N, the DFT +approximates the continuum fourier transform, giving a modulated sine +function. File: gsl-ref.info, Node: Mixed-radix FFT routines for complex data, Next: Overview of real data FFTs, Prev: Radix-2 FFT routines for complex data, Up: Fast Fourier Transforms @@ -2621,7 +2592,7 @@ File: gsl-ref.info, Node: Mixed-radix FFT routines for complex data, Next: Ove This section describes mixed-radix FFT algorithms for complex data. The mixed-radix functions work for FFTs of any length. They are a reimplementation of Paul Swarztrauber's Fortran FFTPACK library. The -theory is explained in the review article `Self-sorting Mixed-radix +theory is explained in the review article 'Self-sorting Mixed-radix FFTs' by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK. @@ -2631,49 +2602,48 @@ There are efficient modules for factors of 2, 3, 4, 5, 6 and 7. The modules for the composite factors of 4 and 6 are faster than combining the modules for 2*2 and 2*3. - For factors which are not implemented as modules there is a -fall-back to a general length-n module which uses Singleton's method for + For factors which are not implemented as modules there is a fall-back +to a general length-n module which uses Singleton's method for efficiently computing a DFT. This module is O(n^2), and slower than a dedicated module would be but works for any length n. Of course, lengths which use the general length-n module will still be factorized as much as possible. For example, a length of 143 will be factorized -into 11*13. Large prime factors are the worst case scenario, e.g. as -found in n=2*3*99991, and should be avoided because their O(n^2) -scaling will dominate the run-time (consult the document `GSL FFT -Algorithms' included in the GSL distribution if you encounter this -problem). +into 11*13. Large prime factors are the worst case scenario, e.g. as +found in n=2*3*99991, and should be avoided because their O(n^2) scaling +will dominate the run-time (consult the document 'GSL FFT Algorithms' +included in the GSL distribution if you encounter this problem). The mixed-radix initialization function -`gsl_fft_complex_wavetable_alloc' returns the list of factors chosen by +'gsl_fft_complex_wavetable_alloc' returns the list of factors chosen by the library for a given length N. It can be used to check how well the length has been factorized, and estimate the run-time. To a first approximation the run-time scales as N \sum f_i, where the f_i are the factors of N. For programs under user control you may wish to issue a warning that the transform will be slow when the length is poorly factorized. If you frequently encounter data lengths which cannot be -factorized using the existing small-prime modules consult `GSL FFT +factorized using the existing small-prime modules consult 'GSL FFT Algorithms' for details on adding support for other factors. All the functions described in this section are declared in the -header file `gsl_fft_complex.h'. +header file 'gsl_fft_complex.h'. -- Function: gsl_fft_complex_wavetable * -gsl_fft_complex_wavetable_alloc (size_t N) + gsl_fft_complex_wavetable_alloc (size_t N) This function prepares a trigonometric lookup table for a complex - FFT of length N. The function returns a pointer to the newly - allocated `gsl_fft_complex_wavetable' if no errors were detected, + FFT of length N. The function returns a pointer to the newly + allocated 'gsl_fft_complex_wavetable' if no errors were detected, and a null pointer in the case of error. The length N is factorized into a product of subtransforms, and the factors and - their trigonometric coefficients are stored in the wavetable. The - trigonometric coefficients are computed using direct calls to - `sin' and `cos', for accuracy. Recursion relations could be used - to compute the lookup table faster, but if an application performs + their trigonometric coefficients are stored in the wavetable. The + trigonometric coefficients are computed using direct calls to 'sin' + and 'cos', for accuracy. Recursion relations could be used to + compute the lookup table faster, but if an application performs many FFTs of the same length then this computation is a one-off overhead which does not affect the final throughput. - The wavetable structure can be used repeatedly for any transform - of the same length. The table is not modified by calls to any of - the other FFT functions. The same wavetable can be used for both + The wavetable structure can be used repeatedly for any transform of + the same length. The table is not modified by calls to any of the + other FFT functions. The same wavetable can be used for both forward and backward (or inverse) transforms of a given length. -- Function: void gsl_fft_complex_wavetable_free @@ -2682,94 +2652,92 @@ gsl_fft_complex_wavetable_alloc (size_t N) WAVETABLE. The wavetable can be freed if no further FFTs of the same length will be needed. -These functions operate on a `gsl_fft_complex_wavetable' structure -which contains internal parameters for the FFT. It is not necessary to -set any of the components directly but it can sometimes be useful to -examine them. For example, the chosen factorization of the FFT length -is given and can be used to provide an estimate of the run-time or -numerical error. The wavetable structure is declared in the header file -`gsl_fft_complex.h'. +These functions operate on a 'gsl_fft_complex_wavetable' structure which +contains internal parameters for the FFT. It is not necessary to set any +of the components directly but it can sometimes be useful to examine +them. For example, the chosen factorization of the FFT length is given +and can be used to provide an estimate of the run-time or numerical +error. The wavetable structure is declared in the header file +'gsl_fft_complex.h'. -- Data Type: gsl_fft_complex_wavetable This is a structure that holds the factorization and trigonometric lookup tables for the mixed radix fft algorithm. It has the following components: - `size_t n' + 'size_t n' This is the number of complex data points - `size_t nf' - This is the number of factors that the length `n' was + 'size_t nf' + This is the number of factors that the length 'n' was decomposed into. - `size_t factor[64]' - This is the array of factors. Only the first `nf' elements + 'size_t factor[64]' + This is the array of factors. Only the first 'nf' elements are used. - `gsl_complex * trig' - This is a pointer to a preallocated trigonometric lookup - table of `n' complex elements. + 'gsl_complex * trig' + This is a pointer to a preallocated trigonometric lookup table + of 'n' complex elements. - `gsl_complex * twiddle[64]' - This is an array of pointers into `trig', giving the twiddle + 'gsl_complex * twiddle[64]' + This is an array of pointers into 'trig', giving the twiddle factors for each pass. The mixed radix algorithms require additional working space to hold the intermediate steps of the transform. -- Function: gsl_fft_complex_workspace * -gsl_fft_complex_workspace_alloc (size_t N) + gsl_fft_complex_workspace_alloc (size_t N) This function allocates a workspace for a complex transform of length N. -- Function: void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * WORKSPACE) This function frees the memory associated with the workspace - WORKSPACE. The workspace can be freed if no further FFTs of the + WORKSPACE. The workspace can be freed if no further FFTs of the same length will be needed. The following functions compute the transform, -- Function: int gsl_fft_complex_forward (gsl_complex_packed_array - DATA, size_t STRIDE, size_t N, const - gsl_fft_complex_wavetable * WAVETABLE, - gsl_fft_complex_workspace * WORK) + DATA, size_t STRIDE, size_t N, const gsl_fft_complex_wavetable + * WAVETABLE, gsl_fft_complex_workspace * WORK) -- Function: int gsl_fft_complex_transform (gsl_complex_packed_array - DATA, size_t STRIDE, size_t N, const - gsl_fft_complex_wavetable * WAVETABLE, - gsl_fft_complex_workspace * WORK, gsl_fft_direction SIGN) + DATA, size_t STRIDE, size_t N, const gsl_fft_complex_wavetable + * WAVETABLE, gsl_fft_complex_workspace * WORK, + gsl_fft_direction SIGN) -- Function: int gsl_fft_complex_backward (gsl_complex_packed_array - DATA, size_t STRIDE, size_t N, const - gsl_fft_complex_wavetable * WAVETABLE, - gsl_fft_complex_workspace * WORK) + DATA, size_t STRIDE, size_t N, const gsl_fft_complex_wavetable + * WAVETABLE, gsl_fft_complex_workspace * WORK) -- Function: int gsl_fft_complex_inverse (gsl_complex_packed_array - DATA, size_t STRIDE, size_t N, const - gsl_fft_complex_wavetable * WAVETABLE, - gsl_fft_complex_workspace * WORK) + DATA, size_t STRIDE, size_t N, const gsl_fft_complex_wavetable + * WAVETABLE, gsl_fft_complex_workspace * WORK) + These functions compute forward, backward and inverse FFTs of length N with stride STRIDE, on the packed complex array DATA, - using a mixed radix decimation-in-frequency algorithm. There is - no restriction on the length N. Efficient modules are provided - for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining - factors are computed with a slow, O(n^2), general-n module. The - caller must supply a WAVETABLE containing the trigonometric lookup - tables and a workspace WORK. For the `transform' version of the - function the SIGN argument can be either `forward' (-1) or - `backward' (+1). - - The functions return a value of `0' if no errors were detected. The - following `gsl_errno' conditions are defined for these functions: - - `GSL_EDOM' - The length of the data N is not a positive integer (i.e. N is + using a mixed radix decimation-in-frequency algorithm. There is no + restriction on the length N. Efficient modules are provided for + subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining factors + are computed with a slow, O(n^2), general-n module. The caller + must supply a WAVETABLE containing the trigonometric lookup tables + and a workspace WORK. For the 'transform' version of the function + the SIGN argument can be either 'forward' (-1) or 'backward' (+1). + + The functions return a value of '0' if no errors were detected. + The following 'gsl_errno' conditions are defined for these + functions: + + 'GSL_EDOM' + The length of the data N is not a positive integer (i.e. N is zero). - `GSL_EINVAL' + 'GSL_EINVAL' The length of the data N and the length used to compute the given WAVETABLE do not match. - Here is an example program which computes the FFT of a short pulse -in a sample of length 630 (=2*3*3*5*7) using the mixed-radix algorithm. + Here is an example program which computes the FFT of a short pulse in +a sample of length 630 (=2*3*3*5*7) using the mixed-radix algorithm. #include <stdio.h> #include <math.h> @@ -2832,10 +2800,10 @@ in a sample of length 630 (=2*3*3*5*7) using the mixed-radix algorithm. return 0; } -Note that we have assumed that the program is using the default `gsl' -error handler (which calls `abort' for any errors). If you are not +Note that we have assumed that the program is using the default 'gsl' +error handler (which calls 'abort' for any errors). If you are not using a safe error handler you would need to check the return status of -all the `gsl' routines. +all the 'gsl' routines. File: gsl-ref.info, Node: Overview of real data FFTs, Next: Radix-2 FFT routines for real data, Prev: Mixed-radix FFT routines for complex data, Up: Fast Fourier Transforms @@ -2854,17 +2822,17 @@ A sequence with this symmetry is called "conjugate-complex" or "half-complex". This different structure requires different storage layouts for the forward transform (from real to half-complex) and inverse transform (from half-complex back to real). As a consequence -the routines are divided into two sets: functions in `gsl_fft_real' -which operate on real sequences and functions in `gsl_fft_halfcomplex' +the routines are divided into two sets: functions in 'gsl_fft_real' +which operate on real sequences and functions in 'gsl_fft_halfcomplex' which operate on half-complex sequences. - Functions in `gsl_fft_real' compute the frequency coefficients of a -real sequence. The half-complex coefficients c of a real sequence x -are given by fourier analysis, + Functions in 'gsl_fft_real' compute the frequency coefficients of a +real sequence. The half-complex coefficients c of a real sequence x are +given by fourier analysis, c_k = \sum_{j=0}^{N-1} x_j \exp(-2 \pi i j k /N) -Functions in `gsl_fft_halfcomplex' compute inverse or backwards +Functions in 'gsl_fft_halfcomplex' compute inverse or backwards transforms. They reconstruct real sequences by fourier synthesis from their half-complex frequency coefficients, c, @@ -2872,7 +2840,7 @@ their half-complex frequency coefficients, c, The symmetry of the half-complex sequence implies that only half of the complex numbers in the output need to be stored. The remaining half can -be reconstructed using the half-complex symmetry condition. This works +be reconstructed using the half-complex symmetry condition. This works for all lengths, even and odd--when the length is even the middle value where k=N/2 is also real. Thus only N real numbers are required to store the half-complex sequence, and the transform of a real sequence @@ -2898,21 +2866,22 @@ the Cooley-Tukey algorithm to compute in-place FFTs for lengths which are a power of 2. The radix-2 FFT functions for real data are declared in the header -files `gsl_fft_real.h' +files 'gsl_fft_real.h' -- Function: int gsl_fft_real_radix2_transform (double DATA[], size_t STRIDE, size_t N) + This function computes an in-place radix-2 FFT of length N and - stride STRIDE on the real array DATA. The output is a - half-complex sequence, which is stored in-place. The arrangement - of the half-complex terms uses the following scheme: for k < N/2 - the real part of the k-th term is stored in location k, and the - corresponding imaginary part is stored in location N-k. Terms - with k > N/2 can be reconstructed using the symmetry z_k = - z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and - count as a special case. Their real parts are stored in locations - 0 and N/2 respectively, while their imaginary parts which are zero - are not stored. + stride STRIDE on the real array DATA. The output is a half-complex + sequence, which is stored in-place. The arrangement of the + half-complex terms uses the following scheme: for k < N/2 the real + part of the k-th term is stored in location k, and the + corresponding imaginary part is stored in location N-k. Terms with + k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. + The terms for k=0 and k=N/2 are both purely real, and count as a + special case. Their real parts are stored in locations 0 and N/2 + respectively, while their imaginary parts which are zero are not + stored. The following table shows the correspondence between the output DATA and the equivalent results obtained by considering the input @@ -2935,22 +2904,21 @@ files `gsl_fft_real.h' complex[N-1].real = data[1] complex[N-1].imag = -data[N-1] Note that the output data can be converted into the full complex - sequence using the function `gsl_fft_halfcomplex_unpack' described + sequence using the function 'gsl_fft_halfcomplex_unpack' described in the next section. The radix-2 FFT functions for halfcomplex data are declared in the -header file `gsl_fft_halfcomplex.h'. +header file 'gsl_fft_halfcomplex.h'. -- Function: int gsl_fft_halfcomplex_radix2_inverse (double DATA[], size_t STRIDE, size_t N) -- Function: int gsl_fft_halfcomplex_radix2_backward (double DATA[], size_t STRIDE, size_t N) - These functions compute the inverse or backwards in-place radix-2 - FFT of length N and stride STRIDE on the half-complex sequence - DATA stored according the output scheme used by - `gsl_fft_real_radix2'. The result is a real array stored in - natural order. + These functions compute the inverse or backwards in-place radix-2 + FFT of length N and stride STRIDE on the half-complex sequence DATA + stored according the output scheme used by 'gsl_fft_real_radix2'. + The result is a real array stored in natural order. File: gsl-ref.info, Node: Mixed-radix FFT routines for real data, Next: FFT References and Further Reading, Prev: Radix-2 FFT routines for real data, Up: Fast Fourier Transforms @@ -2962,7 +2930,7 @@ This section describes mixed-radix FFT algorithms for real data. The mixed-radix functions work for FFTs of any length. They are a reimplementation of the real-FFT routines in the Fortran FFTPACK library by Paul Swarztrauber. The theory behind the algorithm is explained in -the article `Fast Mixed-Radix Real Fourier Transforms' by Clive +the article 'Fast Mixed-Radix Real Fourier Transforms' by Clive Temperton. The routines here use the same indexing scheme and basic algorithms as FFTPACK. @@ -2980,10 +2948,10 @@ since the symmetry z_k = z_{N-k}^* implies that this is purely real too. The storage scheme is best shown by some examples. The table below shows the output for an odd-length sequence, n=5. The two columns give the correspondence between the 5 values in the half-complex sequence -returned by `gsl_fft_real_transform', HALFCOMPLEX[] and the values +returned by 'gsl_fft_real_transform', HALFCOMPLEX[] and the values COMPLEX[] that would be returned if the same real input sequence were -passed to `gsl_fft_complex_backward' as a complex sequence (with -imaginary parts set to `0'), +passed to 'gsl_fft_complex_backward' as a complex sequence (with +imaginary parts set to '0'), complex[0].real = halfcomplex[0] complex[0].imag = 0 @@ -2996,9 +2964,9 @@ imaginary parts set to `0'), complex[4].real = halfcomplex[1] complex[4].imag = -halfcomplex[2] -The upper elements of the COMPLEX array, `complex[3]' and `complex[4]' +The upper elements of the COMPLEX array, 'complex[3]' and 'complex[4]' are filled in using the symmetry condition. The imaginary part of the -zero-frequency term `complex[0].imag' is known to be zero by the +zero-frequency term 'complex[0].imag' is known to be zero by the symmetry. The next table shows the output for an even-length sequence, n=6 In @@ -3017,40 +2985,40 @@ the even case there are two values which are purely real, complex[5].real = halfcomplex[1] complex[5].imag = -halfcomplex[2] -The upper elements of the COMPLEX array, `complex[4]' and `complex[5]' -are filled in using the symmetry condition. Both `complex[0].imag' and -`complex[3].imag' are known to be zero. +The upper elements of the COMPLEX array, 'complex[4]' and 'complex[5]' +are filled in using the symmetry condition. Both 'complex[0].imag' and +'complex[3].imag' are known to be zero. - All these functions are declared in the header files -`gsl_fft_real.h' and `gsl_fft_halfcomplex.h'. + All these functions are declared in the header files 'gsl_fft_real.h' +and 'gsl_fft_halfcomplex.h'. -- Function: gsl_fft_real_wavetable * gsl_fft_real_wavetable_alloc (size_t N) -- Function: gsl_fft_halfcomplex_wavetable * -gsl_fft_halfcomplex_wavetable_alloc (size_t N) + gsl_fft_halfcomplex_wavetable_alloc (size_t N) These functions prepare trigonometric lookup tables for an FFT of size n real elements. The functions return a pointer to the newly allocated struct if no errors were detected, and a null pointer in the case of error. The length N is factorized into a product of - subtransforms, and the factors and their trigonometric - coefficients are stored in the wavetable. The trigonometric - coefficients are computed using direct calls to `sin' and `cos', - for accuracy. Recursion relations could be used to compute the - lookup table faster, but if an application performs many FFTs of - the same length then computing the wavetable is a one-off overhead - which does not affect the final throughput. - - The wavetable structure can be used repeatedly for any transform - of the same length. The table is not modified by calls to any of - the other FFT functions. The appropriate type of wavetable must - be used for forward real or inverse half-complex transforms. - - -- Function: void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable - * WAVETABLE) + subtransforms, and the factors and their trigonometric coefficients + are stored in the wavetable. The trigonometric coefficients are + computed using direct calls to 'sin' and 'cos', for accuracy. + Recursion relations could be used to compute the lookup table + faster, but if an application performs many FFTs of the same length + then computing the wavetable is a one-off overhead which does not + affect the final throughput. + + The wavetable structure can be used repeatedly for any transform of + the same length. The table is not modified by calls to any of the + other FFT functions. The appropriate type of wavetable must be + used for forward real or inverse half-complex transforms. + + -- Function: void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * + WAVETABLE) -- Function: void gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * WAVETABLE) These functions free the memory associated with the wavetable - WAVETABLE. The wavetable can be freed if no further FFTs of the + WAVETABLE. The wavetable can be freed if no further FFTs of the same length will be needed. The mixed radix algorithms require additional working space to hold the @@ -3062,10 +3030,10 @@ intermediate steps of the transform, N. The same workspace can be used for both forward real and inverse halfcomplex transforms. - -- Function: void gsl_fft_real_workspace_free (gsl_fft_real_workspace - * WORKSPACE) + -- Function: void gsl_fft_real_workspace_free (gsl_fft_real_workspace * + WORKSPACE) This function frees the memory associated with the workspace - WORKSPACE. The workspace can be freed if no further FFTs of the + WORKSPACE. The workspace can be freed if no further FFTs of the same length will be needed. The following functions compute the transforms of real and half-complex @@ -3079,8 +3047,8 @@ data, WAVETABLE, gsl_fft_real_workspace * WORK) These functions compute the FFT of DATA, a real or half-complex array of length N, using a mixed radix decimation-in-frequency - algorithm. For `gsl_fft_real_transform' DATA is an array of - time-ordered real data. For `gsl_fft_halfcomplex_transform' DATA + algorithm. For 'gsl_fft_real_transform' DATA is an array of + time-ordered real data. For 'gsl_fft_halfcomplex_transform' DATA contains fourier coefficients in the half-complex ordering described above. There is no restriction on the length N. Efficient modules are provided for subtransforms of length 2, 3, 4 @@ -3091,9 +3059,10 @@ data, -- Function: int gsl_fft_real_unpack (const double REAL_COEFFICIENT[], gsl_complex_packed_array COMPLEX_COEFFICIENT, size_t STRIDE, size_t N) + This function converts a single real array, REAL_COEFFICIENT into an equivalent complex array, COMPLEX_COEFFICIENT, (with imaginary - part set to zero), suitable for `gsl_fft_complex' routines. The + part set to zero), suitable for 'gsl_fft_complex' routines. The algorithm for the conversion is simply, for (i = 0; i < n; i++) @@ -3107,12 +3076,12 @@ data, -- Function: int gsl_fft_halfcomplex_unpack (const double HALFCOMPLEX_COEFFICIENT[], gsl_complex_packed_array COMPLEX_COEFFICIENT, size_t STRIDE, size_t N) + This function converts HALFCOMPLEX_COEFFICIENT, an array of - half-complex coefficients as returned by `gsl_fft_real_transform', + half-complex coefficients as returned by 'gsl_fft_real_transform', into an ordinary complex array, COMPLEX_COEFFICIENT. It fills in - the complex array using the symmetry z_k = z_{N-k}^* to - reconstruct the redundant elements. The algorithm for the - conversion is, + the complex array using the symmetry z_k = z_{N-k}^* to reconstruct + the redundant elements. The algorithm for the conversion is, complex_coefficient[0].real = halfcomplex_coefficient[0]; @@ -3139,11 +3108,11 @@ data, = 0.0; } - Here is an example program using `gsl_fft_real_transform' and -`gsl_fft_halfcomplex_inverse'. It generates a real signal in the shape -of a square pulse. The pulse is fourier transformed to frequency -space, and all but the lowest ten frequency components are removed from -the array of fourier coefficients returned by `gsl_fft_real_transform'. + Here is an example program using 'gsl_fft_real_transform' and +'gsl_fft_halfcomplex_inverse'. It generates a real signal in the shape +of a square pulse. The pulse is fourier transformed to frequency space, +and all but the lowest ten frequency components are removed from the +array of fourier coefficients returned by 'gsl_fft_real_transform'. The remaining fourier coefficients are transformed back to the time-domain, to give a filtered version of the square pulse. Since @@ -3218,66 +3187,66 @@ File: gsl-ref.info, Node: FFT References and Further Reading, Prev: Mixed-radi =================================== A good starting point for learning more about the FFT is the review -article `Fast Fourier Transforms: A Tutorial Review and A State of the +article 'Fast Fourier Transforms: A Tutorial Review and A State of the Art' by Duhamel and Vetterli, P. Duhamel and M. Vetterli. Fast fourier transforms: A tutorial - review and a state of the art. `Signal Processing', 19:259-299, + review and a state of the art. 'Signal Processing', 19:259-299, 1990. To find out about the algorithms used in the GSL routines you may want -to consult the document `GSL FFT Algorithms' (it is included in GSL, as -`doc/fftalgorithms.tex'). This has general information on FFTs and +to consult the document 'GSL FFT Algorithms' (it is included in GSL, as +'doc/fftalgorithms.tex'). This has general information on FFTs and explicit derivations of the implementation for each routine. There are -also references to the relevant literature. For convenience some of -the more important references are reproduced below. +also references to the relevant literature. For convenience some of the +more important references are reproduced below. There are several introductory books on the FFT with example programs, -such as `The Fast Fourier Transform' by Brigham and `DFT/FFT and +such as 'The Fast Fourier Transform' by Brigham and 'DFT/FFT and Convolution Algorithms' by Burrus and Parks, - E. Oran Brigham. `The Fast Fourier Transform'. Prentice Hall, + E. Oran Brigham. 'The Fast Fourier Transform'. Prentice Hall, 1974. - C. S. Burrus and T. W. Parks. `DFT/FFT and Convolution + C. S. Burrus and T. W. Parks. 'DFT/FFT and Convolution Algorithms'. Wiley, 1984. -Both these introductory books cover the radix-2 FFT in some detail. -The mixed-radix algorithm at the heart of the FFTPACK routines is -reviewed in Clive Temperton's paper, +Both these introductory books cover the radix-2 FFT in some detail. The +mixed-radix algorithm at the heart of the FFTPACK routines is reviewed +in Clive Temperton's paper, Clive Temperton. Self-sorting mixed-radix fast fourier transforms. - `Journal of Computational Physics', 52(1):1-23, 1983. + 'Journal of Computational Physics', 52(1):1-23, 1983. The derivation of FFTs for real-valued data is explained in the following two articles, Henrik V. Sorenson, Douglas L. Jones, Michael T. Heideman, and C. Sidney Burrus. Real-valued fast fourier transform algorithms. - `IEEE Transactions on Acoustics, Speech, and Signal Processing', + 'IEEE Transactions on Acoustics, Speech, and Signal Processing', ASSP-35(6):849-863, 1987. Clive Temperton. Fast mixed-radix real fourier transforms. - `Journal of Computational Physics', 52:340-350, 1983. + 'Journal of Computational Physics', 52:340-350, 1983. In 1979 the IEEE published a compendium of carefully-reviewed Fortran -FFT programs in `Programs for Digital Signal Processing'. It is a +FFT programs in 'Programs for Digital Signal Processing'. It is a useful reference for implementations of many different FFT algorithms, - Digital Signal Processing Committee and IEEE Acoustics, Speech, - and Signal Processing Committee, editors. `Programs for Digital - Signal Processing'. IEEE Press, 1979. + Digital Signal Processing Committee and IEEE Acoustics, Speech, and + Signal Processing Committee, editors. 'Programs for Digital Signal + Processing'. IEEE Press, 1979. For large-scale FFT work we recommend the use of the dedicated FFTW library by Frigo and Johnson. The FFTW library is self-optimizing--it automatically tunes itself for each hardware platform in order to achieve maximum performance. It is available under the GNU GPL. - FFTW Website, `http://www.fftw.org/' + FFTW Website, <http://www.fftw.org/> The source code for FFTPACK is available from Netlib, - FFTPACK, `http://www.netlib.org/fftpack/' + FFTPACK, <http://www.netlib.org/fftpack/> File: gsl-ref.info, Node: Numerical Integration, Next: Random Number Generation, Prev: Fast Fourier Transforms, Up: Top @@ -3297,23 +3266,23 @@ Doncker-Kapenga, Uberhuber and Kahaner. Fortran code for QUADPACK is available on Netlib. The functions described in this chapter are declared in the header -file `gsl_integration.h'. +file 'gsl_integration.h'. * Menu: -* Numerical Integration Introduction:: -* QNG non-adaptive Gauss-Kronrod integration:: -* QAG adaptive integration:: -* QAGS adaptive integration with singularities:: -* QAGP adaptive integration with known singular points:: -* QAGI adaptive integration on infinite intervals:: -* QAWC adaptive integration for Cauchy principal values:: -* QAWS adaptive integration for singular functions:: -* QAWO adaptive integration for oscillatory functions:: -* QAWF adaptive integration for Fourier integrals:: -* Numerical integration error codes:: -* Numerical integration examples:: -* Numerical integration References and Further Reading:: +* Numerical Integration Introduction:: +* QNG non-adaptive Gauss-Kronrod integration:: +* QAG adaptive integration:: +* QAGS adaptive integration with singularities:: +* QAGP adaptive integration with known singular points:: +* QAGI adaptive integration on infinite intervals:: +* QAWC adaptive integration for Cauchy principal values:: +* QAWS adaptive integration for singular functions:: +* QAWO adaptive integration for oscillatory functions:: +* QAWF adaptive integration for Fourier integrals:: +* Numerical integration error codes:: +* Numerical integration examples:: +* Numerical integration References and Further Reading:: File: gsl-ref.info, Node: Numerical Integration Introduction, Next: QNG non-adaptive Gauss-Kronrod integration, Up: Numerical Integration @@ -3333,8 +3302,8 @@ specify the following accuracy requirement, |RESULT - I| <= max(epsabs, epsrel |I|) where RESULT is the numerical approximation obtained by the algorithm. -The algorithms attempt to estimate the absolute error ABSERR = |RESULT -- I| in such a way that the following inequality holds, +The algorithms attempt to estimate the absolute error ABSERR = |RESULT - +I| in such a way that the following inequality holds, |RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|) @@ -3345,20 +3314,20 @@ stage. The algorithms in QUADPACK use a naming convention based on the following letters, - `Q' - quadrature routine + 'Q' - quadrature routine - `N' - non-adaptive integrator - `A' - adaptive integrator + 'N' - non-adaptive integrator + 'A' - adaptive integrator - `G' - general integrand (user-defined) - `W' - weight function with integrand + 'G' - general integrand (user-defined) + 'W' - weight function with integrand - `S' - singularities can be more readily integrated - `P' - points of special difficulty can be supplied - `I' - infinite range of integration - `O' - oscillatory weight function, cos or sin - `F' - Fourier integral - `C' - Cauchy principal value + 'S' - singularities can be more readily integrated + 'P' - points of special difficulty can be supplied + 'I' - infinite range of integration + 'O' - oscillatory weight function, cos or sin + 'F' - Fourier integral + 'C' - Cauchy principal value The algorithms are built on pairs of quadrature rules, a higher order rule and a lower order rule. The higher order rule is used to compute @@ -3368,9 +3337,9 @@ order rule gives an estimate of the error in the approximation. * Menu: -* Integrands without weight functions:: -* Integrands with weight functions:: -* Integrands with singular weight functions:: +* Integrands without weight functions:: +* Integrands with weight functions:: +* Integrands with singular weight functions:: File: gsl-ref.info, Node: Integrands without weight functions, Next: Integrands with weight functions, Up: Numerical Integration Introduction @@ -3381,10 +3350,10 @@ File: gsl-ref.info, Node: Integrands without weight functions, Next: Integrand The algorithms for general functions (without a weight function) are based on Gauss-Kronrod rules. - A Gauss-Kronrod rule begins with a classical Gaussian quadrature -rule of order m. This is extended with additional points between each -of the abscissae to give a higher order Kronrod rule of order 2m+1. -The Kronrod rule is efficient because it reuses existing function + A Gauss-Kronrod rule begins with a classical Gaussian quadrature rule +of order m. This is extended with additional points between each of the +abscissae to give a higher order Kronrod rule of order 2m+1. The +Kronrod rule is efficient because it reuses existing function evaluations from the Gaussian rule. The higher order Kronrod rule is used as the best approximation to @@ -3400,12 +3369,11 @@ File: gsl-ref.info, Node: Integrands with weight functions, Next: Integrands w For integrands with weight functions the algorithms use Clenshaw-Curtis quadrature rules. - A Clenshaw-Curtis rule begins with an n-th order Chebyshev -polynomial approximation to the integrand. This polynomial can be -integrated exactly to give an approximation to the integral of the -original function. The Chebyshev expansion can be extended to higher -orders to improve the approximation and provide an estimate of the -error. + A Clenshaw-Curtis rule begins with an n-th order Chebyshev polynomial +approximation to the integrand. This polynomial can be integrated +exactly to give an approximation to the integral of the original +function. The Chebyshev expansion can be extended to higher orders to +improve the approximation and provide an estimate of the error. File: gsl-ref.info, Node: Integrands with singular weight functions, Prev: Integrands with weight functions, Up: Numerical Integration Introduction @@ -3436,9 +3404,10 @@ The QNG algorithm is a non-adaptive procedure which uses fixed Gauss-Kronrod abscissae to sample the integrand at a maximum of 87 points. It is provided for fast integration of smooth functions. - -- Function: int gsl_integration_qng (const gsl_function * F, double - A, double B, double EPSABS, double EPSREL, double * RESULT, + -- Function: int gsl_integration_qng (const gsl_function * F, double A, + double B, double EPSABS, double EPSREL, double * RESULT, double * ABSERR, size_t * NEVAL) + This function applies the Gauss-Kronrod 10-point, 21-point, 43-point and 87-point integration rules in succession until an estimate of the integral of f over (a,b) is achieved within the @@ -3460,12 +3429,12 @@ integration region is divided into subintervals, and on each iteration the subinterval with the largest estimated error is bisected. This reduces the overall error rapidly, as the subintervals become concentrated around local difficulties in the integrand. These -subintervals are managed by a `gsl_integration_workspace' struct, which +subintervals are managed by a 'gsl_integration_workspace' struct, which handles the memory for the subinterval ranges, results and error estimates. -- Function: gsl_integration_workspace * -gsl_integration_workspace_alloc (size_t N) + gsl_integration_workspace_alloc (size_t N) This function allocates a workspace sufficient to hold N double precision intervals, their integration results and error estimates. @@ -3473,17 +3442,18 @@ gsl_integration_workspace_alloc (size_t N) (gsl_integration_workspace * W) This function frees the memory associated with the workspace W. - -- Function: int gsl_integration_qag (const gsl_function * F, double - A, double B, double EPSABS, double EPSREL, size_t LIMIT, int - KEY, gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + -- Function: int gsl_integration_qag (const gsl_function * F, double A, + double B, double EPSABS, double EPSREL, size_t LIMIT, int KEY, + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) + This function applies an integration rule adaptively until an estimate of the integral of f over (a,b) is achieved within the - desired absolute and relative error limits, EPSABS and EPSREL. - The function returns the final approximation, RESULT, and an - estimate of the absolute error, ABSERR. The integration rule is - determined by the value of KEY, which should be chosen from the - following symbolic names, + desired absolute and relative error limits, EPSABS and EPSREL. The + function returns the final approximation, RESULT, and an estimate + of the absolute error, ABSERR. The integration rule is determined + by the value of KEY, which should be chosen from the following + symbolic names, GSL_INTEG_GAUSS15 (key = 1) GSL_INTEG_GAUSS21 (key = 2) @@ -3520,8 +3490,9 @@ integrable singularities. -- Function: int gsl_integration_qags (const gsl_function * F, double A, double B, double EPSABS, double EPSREL, size_t LIMIT, - gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) + This function applies the Gauss-Kronrod 21-point integration rule adaptively until an estimate of the integral of f over (a,b) is achieved within the desired absolute and relative error limits, @@ -3535,24 +3506,24 @@ integrable singularities. is given by LIMIT, which may not exceed the allocated size of the workspace. - File: gsl-ref.info, Node: QAGP adaptive integration with known singular points, Next: QAGI adaptive integration on infinite intervals, Prev: QAGS adaptive integration with singularities, Up: Numerical Integration 16.5 QAGP adaptive integration with known singular points ========================================================= - -- Function: int gsl_integration_qagp (const gsl_function * F, double - * PTS, size_t NPTS, double EPSABS, double EPSREL, size_t - LIMIT, gsl_integration_workspace * WORKSPACE, double * - RESULT, double * ABSERR) + -- Function: int gsl_integration_qagp (const gsl_function * F, double * + PTS, size_t NPTS, double EPSABS, double EPSREL, size_t LIMIT, + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) + This function applies the adaptive integration algorithm QAGS taking account of the user-supplied locations of singular points. The array PTS of length NPTS should contain the endpoints of the integration ranges defined by the integration region and locations of the singularities. For example, to integrate over the region - (a,b) with break-points at x_1, x_2, x_3 (where a < x_1 < x_2 < - x_3 < b) the following PTS array should be used + (a,b) with break-points at x_1, x_2, x_3 (where a < x_1 < x_2 < x_3 + < b) the following PTS array should be used pts[0] = a pts[1] = x_1 @@ -3563,8 +3534,7 @@ File: gsl-ref.info, Node: QAGP adaptive integration with known singular points, with NPTS = 5. If you know the locations of the singular points in the integration - region then this routine will be faster than `QAGS'. - + region then this routine will be faster than 'QAGS'. File: gsl-ref.info, Node: QAGI adaptive integration on infinite intervals, Next: QAWC adaptive integration for Cauchy principal values, Prev: QAGP adaptive integration with known singular points, Up: Numerical Integration @@ -3572,10 +3542,10 @@ File: gsl-ref.info, Node: QAGI adaptive integration on infinite intervals, Nex 16.6 QAGI adaptive integration on infinite intervals ==================================================== - -- Function: int gsl_integration_qagi (gsl_function * F, double - EPSABS, double EPSREL, size_t LIMIT, - gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + -- Function: int gsl_integration_qagi (gsl_function * F, double EPSABS, + double EPSREL, size_t LIMIT, gsl_integration_workspace * + WORKSPACE, double * RESULT, double * ABSERR) + This function computes the integral of the function F over the infinite interval (-\infty,+\infty). The integral is mapped onto the semi-open interval (0,1] using the transformation x = (1-t)/t, @@ -3584,15 +3554,15 @@ File: gsl-ref.info, Node: QAGI adaptive integration on infinite intervals, Nex \int_0^1 dt (f((1-t)/t) + f((-1+t)/t))/t^2. It is then integrated using the QAGS algorithm. The normal - 21-point Gauss-Kronrod rule of QAGS is replaced by a 15-point - rule, because the transformation can generate an integrable - singularity at the origin. In this case a lower-order rule is - more efficient. + 21-point Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, + because the transformation can generate an integrable singularity + at the origin. In this case a lower-order rule is more efficient. -- Function: int gsl_integration_qagiu (gsl_function * F, double A, double EPSABS, double EPSREL, size_t LIMIT, - gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) + This function computes the integral of the function F over the semi-infinite interval (a,+\infty). The integral is mapped onto the semi-open interval (0,1] using the transformation x = a + @@ -3605,8 +3575,8 @@ File: gsl-ref.info, Node: QAGI adaptive integration on infinite intervals, Nex -- Function: int gsl_integration_qagil (gsl_function * F, double B, double EPSABS, double EPSREL, size_t LIMIT, - gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) This function computes the integral of the function F over the semi-infinite interval (-\infty,b). The integral is mapped onto the semi-open interval (0,1] using the transformation x = b - @@ -3625,21 +3595,20 @@ File: gsl-ref.info, Node: QAWC adaptive integration for Cauchy principal values -- Function: int gsl_integration_qawc (gsl_function * F, double A, double B, double C, double EPSABS, double EPSREL, size_t - LIMIT, gsl_integration_workspace * WORKSPACE, double * - RESULT, double * ABSERR) + LIMIT, gsl_integration_workspace * WORKSPACE, double * RESULT, + double * ABSERR) + This function computes the Cauchy principal value of the integral of f over (a,b), with a singularity at C, I = \int_a^b dx f(x) / (x - c) - The adaptive bisection algorithm of QAG is used, with - modifications to ensure that subdivisions do not occur at the - singular point x = c. When a subinterval contains the point x = c - or is close to it then a special 25-point modified Clenshaw-Curtis - rule is used to control the singularity. Further away from the - singularity the algorithm uses an ordinary 15-point Gauss-Kronrod - integration rule. - + The adaptive bisection algorithm of QAG is used, with modifications + to ensure that subdivisions do not occur at the singular point x = + c. When a subinterval contains the point x = c or is close to it + then a special 25-point modified Clenshaw-Curtis rule is used to + control the singularity. Further away from the singularity the + algorithm uses an ordinary 15-point Gauss-Kronrod integration rule. File: gsl-ref.info, Node: QAWS adaptive integration for singular functions, Next: QAWO adaptive integration for oscillatory functions, Prev: QAWC adaptive integration for Cauchy principal values, Up: Numerical Integration @@ -3649,13 +3618,14 @@ File: gsl-ref.info, Node: QAWS adaptive integration for singular functions, Ne The QAWS algorithm is designed for integrands with algebraic-logarithmic singularities at the end-points of an integration region. In order to -work efficiently the algorithm requires a precomputed table of -Chebyshev moments. +work efficiently the algorithm requires a precomputed table of Chebyshev +moments. -- Function: gsl_integration_qaws_table * -gsl_integration_qaws_table_alloc (double ALPHA, double BETA, int MU, - int NU) - This function allocates space for a `gsl_integration_qaws_table' + gsl_integration_qaws_table_alloc (double ALPHA, double BETA, + int MU, int NU) + + This function allocates space for a 'gsl_integration_qaws_table' struct describing a singular weight function W(x) with the parameters (\alpha, \beta, \mu, \nu), @@ -3675,25 +3645,26 @@ gsl_integration_qaws_table_alloc (double ALPHA, double BETA, int MU, integration range. The function returns a pointer to the newly allocated table - `gsl_integration_qaws_table' if no errors were detected, and 0 in + 'gsl_integration_qaws_table' if no errors were detected, and 0 in the case of error. -- Function: int gsl_integration_qaws_table_set (gsl_integration_qaws_table * T, double ALPHA, double BETA, int MU, int NU) This function modifies the parameters (\alpha, \beta, \mu, \nu) of - an existing `gsl_integration_qaws_table' struct T. + an existing 'gsl_integration_qaws_table' struct T. -- Function: void gsl_integration_qaws_table_free (gsl_integration_qaws_table * T) This function frees all the memory associated with the - `gsl_integration_qaws_table' struct T. + 'gsl_integration_qaws_table' struct T. -- Function: int gsl_integration_qaws (gsl_function * F, const double A, const double B, gsl_integration_qaws_table * T, const double EPSABS, const double EPSREL, const size_t LIMIT, - gsl_integration_workspace * WORKSPACE, double * RESULT, - double * ABSERR) + gsl_integration_workspace * WORKSPACE, double * RESULT, double + * ABSERR) + This function computes the integral of the function f(x) over the interval (a,b) with the singular weight function (x-a)^\alpha (b-x)^\beta \log^\mu (x-a) \log^\nu (b-x). The parameters of the @@ -3704,11 +3675,9 @@ gsl_integration_qaws_table_alloc (double ALPHA, double BETA, int MU, The adaptive bisection algorithm of QAG is used. When a subinterval contains one of the endpoints then a special 25-point - modified Clenshaw-Curtis rule is used to control the - singularities. For subintervals which do not include the - endpoints an ordinary 15-point Gauss-Kronrod integration rule is - used. - + modified Clenshaw-Curtis rule is used to control the singularities. + For subintervals which do not include the endpoints an ordinary + 15-point Gauss-Kronrod integration rule is used. File: gsl-ref.info, Node: QAWO adaptive integration for oscillatory functions, Next: QAWF adaptive integration for Fourier integrals, Prev: QAWS adaptive integration for singular functions, Up: Numerical Integration @@ -3722,9 +3691,10 @@ the algorithm requires a table of Chebyshev moments which must be pre-computed with calls to the functions below. -- Function: gsl_integration_qawo_table * -gsl_integration_qawo_table_alloc (double OMEGA, double L, enum + gsl_integration_qawo_table_alloc (double OMEGA, double L, enum gsl_integration_qawo_enum SINE, size_t N) - This function allocates space for a `gsl_integration_qawo_table' + + This function allocates space for a 'gsl_integration_qawo_table' struct and its associated workspace describing a sine or cosine weight function W(x) with the parameters (\omega, L), @@ -3739,16 +3709,15 @@ gsl_integration_qawo_table_alloc (double OMEGA, double L, enum GSL_INTEG_COSINE GSL_INTEG_SINE - The `gsl_integration_qawo_table' is a table of the trigonometric + The 'gsl_integration_qawo_table' is a table of the trigonometric coefficients required in the integration process. The parameter N determines the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that N levels are sufficient for subintervals down to the length L/2^n. - The integration routine `gsl_integration_qawo' returns the error - `GSL_ETABLE' if the number of levels is insufficient for the + The integration routine 'gsl_integration_qawo' returns the error + 'GSL_ETABLE' if the number of levels is insufficient for the requested accuracy. - -- Function: int gsl_integration_qawo_table_set (gsl_integration_qawo_table * T, double OMEGA, double L, enum gsl_integration_qawo_enum SINE) @@ -3769,9 +3738,10 @@ gsl_integration_qawo_table_alloc (double OMEGA, double L, enum LIMIT, gsl_integration_workspace * WORKSPACE, gsl_integration_qawo_table * WF, double * RESULT, double * ABSERR) + This function uses an adaptive algorithm to compute the integral of - f over (a,b) with the weight function \sin(\omega x) or - \cos(\omega x) defined by the table WF, + f over (a,b) with the weight function \sin(\omega x) or \cos(\omega + x) defined by the table WF, I = \int_a^b dx f(x) sin(omega x) I = \int_a^b dx f(x) cos(omega x) @@ -3779,8 +3749,8 @@ gsl_integration_qawo_table_alloc (double OMEGA, double L, enum The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final approximation from the extrapolation, RESULT, and an - estimate of the absolute error, ABSERR. The subintervals and - their results are stored in the memory provided by WORKSPACE. The + estimate of the absolute error, ABSERR. The subintervals and their + results are stored in the memory provided by WORKSPACE. The maximum number of subintervals is given by LIMIT, which may not exceed the allocated size of the workspace. @@ -3790,7 +3760,6 @@ gsl_integration_qawo_table_alloc (double OMEGA, double L, enum widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration. - File: gsl-ref.info, Node: QAWF adaptive integration for Fourier integrals, Next: Numerical integration error codes, Prev: QAWO adaptive integration for oscillatory functions, Up: Numerical Integration @@ -3803,6 +3772,7 @@ File: gsl-ref.info, Node: QAWF adaptive integration for Fourier integrals, Nex gsl_integration_workspace * CYCLE_WORKSPACE, gsl_integration_qawo_table * WF, double * RESULT, double * ABSERR) + This function attempts to compute a Fourier integral of the function F over the semi-infinite interval [a,+\infty). @@ -3845,11 +3815,10 @@ File: gsl-ref.info, Node: QAWF adaptive integration for Fourier integrals, Nex where E_k is the estimated error on the interval C_k. The subintervals and their results are stored in the memory - provided by WORKSPACE. The maximum number of subintervals is - given by LIMIT, which may not exceed the allocated size of the - workspace. The integration over each subinterval uses the memory - provided by CYCLE_WORKSPACE as workspace for the QAWO algorithm. - + provided by WORKSPACE. The maximum number of subintervals is given + by LIMIT, which may not exceed the allocated size of the workspace. + The integration over each subinterval uses the memory provided by + CYCLE_WORKSPACE as workspace for the QAWO algorithm. File: gsl-ref.info, Node: Numerical integration error codes, Next: Numerical integration examples, Prev: QAWF adaptive integration for Fourier integrals, Up: Numerical Integration @@ -3860,18 +3829,15 @@ File: gsl-ref.info, Node: Numerical integration error codes, Next: Numerical i In addition to the standard error codes for invalid arguments the functions can return the following values, -`GSL_EMAXITER' +'GSL_EMAXITER' the maximum number of subdivisions was exceeded. - -`GSL_EROUND' - cannot reach tolerance because of roundoff error, or roundoff - error was detected in the extrapolation table. - -`GSL_ESING' +'GSL_EROUND' + cannot reach tolerance because of roundoff error, or roundoff error + was detected in the extrapolation table. +'GSL_ESING' a non-integrable singularity or other bad integrand behavior was found in the integration interval. - -`GSL_EDIVERGE' +'GSL_EDIVERGE' the integral is divergent, or too slowly convergent to be integrated numerically. @@ -3881,14 +3847,14 @@ File: gsl-ref.info, Node: Numerical integration examples, Next: Numerical inte 16.12 Examples ============== -The integrator `QAGS' will handle a large class of definite integrals. +The integrator 'QAGS' will handle a large class of definite integrals. For example, consider the following integral, which has a algebraic-logarithmic singularity at the origin, \int_0^1 x^{-1/2} log(x) dx = -4 The program below computes this integral to a relative accuracy bound of -`1e-7'. +'1e-7'. #include <stdio.h> #include <math.h> @@ -3938,10 +3904,10 @@ subdivisions. actual error = 0.000000000000026201 intervals = 8 -In fact, the extrapolation procedure used by `QAGS' produces an -accuracy of almost twice as many digits. The error estimate returned by -the extrapolation procedure is larger than the actual error, giving a -margin of safety of one order of magnitude. +In fact, the extrapolation procedure used by 'QAGS' produces an accuracy +of almost twice as many digits. The error estimate returned by the +extrapolation procedure is larger than the actual error, giving a margin +of safety of one order of magnitude. File: gsl-ref.info, Node: Numerical integration References and Further Reading, Prev: Numerical integration examples, Up: Numerical Integration @@ -3956,10 +3922,9 @@ includes useful advice on numerical integration and many references to the numerical integration literature used in developing QUADPACK. R. Piessens, E. de Doncker-Kapenga, C.W. Uberhuber, D.K. Kahaner. - `QUADPACK A subroutine package for automatic integration' Springer + 'QUADPACK A subroutine package for automatic integration' Springer Verlag, 1983. - File: gsl-ref.info, Node: Random Number Generation, Next: Quasi-Random Sequences, Prev: Numerical Integration, Up: Top @@ -3976,25 +3941,25 @@ programs. Additional functions are available for transforming uniform random numbers into samples from continuous or discrete probability distributions such as the Gaussian, log-normal or Poisson distributions. - These functions are declared in the header file `gsl_rng.h'. + These functions are declared in the header file 'gsl_rng.h'. * Menu: -* General comments on random numbers:: -* The Random Number Generator Interface:: -* Random number generator initialization:: -* Sampling from a random number generator:: -* Auxiliary random number generator functions:: -* Random number environment variables:: -* Copying random number generator state:: -* Reading and writing random number generator state:: -* Random number generator algorithms:: -* Unix random number generators:: -* Other random number generators:: -* Random Number Generator Performance:: -* Random Number Generator Examples:: -* Random Number References and Further Reading:: -* Random Number Acknowledgements:: +* General comments on random numbers:: +* The Random Number Generator Interface:: +* Random number generator initialization:: +* Sampling from a random number generator:: +* Auxiliary random number generator functions:: +* Random number environment variables:: +* Copying random number generator state:: +* Reading and writing random number generator state:: +* Random number generator algorithms:: +* Unix random number generators:: +* Other random number generators:: +* Random Number Generator Performance:: +* Random Number Generator Examples:: +* Random Number References and Further Reading:: +* Random Number Acknowledgements:: File: gsl-ref.info, Node: General comments on random numbers, Next: The Random Number Generator Interface, Up: Random Number Generation @@ -4003,7 +3968,7 @@ File: gsl-ref.info, Node: General comments on random numbers, Next: The Random ======================================= In 1988, Park and Miller wrote a paper entitled "Random number -generators: good ones are hard to find." [Commun. ACM, 31, 1192-1201]. +generators: good ones are hard to find." [Commun. ACM, 31, 1192-1201]. Fortunately, some excellent random number generators are available, though poor ones are still in common use. You may be happy with the system-supplied random number generator on your computer, but you should @@ -4014,7 +3979,7 @@ the hall to the coffee machine and back. A very nice review of random number generators was written by Pierre L'Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks, -ed. (Wiley, 1997). The chapter is available in postscript from +ed. (Wiley, 1997). The chapter is available in postscript from L'Ecuyer's ftp site (see references). Knuth's volume on Seminumerical Algorithms (originally published in 1968) devotes 170 pages to random number generators, and has recently been updated in its 3rd edition @@ -4050,12 +4015,12 @@ generators, you do not need to know the details of what comprises the state, and besides that varies from algorithm to algorithm. The random number generator library uses two special structs, -`gsl_rng_type' which holds static information about each type of -generator and `gsl_rng' which describes an instance of a generator -created from a given `gsl_rng_type'. +'gsl_rng_type' which holds static information about each type of +generator and 'gsl_rng' which describes an instance of a generator +created from a given 'gsl_rng_type'. The functions described in this section are declared in the header -file `gsl_rng.h'. +file 'gsl_rng.h'. File: gsl-ref.info, Node: Random number generator initialization, Next: Sampling from a random number generator, Prev: The Random Number Generator Interface, Up: Random Number Generation @@ -4065,34 +4030,34 @@ File: gsl-ref.info, Node: Random number generator initialization, Next: Sampli -- Function: gsl_rng * gsl_rng_alloc (const gsl_rng_type * T) This function returns a pointer to a newly-created instance of a - random number generator of type T. For example, the following - code creates an instance of the Tausworthe generator, + random number generator of type T. For example, the following code + creates an instance of the Tausworthe generator, gsl_rng * r = gsl_rng_alloc (gsl_rng_taus); If there is insufficient memory to create the generator then the function returns a null pointer and the error handler is invoked - with an error code of `GSL_ENOMEM'. + with an error code of 'GSL_ENOMEM'. The generator is automatically initialized with the default seed, - `gsl_rng_default_seed'. This is zero by default but can be changed - either directly or by using the environment variable `GSL_RNG_SEED' + 'gsl_rng_default_seed'. This is zero by default but can be changed + either directly or by using the environment variable 'GSL_RNG_SEED' (*note Random number environment variables::). - The details of the available generator types are described later - in this chapter. + The details of the available generator types are described later in + this chapter. -- Function: void gsl_rng_set (const gsl_rng * R, unsigned long int S) - This function initializes (or `seeds') the random number - generator. If the generator is seeded with the same value of S on - two different runs, the same stream of random numbers will be - generated by successive calls to the routines below. If different - values of S are supplied, then the generated streams of random - numbers should be completely different. If the seed S is zero - then the standard seed from the original implementation is used - instead. For example, the original Fortran source code for the - `ranlux' generator used a seed of 314159265, and so choosing S - equal to zero reproduces this when using `gsl_rng_ranlux'. + This function initializes (or 'seeds') the random number generator. + If the generator is seeded with the same value of S on two + different runs, the same stream of random numbers will be generated + by successive calls to the routines below. If different values of + S are supplied, then the generated streams of random numbers should + be completely different. If the seed S is zero then the standard + seed from the original implementation is used instead. For + example, the original Fortran source code for the 'ranlux' + generator used a seed of 314159265, and so choosing S equal to zero + reproduces this when using 'gsl_rng_ranlux'. -- Function: void gsl_rng_free (gsl_rng * R) This function frees all the memory associated with the generator R. @@ -4112,25 +4077,25 @@ obtain non-uniform distributions *note Random Number Distributions::. minimum and maximum values depend on the algorithm used, but all integers in the range [MIN,MAX] are equally likely. The values of MIN and MAX can determined using the auxiliary functions - `gsl_rng_max (r)' and `gsl_rng_min (r)'. + 'gsl_rng_max (r)' and 'gsl_rng_min (r)'. -- Function: double gsl_rng_uniform (const gsl_rng * R) This function returns a double precision floating point number uniformly distributed in the range [0,1). The range includes 0.0 but excludes 1.0. The value is typically obtained by dividing the - result of `gsl_rng_get(r)' by `gsl_rng_max(r) + 1.0' in double + result of 'gsl_rng_get(r)' by 'gsl_rng_max(r) + 1.0' in double precision. Some generators compute this ratio internally so that they can provide floating point numbers with more than 32 bits of randomness (the maximum number of bits that can be portably - represented in a single `unsigned long int'). + represented in a single 'unsigned long int'). -- Function: double gsl_rng_uniform_pos (const gsl_rng * R) This function returns a positive double precision floating point - number uniformly distributed in the range (0,1), excluding both - 0.0 and 1.0. The number is obtained by sampling the generator - with the algorithm of `gsl_rng_uniform' until a non-zero value is - obtained. You can use this function if you need to avoid a - singularity at 0.0. + number uniformly distributed in the range (0,1), excluding both 0.0 + and 1.0. The number is obtained by sampling the generator with the + algorithm of 'gsl_rng_uniform' until a non-zero value is obtained. + You can use this function if you need to avoid a singularity at + 0.0. -- Function: unsigned long int gsl_rng_uniform_int (const gsl_rng * R, unsigned long int N) @@ -4142,16 +4107,16 @@ obtain non-uniform distributions *note Random Number Distributions::. Note that this function is designed for sampling from ranges smaller than the range of the underlying generator. The parameter - N must be less than or equal to the range of the generator R. If - N is larger than the range of the generator then the function - calls the error handler with an error code of `GSL_EINVAL' and - returns zero. + N must be less than or equal to the range of the generator R. If N + is larger than the range of the generator then the function calls + the error handler with an error code of 'GSL_EINVAL' and returns + zero. In particular, this function is not intended for generating the - full range of unsigned integer values [0,2^32-1]. Instead choose a + full range of unsigned integer values [0,2^32-1]. Instead choose a generator with the maximal integer range and zero mimimum value, - such as `gsl_rng_ranlxd1', `gsl_rng_mt19937' or `gsl_rng_taus', - and sample it directly using `gsl_rng_get'. The range of each + such as 'gsl_rng_ranlxd1', 'gsl_rng_mt19937' or 'gsl_rng_taus', and + sample it directly using 'gsl_rng_get'. The range of each generator can be found using the auxiliary functions described in the next section. @@ -4161,9 +4126,9 @@ File: gsl-ref.info, Node: Auxiliary random number generator functions, Next: R 17.5 Auxiliary random number generator functions ================================================ -The following functions provide information about an existing -generator. You should use them in preference to hard-coding the -generator parameters into your own code. +The following functions provide information about an existing generator. +You should use them in preference to hard-coding the generator +parameters into your own code. -- Function: const char * gsl_rng_name (const gsl_rng * R) This function returns a pointer to the name of the generator. For @@ -4172,14 +4137,14 @@ generator parameters into your own code. printf ("r is a '%s' generator\n", gsl_rng_name (r)); - would print something like `r is a 'taus' generator'. + would print something like 'r is a 'taus' generator'. -- Function: unsigned long int gsl_rng_max (const gsl_rng * R) - `gsl_rng_max' returns the largest value that `gsl_rng_get' can + 'gsl_rng_max' returns the largest value that 'gsl_rng_get' can return. -- Function: unsigned long int gsl_rng_min (const gsl_rng * R) - `gsl_rng_min' returns the smallest value that `gsl_rng_get' can + 'gsl_rng_min' returns the smallest value that 'gsl_rng_get' can return. Usually this value is zero. There are some generators with algorithms that cannot return zero, and for these generators the minimum value is 1. @@ -4197,7 +4162,7 @@ generator parameters into your own code. -- Function: const gsl_rng_type ** gsl_rng_types_setup (void) This function returns a pointer to an array of all the available - generator types, terminated by a null pointer. The function should + generator types, terminated by a null pointer. The function should be called once at the start of the program, if needed. The following code fragment shows how to iterate over the array of generator types to print the names of the available algorithms, @@ -4220,32 +4185,31 @@ File: gsl-ref.info, Node: Random number environment variables, Next: Copying r ======================================== The library allows you to choose a default generator and seed from the -environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED' and the -function `gsl_rng_env_setup'. This makes it easy try out different -generators and seeds without having to recompile your program. +environment variables 'GSL_RNG_TYPE' and 'GSL_RNG_SEED' and the function +'gsl_rng_env_setup'. This makes it easy try out different generators +and seeds without having to recompile your program. -- Function: const gsl_rng_type * gsl_rng_env_setup (void) - This function reads the environment variables `GSL_RNG_TYPE' and - `GSL_RNG_SEED' and uses their values to set the corresponding - library variables `gsl_rng_default' and `gsl_rng_default_seed'. + This function reads the environment variables 'GSL_RNG_TYPE' and + 'GSL_RNG_SEED' and uses their values to set the corresponding + library variables 'gsl_rng_default' and 'gsl_rng_default_seed'. These global variables are defined as follows, extern const gsl_rng_type *gsl_rng_default extern unsigned long int gsl_rng_default_seed - The environment variable `GSL_RNG_TYPE' should be the name of a - generator, such as `taus' or `mt19937'. The environment variable - `GSL_RNG_SEED' should contain the desired seed value. It is - converted to an `unsigned long int' using the C library function - `strtoul'. - - If you don't specify a generator for `GSL_RNG_TYPE' then - `gsl_rng_mt19937' is used as the default. The initial value of - `gsl_rng_default_seed' is zero. + The environment variable 'GSL_RNG_TYPE' should be the name of a + generator, such as 'taus' or 'mt19937'. The environment variable + 'GSL_RNG_SEED' should contain the desired seed value. It is + converted to an 'unsigned long int' using the C library function + 'strtoul'. + If you don't specify a generator for 'GSL_RNG_TYPE' then + 'gsl_rng_mt19937' is used as the default. The initial value of + 'gsl_rng_default_seed' is zero. Here is a short program which shows how to create a global generator -using the environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED', +using the environment variables 'GSL_RNG_TYPE' and 'GSL_RNG_SEED', #include <stdio.h> #include <gsl/gsl_rng.h> @@ -4271,7 +4235,7 @@ using the environment variables `GSL_RNG_TYPE' and `GSL_RNG_SEED', } Running the program without any environment variables uses the initial -defaults, an `mt19937' generator with a seed of 0, +defaults, an 'mt19937' generator with a seed of 0, $ ./a.out generator type: mt19937 @@ -4294,15 +4258,15 @@ File: gsl-ref.info, Node: Copying random number generator state, Next: Reading 17.7 Copying random number generator state ========================================== -The above methods do not expose the random number `state' which changes +The above methods do not expose the random number 'state' which changes from call to call. It is often useful to be able to save and restore the state. To permit these practices, a few somewhat more advanced functions are supplied. These include: -- Function: int gsl_rng_memcpy (gsl_rng * DEST, const gsl_rng * SRC) This function copies the random number generator SRC into the - pre-existing generator DEST, making DEST into an exact copy of - SRC. The two generators must be of the same type. + pre-existing generator DEST, making DEST into an exact copy of SRC. + The two generators must be of the same type. -- Function: gsl_rng * gsl_rng_clone (const gsl_rng * R) This function returns a pointer to a newly created generator which @@ -4314,13 +4278,13 @@ File: gsl-ref.info, Node: Reading and writing random number generator state, N 17.8 Reading and writing random number generator state ====================================================== -The library provides functions for reading and writing the random -number state to a file as binary data or formatted text. +The library provides functions for reading and writing the random number +state to a file as binary data or formatted text. -- Function: int gsl_rng_fwrite (FILE * STREAM, const gsl_rng * R) This function writes the random number state of the random number generator R to the stream STREAM in binary format. The return - value is 0 for success and `GSL_EFAILED' if there was a problem + value is 0 for success and 'GSL_EFAILED' if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures. @@ -4330,7 +4294,7 @@ number state to a file as binary data or formatted text. generator R from the open stream STREAM in binary format. The random number generator R must be preinitialized with the correct random number generator type since type information is not saved. - The return value is 0 for success and `GSL_EFAILED' if there was a + The return value is 0 for success and 'GSL_EFAILED' if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture. @@ -4362,70 +4326,72 @@ randomness. statistical tests. It uses 624 words of state per generator and is comparable in speed to the other generators. The original generator used a default seed of 4357 and choosing S equal to zero - in `gsl_rng_set' reproduces this. Later versions switched to 5489 + in 'gsl_rng_set' reproduces this. Later versions switched to 5489 as the default seed, you can choose this explicitly via - `gsl_rng_set' instead if you require it. + 'gsl_rng_set' instead if you require it. For more information see, Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number - generator". `ACM Transactions on Modeling and Computer - Simulation', Vol. 8, No. 1 (Jan. 1998), Pages 3-30 + generator". 'ACM Transactions on Modeling and Computer + Simulation', Vol. 8, No. 1 (Jan. 1998), Pages 3-30 - The generator `gsl_rng_mt19937' uses the second revision of the + The generator 'gsl_rng_mt19937' uses the second revision of the seeding procedure published by the two authors above in 2002. The - original seeding procedures could cause spurious artifacts for - some seed values. They are still available through the alternative - generators `gsl_rng_mt19937_1999' and `gsl_rng_mt19937_1998'. + original seeding procedures could cause spurious artifacts for some + seed values. They are still available through the alternative + generators 'gsl_rng_mt19937_1999' and 'gsl_rng_mt19937_1998'. -- Generator: gsl_rng_ranlxs0 -- Generator: gsl_rng_ranlxs1 -- Generator: gsl_rng_ranlxs2 - The generator `ranlxs0' is a second-generation version of the + + The generator 'ranlxs0' is a second-generation version of the RANLUX algorithm of Lu"scher, which produces "luxury random numbers". This generator provides single precision output (24 - bits) at three luxury levels `ranlxs0', `ranlxs1' and `ranlxs2', - in increasing order of strength. It uses double-precision - floating point arithmetic internally and can be significantly - faster than the integer version of `ranlux', particularly on - 64-bit architectures. The period of the generator is about - 10^171. The algorithm has mathematically proven properties and - can provide truly decorrelated numbers at a known level of - randomness. The higher luxury levels provide increased - decorrelation between samples as an additional safety margin. + bits) at three luxury levels 'ranlxs0', 'ranlxs1' and 'ranlxs2', in + increasing order of strength. It uses double-precision floating + point arithmetic internally and can be significantly faster than + the integer version of 'ranlux', particularly on 64-bit + architectures. The period of the generator is about 10^171. The + algorithm has mathematically proven properties and can provide + truly decorrelated numbers at a known level of randomness. The + higher luxury levels provide increased decorrelation between + samples as an additional safety margin. -- Generator: gsl_rng_ranlxd1 -- Generator: gsl_rng_ranlxd2 + These generators produce double precision output (48 bits) from the - RANLXS generator. The library provides two luxury levels - `ranlxd1' and `ranlxd2', in increasing order of strength. + RANLXS generator. The library provides two luxury levels 'ranlxd1' + and 'ranlxd2', in increasing order of strength. -- Generator: gsl_rng_ranlux -- Generator: gsl_rng_ranlux389 - The `ranlux' generator is an implementation of the original + + The 'ranlux' generator is an implementation of the original algorithm developed by Lu"scher. It uses a lagged-fibonacci-with-skipping algorithm to produce "luxury random numbers". It is a 24-bit generator, originally designed for single-precision IEEE floating point numbers. This implementation is based on integer arithmetic, while the second-generation versions RANLXS and RANLXD described above provide floating-point - implementations which will be faster on many platforms. The - period of the generator is about 10^171. The algorithm has - mathematically proven properties and it can provide truly - decorrelated numbers at a known level of randomness. The default - level of decorrelation recommended by Lu"scher is provided by - `gsl_rng_ranlux', while `gsl_rng_ranlux389' gives the highest - level of randomness, with all 24 bits decorrelated. Both types of - generator use 24 words of state per generator. + implementations which will be faster on many platforms. The period + of the generator is about 10^171. The algorithm has mathematically + proven properties and it can provide truly decorrelated numbers at + a known level of randomness. The default level of decorrelation + recommended by Lu"scher is provided by 'gsl_rng_ranlux', while + 'gsl_rng_ranlux389' gives the highest level of randomness, with all + 24 bits decorrelated. Both types of generator use 24 words of + state per generator. For more information see, M. Lu"scher, "A portable high-quality random number generator - for lattice field theory calculations", `Computer Physics + for lattice field theory calculations", 'Computer Physics Communications', 79 (1994) 100-110. - F. James, "RANLUX: A Fortran implementation of the high-quality pseudo-random number generator of Lu"scher", - `Computer Physics Communications', 79 (1994) 111-114 + 'Computer Physics Communications', 79 (1994) 111-114 -- Generator: gsl_rng_cmrg This is a combined multiple recursive generator by L'Ecuyer. Its @@ -4447,7 +4413,7 @@ randomness. generator. For more information see, P. L'Ecuyer, "Combined Multiple Recursive Random Number - Generators", `Operations Research', 44, 5 (1996), 816-822. + Generators", 'Operations Research', 44, 5 (1996), 816-822. -- Generator: gsl_rng_mrg This is a fifth-order multiple recursive generator by L'Ecuyer, @@ -4462,7 +4428,7 @@ randomness. state per generator. More information can be found in the following paper, P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good - multiple recursive random number generators", `ACM + multiple recursive random number generators", 'ACM Transactions on Modeling and Computer Simulation' 3, 87-98 (1993). @@ -4482,29 +4448,29 @@ randomness. computed modulo 2^32. In the formulas above ^^ denotes "exclusive-or". Note that the algorithm relies on the properties of 32-bit unsigned integers and has been implemented using a - bitmask of `0xFFFFFFFF' to make it work on 64 bit machines. + bitmask of '0xFFFFFFFF' to make it work on 64 bit machines. The period of this generator is 2^88 (about 10^26). It uses 3 words of state per generator. For more information see, P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe - Generators", `Mathematics of Computation', 65, 213 (1996), + Generators", 'Mathematics of Computation', 65, 213 (1996), 203-213. - The generator `gsl_rng_taus2' uses the same algorithm as - `gsl_rng_taus' but with an improved seeding procedure described in + The generator 'gsl_rng_taus2' uses the same algorithm as + 'gsl_rng_taus' but with an improved seeding procedure described in the paper, P. L'Ecuyer, "Tables of Maximally Equidistributed Combined - LFSR Generators", `Mathematics of Computation', 68, 225 + LFSR Generators", 'Mathematics of Computation', 68, 225 (1999), 261-269 - The generator `gsl_rng_taus2' should now be used in preference to - `gsl_rng_taus'. + The generator 'gsl_rng_taus2' should now be used in preference to + 'gsl_rng_taus'. -- Generator: gsl_rng_gfsr4 - The `gfsr4' generator is like a lagged-fibonacci generator, and - produces each number as an `xor''d sum of four previous values. + The 'gfsr4' generator is like a lagged-fibonacci generator, and + produces each number as an 'xor''d sum of four previous values. r_n = r_{n-A} ^^ r_{n-B} ^^ r_{n-C} ^^ r_{n-D} @@ -4520,11 +4486,11 @@ randomness. of Ziff's article: A=471, B=1586, C=6988, D=9689. If the offsets are appropriately chosen (such as the one ones in - this implementation), then the sequence is said to be maximal; - that means that the period is 2^D - 1, where D is the longest lag. - (It is one less than 2^D because it is not permitted to have all - zeros in the `ra[]' array.) For this implementation with D=9689 - that works out to about 10^2917. + this implementation), then the sequence is said to be maximal; that + means that the period is 2^D - 1, where D is the longest lag. (It + is one less than 2^D because it is not permitted to have all zeros + in the 'ra[]' array.) For this implementation with D=9689 that + works out to about 10^2917. Note that the implementation of this generator using a 32-bit integer amounts to 32 parallel implementations of one-bit @@ -4537,7 +4503,7 @@ randomness. For more information see, Robert M. Ziff, "Four-tap shift-register-sequence - random-number generators", `Computers in Physics', 12(4), + random-number generators", 'Computers in Physics', 12(4), Jul/Aug 1998, pp 385-392. @@ -4546,19 +4512,19 @@ File: gsl-ref.info, Node: Unix random number generators, Next: Other random nu 17.10 Unix random number generators =================================== -The standard Unix random number generators `rand', `random' and -`rand48' are provided as part of GSL. Although these generators are -widely available individually often they aren't all available on the -same platform. This makes it difficult to write portable code using -them and so we have included the complete set of Unix generators in GSL -for convenience. Note that these generators don't produce high-quality +The standard Unix random number generators 'rand', 'random' and 'rand48' +are provided as part of GSL. Although these generators are widely +available individually often they aren't all available on the same +platform. This makes it difficult to write portable code using them and +so we have included the complete set of Unix generators in GSL for +convenience. Note that these generators don't produce high-quality randomness and aren't suitable for work requiring accurate statistics. However, if you won't be measuring statistical quantities and just want to introduce some variation into your program then these generators are quite acceptable. -- Generator: gsl_rng_rand - This is the BSD `rand' generator. Its sequence is + This is the BSD 'rand' generator. Its sequence is x_{n+1} = (a x_n + c) mod m @@ -4569,20 +4535,20 @@ quite acceptable. -- Generator: gsl_rng_random_bsd -- Generator: gsl_rng_random_libc5 -- Generator: gsl_rng_random_glibc2 - These generators implement the `random' family of functions, a set + These generators implement the 'random' family of functions, a set of linear feedback shift register generators originally used in BSD - Unix. There are several versions of `random' in use today: the - original BSD version (e.g. on SunOS4), a libc5 version (found on + Unix. There are several versions of 'random' in use today: the + original BSD version (e.g. on SunOS4), a libc5 version (found on older GNU/Linux systems) and a glibc2 version. Each version uses a different seeding procedure, and thus produces different sequences. The original BSD routines accepted a variable length buffer for the generator state, with longer buffers providing higher-quality - randomness. The `random' function implemented algorithms for + randomness. The 'random' function implemented algorithms for buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm with the largest length that would fit into the user-supplied - buffer was used. To support these algorithms additional - generators are available with the following names, + buffer was used. To support these algorithms additional generators + are available with the following names, gsl_rng_random8_bsd gsl_rng_random32_bsd @@ -4590,31 +4556,30 @@ quite acceptable. gsl_rng_random128_bsd gsl_rng_random256_bsd - where the numeric suffix indicates the buffer length. The - original BSD `random' function used a 128-byte default buffer and - so `gsl_rng_random_bsd' has been made equivalent to - `gsl_rng_random128_bsd'. Corresponding versions of the `libc5' - and `glibc2' generators are also available, with the names - `gsl_rng_random8_libc5', `gsl_rng_random8_glibc2', etc. + where the numeric suffix indicates the buffer length. The original + BSD 'random' function used a 128-byte default buffer and so + 'gsl_rng_random_bsd' has been made equivalent to + 'gsl_rng_random128_bsd'. Corresponding versions of the 'libc5' and + 'glibc2' generators are also available, with the names + 'gsl_rng_random8_libc5', 'gsl_rng_random8_glibc2', etc. -- Generator: gsl_rng_rand48 - This is the Unix `rand48' generator. Its sequence is + This is the Unix 'rand48' generator. Its sequence is x_{n+1} = (a x_n + c) mod m defined on 48-bit unsigned integers with a = 25214903917, c = 11 and m = 2^48. The seed specifies the upper 32 bits of the initial - value, x_1, with the lower 16 bits set to `0x330E'. The function - `gsl_rng_get' returns the upper 32 bits from each term of the + value, x_1, with the lower 16 bits set to '0x330E'. The function + 'gsl_rng_get' returns the upper 32 bits from each term of the sequence. This does not have a direct parallel in the original - `rand48' functions, but forcing the result to type `long int' - reproduces the output of `mrand48'. The function - `gsl_rng_uniform' uses the full 48 bits of internal state to return - the double precision number x_n/m, which is equivalent to the - function `drand48'. Note that some versions of the GNU C Library - contained a bug in `mrand48' function which caused it to produce - different results (only the lower 16-bits of the return value were - set). + 'rand48' functions, but forcing the result to type 'long int' + reproduces the output of 'mrand48'. The function 'gsl_rng_uniform' + uses the full 48 bits of internal state to return the double + precision number x_n/m, which is equivalent to the function + 'drand48'. Note that some versions of the GNU C Library contained + a bug in 'mrand48' function which caused it to produce different + results (only the lower 16-bits of the return value were set). File: gsl-ref.info, Node: Other random number generators, Next: Random Number Generator Performance, Prev: Unix random number generators, Up: Random Number Generation @@ -4639,14 +4604,14 @@ having any randomness. Thus if you want to produce a random bitstream it is best to avoid using the least significant bits. -- Generator: gsl_rng_ranf - This is the CRAY random number generator `RANF'. Its sequence is + This is the CRAY random number generator 'RANF'. Its sequence is x_{n+1} = (a x_n) mod m - defined on 48-bit unsigned integers with a = 44485709377909 and m - = 2^48. The seed specifies the lower 32 bits of the initial value, + defined on 48-bit unsigned integers with a = 44485709377909 and m = + 2^48. The seed specifies the lower 32 bits of the initial value, x_1, with the lowest bit set to prevent the seed taking an even - value. The upper 16 bits of x_1 are set to 0. A consequence of + value. The upper 16 bits of x_1 are set to 0. A consequence of this procedure is that the pairs of seeds 2 and 3, 4 and 5, etc produce the same sequences. @@ -4659,12 +4624,12 @@ it is best to avoid using the least significant bits. modular inverse of a mod m. This is done for compatibility with the original CRAY implementation. - Note that you can only seed the generator with integers up to - 2^32, while the original CRAY implementation uses non-portable - wide integers which can cover all 2^48 states of the generator. + Note that you can only seed the generator with integers up to 2^32, + while the original CRAY implementation uses non-portable wide + integers which can cover all 2^48 states of the generator. - The function `gsl_rng_get' returns the upper 32 bits from each term - of the sequence. The function `gsl_rng_uniform' uses the full 48 + The function 'gsl_rng_get' returns the upper 32 bits from each term + of the sequence. The function 'gsl_rng_uniform' uses the full 48 bits to return the double precision number x_n/m. The period of this generator is 2^46. @@ -4687,7 +4652,7 @@ it is best to avoid using the least significant bits. For more information see, S. Kirkpatrick and E. Stoll, "A very fast shift-register - sequence random number generator", `Journal of Computational + sequence random number generator", 'Journal of Computational Physics', 40, 517-526 (1981) -- Generator: gsl_rng_tt800 @@ -4699,17 +4664,17 @@ it is best to avoid using the least significant bits. For more information see, Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR - Generators II", `ACM Transactions on Modelling and Computer + Generators II", 'ACM Transactions on Modelling and Computer Simulation', Vol. 4, No. 3, 1994, pages 254-266. -- Generator: gsl_rng_vax - This is the VAX generator `MTH$RANDOM'. Its sequence is, + This is the VAX generator 'MTH$RANDOM'. Its sequence is, x_{n+1} = (a x_n + c) mod m - with a = 69069, c = 1 and m = 2^32. The seed specifies the - initial value, x_1. The period of this generator is 2^32 and it - uses 1 word of storage per generator. + with a = 69069, c = 1 and m = 2^32. The seed specifies the initial + value, x_1. The period of this generator is 2^32 and it uses 1 + word of storage per generator. -- Generator: gsl_rng_transputer This is the random number generator from the INMOS Transputer @@ -4721,7 +4686,7 @@ it is best to avoid using the least significant bits. value, x_1. -- Generator: gsl_rng_randu - This is the IBM `RANDU' generator. Its sequence is + This is the IBM 'RANDU' generator. Its sequence is x_{n+1} = (a x_n) mod m @@ -4745,16 +4710,16 @@ it is best to avoid using the least significant bits. acronym "GGL" (I'm not sure what that stands for). For more information see, - Park and Miller, "Random Number Generators: Good ones are - hard to find", `Communications of the ACM', October 1988, - Volume 31, No 10, pages 1192-1201. + Park and Miller, "Random Number Generators: Good ones are hard + to find", 'Communications of the ACM', October 1988, Volume + 31, No 10, pages 1192-1201. -- Generator: gsl_rng_uni -- Generator: gsl_rng_uni32 This is a reimplementation of the 16-bit SLATEC random number generator RUNIF. A generalization of the generator to 32 bits is - provided by `gsl_rng_uni32'. The original source code is - available from NETLIB. + provided by 'gsl_rng_uni32'. The original source code is available + from NETLIB. -- Generator: gsl_rng_slatec This is the SLATEC random number generator RAND. It is ancient. @@ -4767,15 +4732,15 @@ it is best to avoid using the least significant bits. t = u_{n-273} + u_{n-607} u_n = t - floor(t) - The original source code is available from NETLIB. For more + The original source code is available from NETLIB. For more information see, W. Petersen, "Lagged Fibonacci Random Number Generators for - the NEC SX-3", `International Journal of High Speed - Computing' (1994). + the NEC SX-3", 'International Journal of High Speed Computing' + (1994). -- Generator: gsl_rng_knuthran2 This is a second-order multiple recursive generator described by - Knuth in `Seminumerical Algorithms', 3rd Ed., page 108. Its + Knuth in 'Seminumerical Algorithms', 3rd Ed., page 108. Its sequence is, x_n = (a_1 x_{n-1} + a_2 x_{n-2}) mod m @@ -4785,10 +4750,10 @@ it is best to avoid using the least significant bits. -- Generator: gsl_rng_knuthran2002 -- Generator: gsl_rng_knuthran This is a second-order multiple recursive generator described by - Knuth in `Seminumerical Algorithms', 3rd Ed., Section 3.6. Knuth - provides its C code. The updated routine `gsl_rng_knuthran2002' - is from the revised 9th printing and corrects some weaknesses in - the earlier version, which is implemented as `gsl_rng_knuthran'. + Knuth in 'Seminumerical Algorithms', 3rd Ed., Section 3.6. Knuth + provides its C code. The updated routine 'gsl_rng_knuthran2002' is + from the revised 9th printing and corrects some weaknesses in the + earlier version, which is implemented as 'gsl_rng_knuthran'. -- Generator: gsl_rng_borosh13 -- Generator: gsl_rng_fishman18 @@ -4796,7 +4761,7 @@ it is best to avoid using the least significant bits. -- Generator: gsl_rng_lecuyer21 -- Generator: gsl_rng_waterman14 These multiplicative generators are taken from Knuth's - `Seminumerical Algorithms', 3rd Ed., pages 106-108. Their sequence + 'Seminumerical Algorithms', 3rd Ed., pages 106-108. Their sequence is, x_{n+1} = (a x_n) mod m @@ -4808,19 +4773,18 @@ it is best to avoid using the least significant bits. 1566083941, m = 2^32. -- Generator: gsl_rng_fishman2x - This is the L'Ecuyer-Fishman random number generator. It is taken - from Knuth's `Seminumerical Algorithms', 3rd Ed., page 108. Its + This is the L'Ecuyer-Fishman random number generator. It is taken + from Knuth's 'Seminumerical Algorithms', 3rd Ed., page 108. Its sequence is, z_{n+1} = (x_n - y_n) mod m - with m = 2^31 - 1. x_n and y_n are given by the `fishman20' and - `lecuyer21' algorithms. The seed specifies the initial value, x_1. - + with m = 2^31 - 1. x_n and y_n are given by the 'fishman20' and + 'lecuyer21' algorithms. The seed specifies the initial value, x_1. -- Generator: gsl_rng_coveyou - This is the Coveyou random number generator. It is taken from - Knuth's `Seminumerical Algorithms', 3rd Ed., Section 3.2.2. Its + This is the Coveyou random number generator. It is taken from + Knuth's 'Seminumerical Algorithms', 3rd Ed., Section 3.2.2. Its sequence is, x_{n+1} = (x_n (x_n + 1)) mod m @@ -4835,7 +4799,7 @@ File: gsl-ref.info, Node: Random Number Generator Performance, Next: Random Nu The following table shows the relative performance of a selection the available random number generators. The fastest simulation quality -generators are `taus', `gfsr4' and `mt19937'. The generators which +generators are 'taus', 'gfsr4' and 'mt19937'. The generators which offer the best mathematically-proven quality are those based on the RANLUX algorithm. @@ -4908,11 +4872,11 @@ Here is the output of the program, 0.73995 The numbers depend on the seed used by the generator. The default seed -can be changed with the `GSL_RNG_SEED' environment variable to produce -a different stream of numbers. The generator itself can be changed -using the environment variable `GSL_RNG_TYPE'. Here is the output of -the program using a seed value of 123 and the multiple-recursive -generator `mrg', +can be changed with the 'GSL_RNG_SEED' environment variable to produce a +different stream of numbers. The generator itself can be changed using +the environment variable 'GSL_RNG_TYPE'. Here is the output of the +program using a seed value of 123 and the multiple-recursive generator +'mrg', $ GSL_RNG_SEED=123 GSL_RNG_TYPE=mrg ./a.out GSL_RNG_TYPE=mrg @@ -4935,9 +4899,9 @@ File: gsl-ref.info, Node: Random Number References and Further Reading, Next: ==================================== The subject of random number generation and testing is reviewed -extensively in Knuth's `Seminumerical Algorithms'. +extensively in Knuth's 'Seminumerical Algorithms'. - Donald E. Knuth, `The Art of Computer Programming: Seminumerical + Donald E. Knuth, 'The Art of Computer Programming: Seminumerical Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842. Further information is available in the review paper written by Pierre @@ -4946,15 +4910,14 @@ L'Ecuyer, P. L'Ecuyer, "Random Number Generation", Chapter 4 of the Handbook on Simulation, Jerry Banks Ed., Wiley, 1998, 93-137. - `http://www.iro.umontreal.ca/~lecuyer/papers.html' in the file - `handsim.ps'. + <http://www.iro.umontreal.ca/~lecuyer/papers.html> in the file + 'handsim.ps'. The source code for the DIEHARD random number generator tests is also available online, - `DIEHARD source code' G. Marsaglia, - - `http://stat.fsu.edu/pub/diehard/' + 'DIEHARD source code' G. Marsaglia, + <http://stat.fsu.edu/pub/diehard/> A comprehensive set of random number generator tests is available from NIST, @@ -4962,8 +4925,7 @@ NIST, NIST Special Publication 800-22, "A Statistical Test Suite for the Validation of Random Number Generators and Pseudo Random Number Generators for Cryptographic Applications". - - `http://csrc.nist.gov/rng/' + <http://csrc.nist.gov/rng/> File: gsl-ref.info, Node: Random Number Acknowledgements, Prev: Random Number References and Further Reading, Up: Random Number Generation @@ -4985,24 +4947,24 @@ File: gsl-ref.info, Node: Quasi-Random Sequences, Next: Random Number Distribu This chapter describes functions for generating quasi-random sequences in arbitrary dimensions. A quasi-random sequence progressively covers a -d-dimensional space with a set of points that are uniformly -distributed. Quasi-random sequences are also known as low-discrepancy -sequences. The quasi-random sequence generators use an interface that -is similar to the interface for random number generators, except that -seeding is not required--each generator produces a single sequence. +d-dimensional space with a set of points that are uniformly distributed. +Quasi-random sequences are also known as low-discrepancy sequences. The +quasi-random sequence generators use an interface that is similar to the +interface for random number generators, except that seeding is not +required--each generator produces a single sequence. The functions described in this section are declared in the header -file `gsl_qrng.h'. +file 'gsl_qrng.h'. * Menu: -* Quasi-random number generator initialization:: -* Sampling from a quasi-random number generator:: -* Auxiliary quasi-random number generator functions:: -* Saving and resorting quasi-random number generator state:: -* Quasi-random number generator algorithms:: -* Quasi-random number generator examples:: -* Quasi-random number references:: +* Quasi-random number generator initialization:: +* Sampling from a quasi-random number generator:: +* Auxiliary quasi-random number generator functions:: +* Saving and resorting quasi-random number generator state:: +* Quasi-random number generator algorithms:: +* Quasi-random number generator examples:: +* Quasi-random number references:: File: gsl-ref.info, Node: Quasi-random number generator initialization, Next: Sampling from a quasi-random number generator, Up: Quasi-Random Sequences @@ -5016,7 +4978,7 @@ File: gsl-ref.info, Node: Quasi-random number generator initialization, Next: quasi-random sequence generator of type T and dimension D. If there is insufficient memory to create the generator then the function returns a null pointer and the error handler is invoked - with an error code of `GSL_ENOMEM'. + with an error code of 'GSL_ENOMEM'. -- Function: void gsl_qrng_free (gsl_qrng * Q) This function frees all the memory associated with the generator Q. @@ -5034,9 +4996,9 @@ File: gsl-ref.info, Node: Sampling from a quasi-random number generator, Next: -- Function: int gsl_qrng_get (const gsl_qrng * Q, double X[]) This function stores the next point from the sequence generator Q - in the array X. The space available for X must match the - dimension of the generator. The point X will lie in the range 0 < - x_i < 1 for each x_i. + in the array X. The space available for X must match the dimension + of the generator. The point X will lie in the range 0 < x_i < 1 + for each x_i. File: gsl-ref.info, Node: Auxiliary quasi-random number generator functions, Next: Saving and resorting quasi-random number generator state, Prev: Sampling from a quasi-random number generator, Up: Quasi-Random Sequences @@ -5084,12 +5046,12 @@ The following quasi-random sequence algorithms are available, -- Generator: gsl_qrng_niederreiter_2 This generator uses the algorithm described in Bratley, Fox, - Niederreiter, `ACM Trans. Model. Comp. Sim.' 2, 195 (1992). It is + Niederreiter, 'ACM Trans. Model. Comp. Sim.' 2, 195 (1992). It is valid up to 12 dimensions. -- Generator: gsl_qrng_sobol This generator uses the Sobol sequence described in Antonov, - Saleev, `USSR Comput. Maths. Math. Phys.' 19, 252 (1980). It is + Saleev, 'USSR Comput. Maths. Math. Phys.' 19, 252 (1980). It is valid up to 40 dimensions. @@ -5147,7 +5109,7 @@ the algorithms described in the following paper, P. Bratley and B.L. Fox and H. Niederreiter, "Algorithm 738: Programs to Generate Niederreiter's Low-discrepancy Sequences", - `ACM Transactions on Mathematical Software', Vol. 20, No. 4, + 'ACM Transactions on Mathematical Software', Vol. 20, No. 4, December, 1994, p. 494-495. @@ -5177,59 +5139,58 @@ inverses are computed separately for the upper and lower tails of the distribution, allowing full accuracy to be retained for small results. The functions for random variates and probability density functions -described in this section are declared in `gsl_randist.h'. The +described in this section are declared in 'gsl_randist.h'. The corresponding cumulative distribution functions are declared in -`gsl_cdf.h'. +'gsl_cdf.h'. - Note that the discrete random variate functions always return a -value of type `unsigned int', and on most platforms this has a maximum -value of 2^32-1 ~=~ 4.29e9. They should only be called with a safe -range of parameters (where there is a negligible probability of a -variate exceeding this limit) to prevent incorrect results due to -overflow. + Note that the discrete random variate functions always return a value +of type 'unsigned int', and on most platforms this has a maximum value +of 2^32-1 ~=~ 4.29e9. They should only be called with a safe range of +parameters (where there is a negligible probability of a variate +exceeding this limit) to prevent incorrect results due to overflow. * Menu: -* Random Number Distribution Introduction:: -* The Gaussian Distribution:: -* The Gaussian Tail Distribution:: -* The Bivariate Gaussian Distribution:: -* The Exponential Distribution:: -* The Laplace Distribution:: -* The Exponential Power Distribution:: -* The Cauchy Distribution:: -* The Rayleigh Distribution:: -* The Rayleigh Tail Distribution:: -* The Landau Distribution:: -* The Levy alpha-Stable Distributions:: -* The Levy skew alpha-Stable Distribution:: -* The Gamma Distribution:: -* The Flat (Uniform) Distribution:: -* The Lognormal Distribution:: -* The Chi-squared Distribution:: -* The F-distribution:: -* The t-distribution:: -* The Beta Distribution:: -* The Logistic Distribution:: -* The Pareto Distribution:: -* Spherical Vector Distributions:: -* The Weibull Distribution:: -* The Type-1 Gumbel Distribution:: -* The Type-2 Gumbel Distribution:: -* The Dirichlet Distribution:: -* General Discrete Distributions:: -* The Poisson Distribution:: -* The Bernoulli Distribution:: -* The Binomial Distribution:: -* The Multinomial Distribution:: -* The Negative Binomial Distribution:: -* The Pascal Distribution:: -* The Geometric Distribution:: -* The Hypergeometric Distribution:: -* The Logarithmic Distribution:: -* Shuffling and Sampling:: -* Random Number Distribution Examples:: -* Random Number Distribution References and Further Reading:: +* Random Number Distribution Introduction:: +* The Gaussian Distribution:: +* The Gaussian Tail Distribution:: +* The Bivariate Gaussian Distribution:: +* The Exponential Distribution:: +* The Laplace Distribution:: +* The Exponential Power Distribution:: +* The Cauchy Distribution:: +* The Rayleigh Distribution:: +* The Rayleigh Tail Distribution:: +* The Landau Distribution:: +* The Levy alpha-Stable Distributions:: +* The Levy skew alpha-Stable Distribution:: +* The Gamma Distribution:: +* The Flat (Uniform) Distribution:: +* The Lognormal Distribution:: +* The Chi-squared Distribution:: +* The F-distribution:: +* The t-distribution:: +* The Beta Distribution:: +* The Logistic Distribution:: +* The Pareto Distribution:: +* Spherical Vector Distributions:: +* The Weibull Distribution:: +* The Type-1 Gumbel Distribution:: +* The Type-2 Gumbel Distribution:: +* The Dirichlet Distribution:: +* General Discrete Distributions:: +* The Poisson Distribution:: +* The Bernoulli Distribution:: +* The Binomial Distribution:: +* The Multinomial Distribution:: +* The Negative Binomial Distribution:: +* The Pascal Distribution:: +* The Geometric Distribution:: +* The Hypergeometric Distribution:: +* The Logarithmic Distribution:: +* Shuffling and Sampling:: +* Random Number Distribution Examples:: +* Random Number Distribution References and Further Reading:: File: gsl-ref.info, Node: Random Number Distribution Introduction, Next: The Gaussian Distribution, Up: Random Number Distributions @@ -5269,8 +5230,8 @@ defined as, P(k) = \sum_{i <= k} p(i) -where the sum is over the allowed range of the distribution less than -or equal to k. +where the sum is over the allowed range of the distribution less than or +equal to k. The cumulative distribution for the upper tail of a discrete distribution Q(k) is defined as @@ -5297,10 +5258,10 @@ File: gsl-ref.info, Node: The Gaussian Distribution, Next: The Gaussian Tail D p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx for x in the range -\infty to +\infty. Use the transformation z = - \mu + x on the numbers returned by `gsl_ran_gaussian' to obtain a + \mu + x on the numbers returned by 'gsl_ran_gaussian' to obtain a Gaussian distribution with mean \mu. This function uses the - Box-Mueller algorithm which requires two calls to the random - number generator R. + Box-Mueller algorithm which requires two calls to the random number + generator R. -- Function: double gsl_ran_gaussian_pdf (double X, double SIGMA) This function computes the probability density p(x) at X for a @@ -5320,24 +5281,24 @@ File: gsl-ref.info, Node: The Gaussian Distribution, Next: The Gaussian Tail D -- Function: double gsl_ran_ugaussian (const gsl_rng * R) -- Function: double gsl_ran_ugaussian_pdf (double X) -- Function: double gsl_ran_ugaussian_ratio_method (const gsl_rng * R) - These functions compute results for the unit Gaussian - distribution. They are equivalent to the functions above with a - standard deviation of one, SIGMA = 1. + These functions compute results for the unit Gaussian distribution. + They are equivalent to the functions above with a standard + deviation of one, SIGMA = 1. -- Function: double gsl_cdf_gaussian_P (double X, double SIGMA) -- Function: double gsl_cdf_gaussian_Q (double X, double SIGMA) -- Function: double gsl_cdf_gaussian_Pinv (double P, double SIGMA) -- Function: double gsl_cdf_gaussian_Qinv (double Q, double SIGMA) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Gaussian distribution with - standard deviation SIGMA. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Gaussian distribution with standard + deviation SIGMA. -- Function: double gsl_cdf_ugaussian_P (double X) -- Function: double gsl_cdf_ugaussian_Q (double X) -- Function: double gsl_cdf_ugaussian_Pinv (double P) -- Function: double gsl_cdf_ugaussian_Qinv (double Q) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the unit Gaussian distribution. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the unit Gaussian distribution. File: gsl-ref.info, Node: The Gaussian Tail Distribution, Next: The Bivariate Gaussian Distribution, Prev: The Gaussian Distribution, Up: Random Number Distributions @@ -5345,15 +5306,14 @@ File: gsl-ref.info, Node: The Gaussian Tail Distribution, Next: The Bivariate 19.3 The Gaussian Tail Distribution =================================== - -- Function: double gsl_ran_gaussian_tail (const gsl_rng * R, double - A, double SIGMA) + -- Function: double gsl_ran_gaussian_tail (const gsl_rng * R, double A, + double SIGMA) This function provides random variates from the upper tail of a Gaussian distribution with standard deviation SIGMA. The values - returned are larger than the lower limit A, which must be - positive. The method is based on Marsaglia's famous - rectangle-wedge-tail algorithm (Ann. Math. Stat. 32, 894-899 - (1961)), with this aspect explained in Knuth, v2, 3rd ed, p139,586 - (exercise 11). + returned are larger than the lower limit A, which must be positive. + The method is based on Marsaglia's famous rectangle-wedge-tail + algorithm (Ann. Math. Stat. 32, 894-899 (1961)), with this + aspect explained in Knuth, v2, 3rd ed, p139,586 (exercise 11). The probability distribution for Gaussian tail random variates is, @@ -5363,7 +5323,6 @@ File: gsl-ref.info, Node: The Gaussian Tail Distribution, Next: The Bivariate N(a;\sigma) = (1/2) erfc(a / sqrt(2 sigma^2)). - -- Function: double gsl_ran_gaussian_tail_pdf (double X, double A, double SIGMA) This function computes the probability density p(x) at X for a @@ -5384,9 +5343,8 @@ File: gsl-ref.info, Node: The Bivariate Gaussian Distribution, Next: The Expon 19.4 The Bivariate Gaussian Distribution ======================================== - -- Function: void gsl_ran_bivariate_gaussian (const gsl_rng * R, - double SIGMA_X, double SIGMA_Y, double RHO, double * X, - double * Y) + -- Function: void gsl_ran_bivariate_gaussian (const gsl_rng * R, double + SIGMA_X, double SIGMA_Y, double RHO, double * X, double * Y) This function generates a pair of correlated Gaussian variates, with mean zero, correlation coefficient RHO and standard deviations SIGMA_X and SIGMA_Y in the x and y directions. The probability @@ -5397,12 +5355,12 @@ File: gsl-ref.info, Node: The Bivariate Gaussian Distribution, Next: The Expon for x,y in the range -\infty to +\infty. The correlation coefficient RHO should lie between 1 and -1. - -- Function: double gsl_ran_bivariate_gaussian_pdf (double X, double - Y, double SIGMA_X, double SIGMA_Y, double RHO) + -- Function: double gsl_ran_bivariate_gaussian_pdf (double X, double Y, + double SIGMA_X, double SIGMA_Y, double RHO) This function computes the probability density p(x,y) at (X,Y) for - a bivariate Gaussian distribution with standard deviations - SIGMA_X, SIGMA_Y and correlation coefficient RHO, using the - formula given above. + a bivariate Gaussian distribution with standard deviations SIGMA_X, + SIGMA_Y and correlation coefficient RHO, using the formula given + above. @@ -5413,7 +5371,7 @@ File: gsl-ref.info, Node: The Exponential Distribution, Next: The Laplace Dist -- Function: double gsl_ran_exponential (const gsl_rng * R, double MU) This function returns a random variate from the exponential - distribution with mean MU. The distribution is, + distribution with mean MU. The distribution is, p(x) dx = {1 \over \mu} \exp(-x/\mu) dx @@ -5429,9 +5387,9 @@ File: gsl-ref.info, Node: The Exponential Distribution, Next: The Laplace Dist -- Function: double gsl_cdf_exponential_Q (double X, double MU) -- Function: double gsl_cdf_exponential_Pinv (double P, double MU) -- Function: double gsl_cdf_exponential_Qinv (double Q, double MU) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the exponential distribution - with mean MU. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the exponential distribution with mean + MU. File: gsl-ref.info, Node: The Laplace Distribution, Next: The Exponential Power Distribution, Prev: The Exponential Distribution, Up: Random Number Distributions @@ -5456,9 +5414,8 @@ File: gsl-ref.info, Node: The Laplace Distribution, Next: The Exponential Powe -- Function: double gsl_cdf_laplace_Q (double X, double A) -- Function: double gsl_cdf_laplace_Pinv (double P, double A) -- Function: double gsl_cdf_laplace_Qinv (double Q, double A) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Laplace distribution with - width A. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Laplace distribution with width A. File: gsl-ref.info, Node: The Exponential Power Distribution, Next: The Cauchy Distribution, Prev: The Laplace Distribution, Up: Random Number Distributions @@ -5466,8 +5423,8 @@ File: gsl-ref.info, Node: The Exponential Power Distribution, Next: The Cauchy 19.7 The Exponential Power Distribution ======================================= - -- Function: double gsl_ran_exppow (const gsl_rng * R, double A, - double B) + -- Function: double gsl_ran_exppow (const gsl_rng * R, double A, double + B) This function returns a random variate from the exponential power distribution with scale parameter A and exponent B. The distribution is, @@ -5475,8 +5432,8 @@ File: gsl-ref.info, Node: The Exponential Power Distribution, Next: The Cauchy p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx for x >= 0. For b = 1 this reduces to the Laplace distribution. - For b = 2 it has the same form as a gaussian distribution, but - with a = \sqrt{2} \sigma. + For b = 2 it has the same form as a gaussian distribution, but with + a = \sqrt{2} \sigma. -- Function: double gsl_ran_exppow_pdf (double X, double A, double B) This function computes the probability density p(x) at X for an @@ -5486,9 +5443,9 @@ File: gsl-ref.info, Node: The Exponential Power Distribution, Next: The Cauchy -- Function: double gsl_cdf_exppow_P (double X, double A, double B) -- Function: double gsl_cdf_exppow_Q (double X, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) for the exponential power distribution with parameters - A and B. + These functions compute the cumulative distribution functions P(x), + Q(x) for the exponential power distribution with parameters A and + B. File: gsl-ref.info, Node: The Cauchy Distribution, Next: The Rayleigh Distribution, Prev: The Exponential Power Distribution, Up: Random Number Distributions @@ -5497,9 +5454,9 @@ File: gsl-ref.info, Node: The Cauchy Distribution, Next: The Rayleigh Distribu ============================ -- Function: double gsl_ran_cauchy (const gsl_rng * R, double A) - This function returns a random variate from the Cauchy - distribution with scale parameter A. The probability distribution - for Cauchy random variates is, + This function returns a random variate from the Cauchy distribution + with scale parameter A. The probability distribution for Cauchy + random variates is, p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx @@ -5508,17 +5465,17 @@ File: gsl-ref.info, Node: The Cauchy Distribution, Next: The Rayleigh Distribu -- Function: double gsl_ran_cauchy_pdf (double X, double A) This function computes the probability density p(x) at X for a - Cauchy distribution with scale parameter A, using the formula - given above. + Cauchy distribution with scale parameter A, using the formula given + above. -- Function: double gsl_cdf_cauchy_P (double X, double A) -- Function: double gsl_cdf_cauchy_Q (double X, double A) -- Function: double gsl_cdf_cauchy_Pinv (double P, double A) -- Function: double gsl_cdf_cauchy_Qinv (double Q, double A) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Cauchy distribution with - scale parameter A. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Cauchy distribution with scale + parameter A. File: gsl-ref.info, Node: The Rayleigh Distribution, Next: The Rayleigh Tail Distribution, Prev: The Cauchy Distribution, Up: Random Number Distributions @@ -5536,17 +5493,17 @@ File: gsl-ref.info, Node: The Rayleigh Distribution, Next: The Rayleigh Tail D -- Function: double gsl_ran_rayleigh_pdf (double X, double SIGMA) This function computes the probability density p(x) at X for a - Rayleigh distribution with scale parameter SIGMA, using the - formula given above. + Rayleigh distribution with scale parameter SIGMA, using the formula + given above. -- Function: double gsl_cdf_rayleigh_P (double X, double SIGMA) -- Function: double gsl_cdf_rayleigh_Q (double X, double SIGMA) -- Function: double gsl_cdf_rayleigh_Pinv (double P, double SIGMA) -- Function: double gsl_cdf_rayleigh_Qinv (double Q, double SIGMA) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Rayleigh distribution with - scale parameter SIGMA. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Rayleigh distribution with scale + parameter SIGMA. File: gsl-ref.info, Node: The Rayleigh Tail Distribution, Next: The Landau Distribution, Prev: The Rayleigh Distribution, Up: Random Number Distributions @@ -5554,8 +5511,8 @@ File: gsl-ref.info, Node: The Rayleigh Tail Distribution, Next: The Landau Dis 19.10 The Rayleigh Tail Distribution ==================================== - -- Function: double gsl_ran_rayleigh_tail (const gsl_rng * R, double - A, double SIGMA) + -- Function: double gsl_ran_rayleigh_tail (const gsl_rng * R, double A, + double SIGMA) This function returns a random variate from the tail of the Rayleigh distribution with scale parameter SIGMA and a lower limit of A. The distribution is, @@ -5609,10 +5566,10 @@ File: gsl-ref.info, Node: The Levy alpha-Stable Distributions, Next: The Levy p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^alpha) There is no explicit solution for the form of p(x) and the library - does not define a corresponding `pdf' function. For \alpha = 1 - the distribution reduces to the Cauchy distribution. For \alpha = - 2 it is a Gaussian distribution with \sigma = \sqrt{2} c. For - \alpha < 1 the tails of the distribution become extremely wide. + does not define a corresponding 'pdf' function. For \alpha = 1 the + distribution reduces to the Cauchy distribution. For \alpha = 2 it + is a Gaussian distribution with \sigma = \sqrt{2} c. For \alpha < + 1 the tails of the distribution become extremely wide. The algorithm only works for 0 < alpha <= 2. @@ -5635,20 +5592,19 @@ File: gsl-ref.info, Node: The Levy skew alpha-Stable Distribution, Next: The G When \alpha = 1 the term \tan(\pi \alpha/2) is replaced by -(2/\pi)\log|t|. There is no explicit solution for the form of - p(x) and the library does not define a corresponding `pdf' + p(x) and the library does not define a corresponding 'pdf' function. For \alpha = 2 the distribution reduces to a Gaussian distribution with \sigma = \sqrt{2} c and the skewness parameter - has no effect. For \alpha < 1 the tails of the distribution - become extremely wide. The symmetric distribution corresponds to - \beta = 0. + has no effect. For \alpha < 1 the tails of the distribution become + extremely wide. The symmetric distribution corresponds to \beta = + 0. The algorithm only works for 0 < alpha <= 2. The Levy alpha-stable distributions have the property that if N alpha-stable variates are drawn from the distribution p(c, \alpha, -\beta) then the sum Y = X_1 + X_2 + \dots + X_N will also be -distributed as an alpha-stable variate, p(N^(1/\alpha) c, \alpha, -\beta). +\beta) then the sum Y = X_1 + X_2 + \dots + X_N will also be distributed +as an alpha-stable variate, p(N^(1/\alpha) c, \alpha, \beta). @@ -5659,8 +5615,8 @@ File: gsl-ref.info, Node: The Gamma Distribution, Next: The Flat (Uniform) Dis -- Function: double gsl_ran_gamma (const gsl_rng * R, double A, double B) - This function returns a random variate from the gamma - distribution. The distribution function is, + This function returns a random variate from the gamma distribution. + The distribution function is, p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx @@ -5671,7 +5627,7 @@ File: gsl-ref.info, Node: The Gamma Distribution, Next: The Flat (Uniform) Dis The variates are computed using the Marsaglia-Tsang fast gamma method. This function for this method was previously called - `gsl_ran_gamma_mt' and can still be accessed using this name. + 'gsl_ran_gamma_mt' and can still be accessed using this name. -- Function: double gsl_ran_gamma_knuth (const gsl_rng * R, double A, double B) @@ -5680,17 +5636,17 @@ File: gsl-ref.info, Node: The Gamma Distribution, Next: The Flat (Uniform) Dis -- Function: double gsl_ran_gamma_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a - gamma distribution with parameters A and B, using the formula - given above. + gamma distribution with parameters A and B, using the formula given + above. -- Function: double gsl_cdf_gamma_P (double X, double A, double B) -- Function: double gsl_cdf_gamma_Q (double X, double A, double B) -- Function: double gsl_cdf_gamma_Pinv (double P, double A, double B) -- Function: double gsl_cdf_gamma_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the gamma distribution with - parameters A and B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the gamma distribution with parameters + A and B. File: gsl-ref.info, Node: The Flat (Uniform) Distribution, Next: The Lognormal Distribution, Prev: The Gamma Distribution, Up: Random Number Distributions @@ -5701,7 +5657,7 @@ File: gsl-ref.info, Node: The Flat (Uniform) Distribution, Next: The Lognormal -- Function: double gsl_ran_flat (const gsl_rng * R, double A, double B) This function returns a random variate from the flat (uniform) - distribution from A to B. The distribution is, + distribution from A to B. The distribution is, p(x) dx = {1 \over (b-a)} dx @@ -5716,9 +5672,8 @@ File: gsl-ref.info, Node: The Flat (Uniform) Distribution, Next: The Lognormal -- Function: double gsl_cdf_flat_Q (double X, double A, double B) -- Function: double gsl_cdf_flat_Pinv (double P, double A, double B) -- Function: double gsl_cdf_flat_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for a uniform distribution from A to - B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for a uniform distribution from A to B. File: gsl-ref.info, Node: The Lognormal Distribution, Next: The Chi-squared Distribution, Prev: The Flat (Uniform) Distribution, Up: Random Number Distributions @@ -5750,8 +5705,8 @@ File: gsl-ref.info, Node: The Lognormal Distribution, Next: The Chi-squared Di double SIGMA) -- Function: double gsl_cdf_lognormal_Qinv (double Q, double ZETA, double SIGMA) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the lognormal distribution with + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the lognormal distribution with parameters ZETA and SIGMA. @@ -5770,7 +5725,7 @@ has a chi-squared distribution with n degrees of freedom. -- Function: double gsl_ran_chisq (const gsl_rng * R, double NU) This function returns a random variate from the chi-squared - distribution with NU degrees of freedom. The distribution function + distribution with NU degrees of freedom. The distribution function is, p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx @@ -5787,9 +5742,9 @@ has a chi-squared distribution with n degrees of freedom. -- Function: double gsl_cdf_chisq_Q (double X, double NU) -- Function: double gsl_cdf_chisq_Pinv (double P, double NU) -- Function: double gsl_cdf_chisq_Qinv (double Q, double NU) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the chi-squared distribution - with NU degrees of freedom. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the chi-squared distribution with NU + degrees of freedom. File: gsl-ref.info, Node: The F-distribution, Next: The t-distribution, Prev: The Chi-squared Distribution, Up: Random Number Distributions @@ -5797,9 +5752,8 @@ File: gsl-ref.info, Node: The F-distribution, Next: The t-distribution, Prev: 19.18 The F-distribution ======================== -The F-distribution arises in statistics. If Y_1 and Y_2 are -chi-squared deviates with \nu_1 and \nu_2 degrees of freedom then the -ratio, +The F-distribution arises in statistics. If Y_1 and Y_2 are chi-squared +deviates with \nu_1 and \nu_2 degrees of freedom then the ratio, X = { (Y_1 / \nu_1) \over (Y_2 / \nu_2) } @@ -5807,8 +5761,8 @@ has an F-distribution F(x;\nu_1,\nu_2). -- Function: double gsl_ran_fdist (const gsl_rng * R, double NU1, double NU2) - This function returns a random variate from the F-distribution - with degrees of freedom NU1 and NU2. The distribution function is, + This function returns a random variate from the F-distribution with + degrees of freedom NU1 and NU2. The distribution function is, p(x) dx = { \Gamma((\nu_1 + \nu_2)/2) @@ -5831,9 +5785,9 @@ has an F-distribution F(x;\nu_1,\nu_2). NU2) -- Function: double gsl_cdf_fdist_Qinv (double Q, double NU1, double NU2) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the F-distribution with NU1 and - NU2 degrees of freedom. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the F-distribution with NU1 and NU2 + degrees of freedom. File: gsl-ref.info, Node: The t-distribution, Next: The Beta Distribution, Prev: The F-distribution, Up: Random Number Distributions @@ -5868,9 +5822,9 @@ has a t-distribution t(x;\nu) with \nu degrees of freedom. -- Function: double gsl_cdf_tdist_Q (double X, double NU) -- Function: double gsl_cdf_tdist_Pinv (double P, double NU) -- Function: double gsl_cdf_tdist_Qinv (double Q, double NU) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the t-distribution with NU - degrees of freedom. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the t-distribution with NU degrees of + freedom. File: gsl-ref.info, Node: The Beta Distribution, Next: The Logistic Distribution, Prev: The t-distribution, Up: Random Number Distributions @@ -5888,8 +5842,8 @@ File: gsl-ref.info, Node: The Beta Distribution, Next: The Logistic Distributi for 0 <= x <= 1. -- Function: double gsl_ran_beta_pdf (double X, double A, double B) - This function computes the probability density p(x) at X for a - beta distribution with parameters A and B, using the formula given + This function computes the probability density p(x) at X for a beta + distribution with parameters A and B, using the formula given above. @@ -5897,9 +5851,9 @@ File: gsl-ref.info, Node: The Beta Distribution, Next: The Logistic Distributi -- Function: double gsl_cdf_beta_Q (double X, double A, double B) -- Function: double gsl_cdf_beta_Pinv (double P, double A, double B) -- Function: double gsl_cdf_beta_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the beta distribution with - parameters A and B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the beta distribution with parameters A + and B. File: gsl-ref.info, Node: The Logistic Distribution, Next: The Pareto Distribution, Prev: The Beta Distribution, Up: Random Number Distributions @@ -5925,9 +5879,9 @@ File: gsl-ref.info, Node: The Logistic Distribution, Next: The Pareto Distribu -- Function: double gsl_cdf_logistic_Q (double X, double A) -- Function: double gsl_cdf_logistic_Pinv (double P, double A) -- Function: double gsl_cdf_logistic_Qinv (double Q, double A) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the logistic distribution with - scale parameter A. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the logistic distribution with scale + parameter A. File: gsl-ref.info, Node: The Pareto Distribution, Next: Spherical Vector Distributions, Prev: The Logistic Distribution, Up: Random Number Distributions @@ -5935,10 +5889,10 @@ File: gsl-ref.info, Node: The Pareto Distribution, Next: Spherical Vector Dist 19.22 The Pareto Distribution ============================= - -- Function: double gsl_ran_pareto (const gsl_rng * R, double A, - double B) - This function returns a random variate from the Pareto - distribution of order A. The distribution function is, + -- Function: double gsl_ran_pareto (const gsl_rng * R, double A, double + B) + This function returns a random variate from the Pareto distribution + of order A. The distribution function is, p(x) dx = (a/b) / (x/b)^{a+1} dx @@ -5954,9 +5908,9 @@ File: gsl-ref.info, Node: The Pareto Distribution, Next: Spherical Vector Dist -- Function: double gsl_cdf_pareto_Q (double X, double A, double B) -- Function: double gsl_cdf_pareto_Pinv (double P, double A, double B) -- Function: double gsl_cdf_pareto_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Pareto distribution with - exponent A and scale B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Pareto distribution with exponent A + and scale B. File: gsl-ref.info, Node: Spherical Vector Distributions, Next: The Weibull Distribution, Prev: The Pareto Distribution, Up: Random Number Distributions @@ -5968,14 +5922,14 @@ The spherical distributions generate random vectors, located on a spherical surface. They can be used as random directions, for example in the steps of a random walk. - -- Function: void gsl_ran_dir_2d (const gsl_rng * R, double * X, - double * Y) - -- Function: void gsl_ran_dir_2d_trig_method (const gsl_rng * R, - double * X, double * Y) + -- Function: void gsl_ran_dir_2d (const gsl_rng * R, double * X, double + * Y) + -- Function: void gsl_ran_dir_2d_trig_method (const gsl_rng * R, double + * X, double * Y) This function returns a random direction vector v = (X,Y) in two - dimensions. The vector is normalized such that |v|^2 = x^2 + y^2 - = 1. The obvious way to do this is to take a uniform random - number between 0 and 2\pi and let X and Y be the sine and cosine + dimensions. The vector is normalized such that |v|^2 = x^2 + y^2 = + 1. The obvious way to do this is to take a uniform random number + between 0 and 2\pi and let X and Y be the sine and cosine respectively. Two trig functions would have been expensive in the old days, but with modern hardware implementations, this is sometimes the fastest way to go. This is the case for the Pentium @@ -5989,17 +5943,18 @@ in the steps of a random walk. and V are chosen at random from the interior of a unit circle, and then x=(u^2-v^2)/(u^2+v^2) and y=2uv/(u^2+v^2). - -- Function: void gsl_ran_dir_3d (const gsl_rng * R, double * X, - double * Y, double * Z) + -- Function: void gsl_ran_dir_3d (const gsl_rng * R, double * X, double + * Y, double * Z) This function returns a random direction vector v = (X,Y,Z) in - three dimensions. The vector is normalized such that |v|^2 = x^2 - + y^2 + z^2 = 1. The method employed is due to Robert E. Knop - (CACM 13, 326 (1970)), and explained in Knuth, v2, 3rd ed, p136. - It uses the surprising fact that the distribution projected along - any axis is actually uniform (this is only true for 3 dimensions). + three dimensions. The vector is normalized such that |v|^2 = x^2 + + y^2 + z^2 = 1. The method employed is due to Robert E. Knop (CACM + 13, 326 (1970)), and explained in Knuth, v2, 3rd ed, p136. It uses + the surprising fact that the distribution projected along any axis + is actually uniform (this is only true for 3 dimensions). + + -- Function: void gsl_ran_dir_nd (const gsl_rng * R, size_t N, double * + X) - -- Function: void gsl_ran_dir_nd (const gsl_rng * R, size_t N, double - * X) This function returns a random direction vector v = (x_1,x_2,...,x_n) in N dimensions. The vector is normalized such that |v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1. The method uses the @@ -6026,17 +5981,17 @@ File: gsl-ref.info, Node: The Weibull Distribution, Next: The Type-1 Gumbel Di -- Function: double gsl_ran_weibull_pdf (double X, double A, double B) This function computes the probability density p(x) at X for a - Weibull distribution with scale A and exponent B, using the - formula given above. + Weibull distribution with scale A and exponent B, using the formula + given above. -- Function: double gsl_cdf_weibull_P (double X, double A, double B) -- Function: double gsl_cdf_weibull_Q (double X, double A, double B) -- Function: double gsl_cdf_weibull_Pinv (double P, double A, double B) -- Function: double gsl_cdf_weibull_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Weibull distribution with - scale A and exponent B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Weibull distribution with scale A + and exponent B. File: gsl-ref.info, Node: The Type-1 Gumbel Distribution, Next: The Type-2 Gumbel Distribution, Prev: The Weibull Distribution, Up: Random Number Distributions @@ -6046,7 +6001,7 @@ File: gsl-ref.info, Node: The Type-1 Gumbel Distribution, Next: The Type-2 Gum -- Function: double gsl_ran_gumbel1 (const gsl_rng * R, double A, double B) - This function returns a random variate from the Type-1 Gumbel + This function returns a random variate from the Type-1 Gumbel distribution. The Type-1 Gumbel distribution function is, p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx @@ -6063,9 +6018,9 @@ File: gsl-ref.info, Node: The Type-1 Gumbel Distribution, Next: The Type-2 Gum -- Function: double gsl_cdf_gumbel1_Q (double X, double A, double B) -- Function: double gsl_cdf_gumbel1_Pinv (double P, double A, double B) -- Function: double gsl_cdf_gumbel1_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Type-1 Gumbel distribution - with parameters A and B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Type-1 Gumbel distribution with + parameters A and B. File: gsl-ref.info, Node: The Type-2 Gumbel Distribution, Next: The Dirichlet Distribution, Prev: The Type-1 Gumbel Distribution, Up: Random Number Distributions @@ -6092,9 +6047,9 @@ File: gsl-ref.info, Node: The Type-2 Gumbel Distribution, Next: The Dirichlet -- Function: double gsl_cdf_gumbel2_Q (double X, double A, double B) -- Function: double gsl_cdf_gumbel2_Pinv (double P, double A, double B) -- Function: double gsl_cdf_gumbel2_Qinv (double Q, double A, double B) - These functions compute the cumulative distribution functions - P(x), Q(x) and their inverses for the Type-2 Gumbel distribution - with parameters A and B. + These functions compute the cumulative distribution functions P(x), + Q(x) and their inverses for the Type-2 Gumbel distribution with + parameters A and B. File: gsl-ref.info, Node: The Dirichlet Distribution, Next: General Discrete Distributions, Prev: The Type-2 Gumbel Distribution, Up: Random Number Distributions @@ -6102,22 +6057,22 @@ File: gsl-ref.info, Node: The Dirichlet Distribution, Next: General Discrete D 19.27 The Dirichlet Distribution ================================ - -- Function: void gsl_ran_dirichlet (const gsl_rng * R, size_t K, - const double ALPHA[], double THETA[]) + -- Function: void gsl_ran_dirichlet (const gsl_rng * R, size_t K, const + double ALPHA[], double THETA[]) This function returns an array of K random variates from a - Dirichlet distribution of order K-1. The distribution function is + Dirichlet distribution of order K-1. The distribution function is p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K = (1/Z) \prod_{i=1}^K \theta_i^{\alpha_i - 1} \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 ... d\theta_K - for theta_i >= 0 and alpha_i >= 0. The delta function ensures - that \sum \theta_i = 1. The normalization factor Z is + for theta_i >= 0 and alpha_i >= 0. The delta function ensures that + \sum \theta_i = 1. The normalization factor Z is Z = {\prod_{i=1}^K \Gamma(\alpha_i)} / {\Gamma( \sum_{i=1}^K \alpha_i)} The random variates are generated by sampling K values from gamma distributions with parameters a=alpha_i, b=1, and renormalizing. - See A.M. Law, W.D. Kelton, `Simulation Modeling and Analysis' + See A.M. Law, W.D. Kelton, 'Simulation Modeling and Analysis' (1991). -- Function: double gsl_ran_dirichlet_pdf (size_t K, const double @@ -6148,8 +6103,8 @@ generating a cumulative probability array with K+1 elements: C[k+1] = C[k]+P[k]. Note that this construction produces C[K]=1. Now choose a uniform -deviate u between 0 and 1, and find the value of k such that C[k] <= u -< C[k+1]. Although this in principle requires of order \log K steps per +deviate u between 0 and 1, and find the value of k such that C[k] <= u < +C[k+1]. Although this in principle requires of order \log K steps per random number generation, they are fast steps, and if you use something like \lfloor uK \rfloor as a starting point, you can often do pretty well. @@ -6159,23 +6114,22 @@ preprocess the probability list, and save the result in some form of lookup table; then the individual calls for a random discrete event can go rapidly. An approach invented by G. Marsaglia (Generating discrete random numbers in a computer, Comm ACM 6, 37-38 (1963)) is very clever, -and readers interested in examples of good algorithm design are -directed to this short and well-written paper. Unfortunately, for -large K, Marsaglia's lookup table can be quite large. +and readers interested in examples of good algorithm design are directed +to this short and well-written paper. Unfortunately, for large K, +Marsaglia's lookup table can be quite large. A much better approach is due to Alastair J. Walker (An efficient method for generating discrete random variables with general -distributions, ACM Trans on Mathematical Software 3, 253-256 (1977); -see also Knuth, v2, 3rd ed, p120-121,139). This requires two lookup -tables, one floating point and one integer, but both only of size K. -After preprocessing, the random numbers are generated in O(1) time, -even for large K. The preprocessing suggested by Walker requires -O(K^2) effort, but that is not actually necessary, and the -implementation provided here only takes O(K) effort. In general, more -preprocessing leads to faster generation of the individual random -numbers, but a diminishing return is reached pretty early. Knuth points -out that the optimal preprocessing is combinatorially difficult for -large K. +distributions, ACM Trans on Mathematical Software 3, 253-256 (1977); see +also Knuth, v2, 3rd ed, p120-121,139). This requires two lookup tables, +one floating point and one integer, but both only of size K. After +preprocessing, the random numbers are generated in O(1) time, even for +large K. The preprocessing suggested by Walker requires O(K^2) effort, +but that is not actually necessary, and the implementation provided here +only takes O(K) effort. In general, more preprocessing leads to faster +generation of the individual random numbers, but a diminishing return is +reached pretty early. Knuth points out that the optimal preprocessing +is combinatorially difficult for large K. This method can be used to speed up some of the discrete random number generators below, such as the binomial distribution. To use it @@ -6190,7 +6144,7 @@ to be made, since it only takes a finite set of K outcomes. elements must all be positive, but they needn't add up to one (so you can think of them more generally as "weights")--the preprocessor will normalize appropriately. This return value is - used as an argument for the `gsl_ran_discrete' function below. + used as an argument for the 'gsl_ran_discrete' function below. -- Function: size_t gsl_ran_discrete (const gsl_rng * R, const gsl_ran_discrete_t * G) @@ -6202,8 +6156,8 @@ to be made, since it only takes a finite set of K outcomes. Returns the probability P[k] of observing the variable K. Since P[k] is not stored as part of the lookup table, it must be recomputed; this computation takes O(K), so if K is large and you - care about the original array P[k] used to create the lookup - table, then you should just keep this original array P[k] around. + care about the original array P[k] used to create the lookup table, + then you should just keep this original array P[k] around. -- Function: void gsl_ran_discrete_free (gsl_ran_discrete_t * G) De-allocates the lookup table pointed to by G. @@ -6225,14 +6179,14 @@ File: gsl-ref.info, Node: The Poisson Distribution, Next: The Bernoulli Distri for k >= 0. -- Function: double gsl_ran_poisson_pdf (unsigned int K, double MU) - This function computes the probability p(k) of obtaining K from a + This function computes the probability p(k) of obtaining K from a Poisson distribution with mean MU, using the formula given above. -- Function: double gsl_cdf_poisson_P (unsigned int K, double MU) -- Function: double gsl_cdf_poisson_Q (unsigned int K, double MU) - These functions compute the cumulative distribution functions - P(k), Q(k) for the Poisson distribution with parameter MU. + These functions compute the cumulative distribution functions P(k), + Q(k) for the Poisson distribution with parameter MU. File: gsl-ref.info, Node: The Bernoulli Distribution, Next: The Binomial Distribution, Prev: The Poisson Distribution, Up: Random Number Distributions @@ -6249,7 +6203,6 @@ File: gsl-ref.info, Node: The Bernoulli Distribution, Next: The Binomial Distr p(0) = 1 - p p(1) = p - -- Function: double gsl_ran_bernoulli_pdf (unsigned int K, double P) This function computes the probability p(k) of obtaining K from a Bernoulli distribution with probability parameter P, using the @@ -6284,8 +6237,8 @@ File: gsl-ref.info, Node: The Binomial Distribution, Next: The Multinomial Dis unsigned int N) -- Function: double gsl_cdf_binomial_Q (unsigned int K, double P, unsigned int N) - These functions compute the cumulative distribution functions - P(k), Q(k) for the binomial distribution with parameters P and N. + These functions compute the cumulative distribution functions P(k), + Q(k) for the binomial distribution with parameters P and N. File: gsl-ref.info, Node: The Multinomial Distribution, Next: The Negative Binomial Distribution, Prev: The Binomial Distribution, Up: Random Number Distributions @@ -6295,9 +6248,10 @@ File: gsl-ref.info, Node: The Multinomial Distribution, Next: The Negative Bin -- Function: void gsl_ran_multinomial (const gsl_rng * R, size_t K, unsigned int N, const double P[], unsigned int N[]) + This function computes a random sample N[] from the multinomial distribution formed by N trials from an underlying distribution - P[K]. The distribution function for N[] is, + P[K]. The distribution function for N[] is, P(n_1, n_2, ..., n_K) = (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K @@ -6309,16 +6263,16 @@ File: gsl-ref.info, Node: The Multinomial Distribution, Next: The Negative Bin normalized appropriately. The arrays N[] and P[] must both be of length K. - Random variates are generated using the conditional binomial - method (see C.S. David, `The computer generation of multinomial - random variates', Comp. Stat. Data Anal. 16 (1993) 205-217 for + Random variates are generated using the conditional binomial method + (see C.S. David, 'The computer generation of multinomial random + variates', Comp. Stat. Data Anal. 16 (1993) 205-217 for details). -- Function: double gsl_ran_multinomial_pdf (size_t K, const double P[], const unsigned int N[]) This function computes the probability P(n_1, n_2, ..., n_K) of - sampling N[K] from a multinomial distribution with parameters - P[K], using the formula given above. + sampling N[K] from a multinomial distribution with parameters P[K], + using the formula given above. -- Function: double gsl_ran_multinomial_lnpdf (size_t K, const double P[], const unsigned int N[]) @@ -6332,8 +6286,8 @@ File: gsl-ref.info, Node: The Negative Binomial Distribution, Next: The Pascal 19.33 The Negative Binomial Distribution ======================================== - -- Function: unsigned int gsl_ran_negative_binomial (const gsl_rng * - R, double P, double N) + -- Function: unsigned int gsl_ran_negative_binomial (const gsl_rng * R, + double P, double N) This function returns a random integer from the negative binomial distribution, the number of failures occurring before N successes in independent trials with probability P of success. The @@ -6350,13 +6304,13 @@ File: gsl-ref.info, Node: The Negative Binomial Distribution, Next: The Pascal formula given above. - -- Function: double gsl_cdf_negative_binomial_P (unsigned int K, - double P, double N) - -- Function: double gsl_cdf_negative_binomial_Q (unsigned int K, - double P, double N) - These functions compute the cumulative distribution functions - P(k), Q(k) for the negative binomial distribution with parameters - P and N. + -- Function: double gsl_cdf_negative_binomial_P (unsigned int K, double + P, double N) + -- Function: double gsl_cdf_negative_binomial_Q (unsigned int K, double + P, double N) + These functions compute the cumulative distribution functions P(k), + Q(k) for the negative binomial distribution with parameters P and + N. File: gsl-ref.info, Node: The Pascal Distribution, Next: The Geometric Distribution, Prev: The Negative Binomial Distribution, Up: Random Number Distributions @@ -6385,8 +6339,8 @@ File: gsl-ref.info, Node: The Pascal Distribution, Next: The Geometric Distrib unsigned int N) -- Function: double gsl_cdf_pascal_Q (unsigned int K, double P, unsigned int N) - These functions compute the cumulative distribution functions - P(k), Q(k) for the Pascal distribution with parameters P and N. + These functions compute the cumulative distribution functions P(k), + Q(k) for the Pascal distribution with parameters P and N. File: gsl-ref.info, Node: The Geometric Distribution, Next: The Hypergeometric Distribution, Prev: The Pascal Distribution, Up: Random Number Distributions @@ -6415,8 +6369,8 @@ File: gsl-ref.info, Node: The Geometric Distribution, Next: The Hypergeometric -- Function: double gsl_cdf_geometric_P (unsigned int K, double P) -- Function: double gsl_cdf_geometric_Q (unsigned int K, double P) - These functions compute the cumulative distribution functions - P(k), Q(k) for the geometric distribution with parameter P. + These functions compute the cumulative distribution functions P(k), + Q(k) for the geometric distribution with parameter P. File: gsl-ref.info, Node: The Hypergeometric Distribution, Next: The Logarithmic Distribution, Prev: The Geometric Distribution, Up: Random Number Distributions @@ -6451,9 +6405,9 @@ File: gsl-ref.info, Node: The Hypergeometric Distribution, Next: The Logarithm int N1, unsigned int N2, unsigned int T) -- Function: double gsl_cdf_hypergeometric_Q (unsigned int K, unsigned int N1, unsigned int N2, unsigned int T) - These functions compute the cumulative distribution functions - P(k), Q(k) for the hypergeometric distribution with parameters N1, - N2 and T. + These functions compute the cumulative distribution functions P(k), + Q(k) for the hypergeometric distribution with parameters N1, N2 and + T. File: gsl-ref.info, Node: The Logarithmic Distribution, Next: Shuffling and Sampling, Prev: The Hypergeometric Distribution, Up: Random Number Distributions @@ -6492,10 +6446,11 @@ short period. For more information see Knuth, v2, 3rd ed, Section -- Function: void gsl_ran_shuffle (const gsl_rng * R, void * BASE, size_t N, size_t SIZE) + This function randomly shuffles the order of N objects, each of size SIZE, stored in the array BASE[0..N-1]. The output of the random number generator R is used to produce the permutation. The - algorithm generates all possible n! permutations with equal + algorithm generates all possible n! permutations with equal probability, assuming a perfect source of random numbers. The following code shows how to shuffle the numbers from 0 to 51, @@ -6509,20 +6464,18 @@ short period. For more information see Knuth, v2, 3rd ed, Section gsl_ran_shuffle (r, a, 52, sizeof (int)); - - -- Function: int gsl_ran_choose (const gsl_rng * R, void * DEST, - size_t K, void * SRC, size_t N, size_t SIZE) - This function fills the array DEST[k] with K objects taken - randomly from the N elements of the array SRC[0..N-1]. The - objects are each of size SIZE. The output of the random number - generator R is used to make the selection. The algorithm ensures - all possible samples are equally likely, assuming a perfect source - of randomness. + -- Function: int gsl_ran_choose (const gsl_rng * R, void * DEST, size_t + K, void * SRC, size_t N, size_t SIZE) + This function fills the array DEST[k] with K objects taken randomly + from the N elements of the array SRC[0..N-1]. The objects are each + of size SIZE. The output of the random number generator R is used + to make the selection. The algorithm ensures all possible samples + are equally likely, assuming a perfect source of randomness. The objects are sampled _without_ replacement, thus each object can - only appear once in DEST[k]. It is required that K be less than - or equal to `n'. The objects in DEST will be in the same relative - order as those in SRC. You will need to call `gsl_ran_shuffle(r, + only appear once in DEST[k]. It is required that K be less than or + equal to 'n'. The objects in DEST will be in the same relative + order as those in SRC. You will need to call 'gsl_ran_shuffle(r, dest, n, size)' if you want to randomize the order. The following code shows how to select a random sample of three @@ -6537,13 +6490,12 @@ short period. For more information see Knuth, v2, 3rd ed, Section gsl_ran_choose (r, a, 3, b, 100, sizeof (double)); - -- Function: void gsl_ran_sample (const gsl_rng * R, void * DEST, size_t K, void * SRC, size_t N, size_t SIZE) - This function is like `gsl_ran_choose' but samples K items from - the original array of N items SRC with replacement, so the same - object can appear more than once in the output sequence DEST. - There is no requirement that K be less than N in this case. + This function is like 'gsl_ran_choose' but samples K items from the + original array of N items SRC with replacement, so the same object + can appear more than once in the output sequence DEST. There is no + requirement that K be less than N in this case. File: gsl-ref.info, Node: Random Number Distribution Examples, Next: Random Number Distribution References and Further Reading, Prev: Shuffling and Sampling, Up: Random Number Distributions @@ -6591,7 +6543,7 @@ Poisson distribution with a mean of 3. return 0; } -If the library and header files are installed under `/usr/local' (the +If the library and header files are installed under '/usr/local' (the default location) then the program can be compiled with these options, $ gcc -Wall demo.c -lgsl -lgslcblas -lm @@ -6602,8 +6554,8 @@ Here is the output of the program, 2 5 5 2 1 0 3 4 1 1 The variates depend on the seed used by the generator. The seed for the -default generator type `gsl_rng_default' can be changed with the -`GSL_RNG_SEED' environment variable to produce a different stream of +default generator type 'gsl_rng_default' can be changed with the +'GSL_RNG_SEED' environment variable to produce a different stream of variates, $ GSL_RNG_SEED=123 ./a.out @@ -6642,8 +6594,8 @@ The following program generates a random walk in two dimensions. return 0; } -Here is the output from the program, three 10-step random walks from -the origin, +Here is the output from the program, three 10-step random walks from the +origin, The following program computes the upper and lower cumulative distribution functions for the standard normal distribution at x=2. @@ -6686,49 +6638,49 @@ File: gsl-ref.info, Node: Random Number Distribution References and Further Rea ==================================== For an encyclopaedic coverage of the subject readers are advised to -consult the book `Non-Uniform Random Variate Generation' by Luc -Devroye. It covers every imaginable distribution and provides hundreds -of algorithms. +consult the book 'Non-Uniform Random Variate Generation' by Luc Devroye. +It covers every imaginable distribution and provides hundreds of +algorithms. - Luc Devroye, `Non-Uniform Random Variate Generation', + Luc Devroye, 'Non-Uniform Random Variate Generation', Springer-Verlag, ISBN 0-387-96305-7. The subject of random variate generation is also reviewed by Knuth, who describes algorithms for all the major distributions. - Donald E. Knuth, `The Art of Computer Programming: Seminumerical + Donald E. Knuth, 'The Art of Computer Programming: Seminumerical Algorithms' (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842. The Particle Data Group provides a short review of techniques for generating distributions of random numbers in the "Monte Carlo" section of its Annual Review of Particle Physics. - `Review of Particle Properties' R.M. Barnett et al., Physical - Review D54, 1 (1996) `http://pdg.lbl.gov/'. + 'Review of Particle Properties' R.M. Barnett et al., Physical + Review D54, 1 (1996) <http://pdg.lbl.gov/>. The Review of Particle Physics is available online in postscript and pdf format. An overview of methods used to compute cumulative distribution functions -can be found in `Statistical Computing' by W.J. Kennedy and J.E. -Gentle. Another general reference is `Elements of Statistical -Computing' by R.A. Thisted. +can be found in 'Statistical Computing' by W.J. Kennedy and J.E. Gentle. +Another general reference is 'Elements of Statistical Computing' by R.A. +Thisted. - William E. Kennedy and James E. Gentle, `Statistical Computing' + William E. Kennedy and James E. Gentle, 'Statistical Computing' (1980), Marcel Dekker, ISBN 0-8247-6898-1. - Ronald A. Thisted, `Elements of Statistical Computing' (1988), + Ronald A. Thisted, 'Elements of Statistical Computing' (1988), Chapman & Hall, ISBN 0-412-01371-1. The cumulative distribution functions for the Gaussian distribution are based on the following papers, - `Rational Chebyshev Approximations Using Linear Equations', W.J. - Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242-251 + 'Rational Chebyshev Approximations Using Linear Equations', W.J. + Cody, W. Fraser, J.F. Hart. Numerische Mathematik 12, 242-251 (1968). - `Rational Chebyshev Approximations for the Error Function', W.J. - Cody. Mathematics of Computation 23, n107, 631-637 (July 1969). + 'Rational Chebyshev Approximations for the Error Function', W.J. + Cody. Mathematics of Computation 23, n107, 631-637 (July 1969). File: gsl-ref.info, Node: Statistics, Next: Histograms, Prev: Random Number Distributions, Up: Top @@ -6746,23 +6698,23 @@ intermediate values that might overflow. The functions are available in versions for datasets in the standard floating-point and integer types. The versions for double precision -floating-point data have the prefix `gsl_stats' and are declared in the -header file `gsl_statistics_double.h'. The versions for integer data -have the prefix `gsl_stats_int' and are declared in the header file -`gsl_statistics_int.h'. +floating-point data have the prefix 'gsl_stats' and are declared in the +header file 'gsl_statistics_double.h'. The versions for integer data +have the prefix 'gsl_stats_int' and are declared in the header file +'gsl_statistics_int.h'. * Menu: -* Mean and standard deviation and variance:: -* Absolute deviation:: -* Higher moments (skewness and kurtosis):: -* Autocorrelation:: -* Covariance:: -* Weighted Samples:: -* Maximum and Minimum values:: -* Median and Percentiles:: -* Example statistical programs:: -* Statistics References and Further Reading:: +* Mean and standard deviation and variance:: +* Absolute deviation:: +* Higher moments (skewness and kurtosis):: +* Autocorrelation:: +* Covariance:: +* Weighted Samples:: +* Maximum and Minimum values:: +* Median and Percentiles:: +* Example statistical programs:: +* Statistics References and Further Reading:: File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Absolute deviation, Up: Statistics @@ -6770,8 +6722,8 @@ File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Abso 20.1 Mean, Standard Deviation and Variance ========================================== - -- Function: double gsl_stats_mean (const double DATA[], size_t - STRIDE, size_t N) + -- Function: double gsl_stats_mean (const double DATA[], size_t STRIDE, + size_t N) This function returns the arithmetic mean of DATA, a dataset of length N with stride STRIDE. The arithmetic mean, or "sample mean", is denoted by \Hat\mu and defined as, @@ -6784,9 +6736,9 @@ File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Abso -- Function: double gsl_stats_variance (const double DATA[], size_t STRIDE, size_t N) - This function returns the estimated, or "sample", variance of - DATA, a dataset of length N with stride STRIDE. The estimated - variance is denoted by \Hat\sigma^2 and is defined by, + This function returns the estimated, or "sample", variance of DATA, + a dataset of length N with stride STRIDE. The estimated variance + is denoted by \Hat\sigma^2 and is defined by, \Hat\sigma^2 = (1/(N-1)) \sum (x_i - \Hat\mu)^2 @@ -6796,9 +6748,9 @@ File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Abso \sigma^2. For samples drawn from a gaussian distribution the variance of \Hat\sigma^2 itself is 2 \sigma^4 / N. - This function computes the mean via a call to `gsl_stats_mean'. If + This function computes the mean via a call to 'gsl_stats_mean'. If you have already computed the mean then you can pass it directly to - `gsl_stats_variance_m'. + 'gsl_stats_variance_m'. -- Function: double gsl_stats_variance_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) @@ -6810,26 +6762,26 @@ File: gsl-ref.info, Node: Mean and standard deviation and variance, Next: Abso -- Function: double gsl_stats_sd (const double DATA[], size_t STRIDE, size_t N) - -- Function: double gsl_stats_sd_m (const double DATA[], size_t - STRIDE, size_t N, double MEAN) + -- Function: double gsl_stats_sd_m (const double DATA[], size_t STRIDE, + size_t N, double MEAN) The standard deviation is defined as the square root of the variance. These functions return the square root of the corresponding variance functions above. -- Function: double gsl_stats_variance_with_fixed_mean (const double DATA[], size_t STRIDE, size_t N, double MEAN) - This function computes an unbiased estimate of the variance of - DATA when the population mean MEAN of the underlying distribution - is known _a priori_. In this case the estimator for the variance - uses the factor 1/N and the sample mean \Hat\mu is replaced by the - known population mean \mu, + This function computes an unbiased estimate of the variance of DATA + when the population mean MEAN of the underlying distribution is + known _a priori_. In this case the estimator for the variance uses + the factor 1/N and the sample mean \Hat\mu is replaced by the known + population mean \mu, \Hat\sigma^2 = (1/N) \sum (x_i - \mu)^2 -- Function: double gsl_stats_sd_with_fixed_mean (const double DATA[], size_t STRIDE, size_t N, double MEAN) - This function calculates the standard deviation of DATA for a - fixed population mean MEAN. The result is the square root of the + This function calculates the standard deviation of DATA for a fixed + population mean MEAN. The result is the square root of the corresponding variance function. @@ -6847,9 +6799,9 @@ File: gsl-ref.info, Node: Absolute deviation, Next: Higher moments (skewness a absdev = (1/N) \sum |x_i - \Hat\mu| where x_i are the elements of the dataset DATA. The absolute - deviation from the mean provides a more robust measure of the - width of a distribution than the variance. This function computes - the mean of DATA via a call to `gsl_stats_mean'. + deviation from the mean provides a more robust measure of the width + of a distribution than the variance. This function computes the + mean of DATA via a call to 'gsl_stats_mean'. -- Function: double gsl_stats_absdev_m (const double DATA[], size_t STRIDE, size_t N, double MEAN) @@ -6869,8 +6821,8 @@ File: gsl-ref.info, Node: Higher moments (skewness and kurtosis), Next: Autoco 20.3 Higher moments (skewness and kurtosis) =========================================== - -- Function: double gsl_stats_skew (const double DATA[], size_t - STRIDE, size_t N) + -- Function: double gsl_stats_skew (const double DATA[], size_t STRIDE, + size_t N) This function computes the skewness of DATA, a dataset of length N with stride STRIDE. The skewness is defined as, @@ -6880,7 +6832,7 @@ File: gsl-ref.info, Node: Higher moments (skewness and kurtosis), Next: Autoco measures the asymmetry of the tails of a distribution. The function computes the mean and estimated standard deviation of - DATA via calls to `gsl_stats_mean' and `gsl_stats_sd'. + DATA via calls to 'gsl_stats_mean' and 'gsl_stats_sd'. -- Function: double gsl_stats_skew_m_sd (const double DATA[], size_t STRIDE, size_t N, double MEAN, double SD) @@ -6934,7 +6886,6 @@ File: gsl-ref.info, Node: Autocorrelation, Next: Covariance, Prev: Higher mom This function computes the lag-1 autocorrelation of the dataset DATA using the given value of the mean MEAN. - File: gsl-ref.info, Node: Covariance, Next: Weighted Samples, Prev: Autocorrelation, Up: Statistics @@ -6949,14 +6900,13 @@ File: gsl-ref.info, Node: Covariance, Next: Weighted Samples, Prev: Autocorre covar = (1/(n - 1)) \sum_{i = 1}^{n} (x_i - \Hat x) (y_i - \Hat y) - -- Function: double gsl_stats_covariance_m (const double DATA1[], - const size_t STRIDE1, const double DATA2[], const size_t - STRIDE2, const size_t N, const double MEAN1, const double - MEAN2) + -- Function: double gsl_stats_covariance_m (const double DATA1[], const + size_t STRIDE1, const double DATA2[], const size_t STRIDE2, + const size_t N, const double MEAN1, const double MEAN2) This function computes the covariance of the datasets DATA1 and DATA2 using the given values of the means, MEAN1 and MEAN2. This - is useful if you have already computed the means of DATA1 and - DATA2 and want to avoid recomputing them. + is useful if you have already computed the means of DATA1 and DATA2 + and want to avoid recomputing them. File: gsl-ref.info, Node: Weighted Samples, Next: Maximum and Minimum values, Prev: Covariance, Up: Statistics @@ -6968,9 +6918,9 @@ The functions described in this section allow the computation of statistics for weighted samples. The functions accept an array of samples, x_i, with associated weights, w_i. Each sample x_i is considered as having been drawn from a Gaussian distribution with -variance \sigma_i^2. The sample weight w_i is defined as the -reciprocal of this variance, w_i = 1/\sigma_i^2. Setting a weight to -zero corresponds to removing a sample from a dataset. +variance \sigma_i^2. The sample weight w_i is defined as the reciprocal +of this variance, w_i = 1/\sigma_i^2. Setting a weight to zero +corresponds to removing a sample from a dataset. -- Function: double gsl_stats_wmean (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N) @@ -7004,12 +6954,12 @@ zero corresponds to removing a sample from a dataset. const double DATA[], size_t STRIDE, size_t N) The standard deviation is defined as the square root of the variance. This function returns the square root of the - corresponding variance function `gsl_stats_wvariance' above. + corresponding variance function 'gsl_stats_wvariance' above. -- Function: double gsl_stats_wsd_m (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, size_t N, double WMEAN) This function returns the square root of the corresponding variance - function `gsl_stats_wvariance_m' above. + function 'gsl_stats_wvariance_m' above. -- Function: double gsl_stats_wvariance_with_fixed_mean (const double W[], size_t WSTRIDE, const double DATA[], size_t STRIDE, @@ -7069,3 +7019,61 @@ zero corresponds to removing a sample from a dataset. using the given values of the weighted mean and weighted standard deviation, WMEAN and WSD. + +File: gsl-ref.info, Node: Maximum and Minimum values, Next: Median and Percentiles, Prev: Weighted Samples, Up: Statistics + +20.7 Maximum and Minimum values +=============================== + +The following functions find the maximum and minimum values of a dataset +(or their indices). If the data contains 'NaN's then a 'NaN' will be +returned, since the maximum or minimum value is undefined. For +functions which return an index, the location of the first 'NaN' in the +array is returned. + + -- Function: double gsl_stats_max (const double DATA[], size_t STRIDE, + size_t N) + This function returns the maximum value in DATA, a dataset of + length N with stride STRIDE. The maximum value is defined as the + value of the element x_i which satisfies x_i >= x_j for all j. + + If you want instead to find the element with the largest absolute + magnitude you will need to apply 'fabs' or 'abs' to your data + before calling this function. + + -- Function: double gsl_stats_min (const double DATA[], size_t STRIDE, + size_t N) + This function returns the minimum value in DATA, a dataset of + length N with stride STRIDE. The minimum value is defined as the + value of the element x_i which satisfies x_i <= x_j for all j. + + If you want instead to find the element with the smallest absolute + magnitude you will need to apply 'fabs' or 'abs' to your data + before calling this function. + + -- Function: void gsl_stats_minmax (double * MIN, double * MAX, const + double DATA[], size_t STRIDE, size_t N) + This function finds both the minimum and maximum values MIN, MAX in + DATA in a single pass. + + -- Function: size_t gsl_stats_max_index (const double DATA[], size_t + STRIDE, size_t N) + This function returns the index of the maximum value in DATA, a + dataset of length N with stride STRIDE. The maximum value is + defined as the value of the element x_i which satisfies x_i >= x_j + for all j. When there are several equal maximum elements then the + first one is chosen. + + -- Function: size_t gsl_stats_min_index (const double DATA[], size_t + STRIDE, size_t N) + This function returns the index of the minimum value in DATA, a + dataset of length N with stride STRIDE. The minimum value is + defined as the value of the element x_i which satisfies x_i >= x_j + for all j. When there are several equal minimum elements then the + first one is chosen. + + -- Function: void gsl_stats_minmax_index (size_t * MIN_INDEX, size_t * + MAX_INDEX, const double DATA[], size_t STRIDE, size_t N) + This function returns the indexes MIN_INDEX, MAX_INDEX of the + minimum and maximum values in DATA in a single pass. + |