Numerical Recipes Second Edition (Legacy)
C Routines by Chapter and Section

Users are strongly urged to use Numerical Recipes Third Edition C++ routines rather than the legacy C routines linked here. The files here have not been significantly updated since 1992 and contain many deprecated features, for example one- (not zero-) based arrays, etc.

Click on any link to download the source code file for the recipe or example.

Header and Utility Files (rename to _h to .h before using)

nr_h (nr.h)

nrutil_h (nrutil.h)

nrutil.c

complex_h (complex.h)

complex.c

Chapter 1

[1.0] flmoon calculate phases of the moon by date (example)

[1.1] julday Julian Day number from calendar date (example)

[1.1] badluk Friday the 13th when the moon is full

[1.1] caldat calendar date from Julian day number (example)

Chapter 2

[2.1] gaussj Gauss-Jordan matrix inversion and linear equation solution (example)

[2.3] ludcmp linear equation solution, LU decomposition (example)

[2.3] lubksb linear equation solution, backsubstitution (example)

[2.4] tridag solution of tridiagonal systems (example)

[2.4] banmul multiply vector by band diagonal matrix (example)

[2.4] bandec band diagonal systems, decomposition (example)

[2.4] banbks band diagonal systems, backsubstitution

[2.5] mprove linear equation solution, iterative improvement (example)

[2.6] svbksb singular value backsubstitution (example)

[2.6] svdcmp singular value decomposition of a matrix (example)

[2.6] pythag calculate (a^2+b^2)^{1/2} without overflow

[2.7] cyclic solution of cyclic tridiagonal systems (example)

[2.7] sprsin convert matrix to sparse format (example)

[2.7] sprsax product of sparse matrix and vector (example)

[2.7] sprstx product of transpose sparse matrix and vector (example)

[2.7] sprstp transpose of sparse matrix (example)

[2.7] sprspm pattern multiply two sparse matrices (example)

[2.7] sprstm threshold multiply two sparse matrices (example)

[2.7] linbcg biconjugate gradient solution of sparse systems (example)

[2.7] snrm used by linbcg for vector norm

[2.7] atimes used by linbcg for sparse multiplication

[2.7] asolve used by linbcg for preconditioner

[2.8] vander solve Vandermonde systems (example)

[2.8] toeplz solve Toeplitz systems (example)

[2.9] choldc Cholesky decomposition

[2.9] cholsl Cholesky backsubstitution (example)

[2.10] qrdcmp QR decomposition (example)

[2.10] qrsolv QR backsubstitution (example)

[2.10] rsolv right triangular backsubstitution

[2.10] qrupdt update a QR decomposition (example)

[2.10] rotate Jacobi rotation used by qrupdt

Chapter 3

[3.1] polint polynomial interpolation (example)

[3.2] ratint rational function interpolation (example)

[3.3] spline construct a cubic spline (example)

[3.3] splint cubic spline interpolation (example)

[3.4] locate search an ordered table by bisection (example)

[3.4] hunt search a table when calls are correlated (example)

[3.5] polcoe polynomial coefficients from table of values (example)

[3.5] polcof polynomial coefficients from table of values (example)

[3.6] polin2 two-dimensional polynomial interpolation (example)

[3.6] bcucof construct two-dimensional bicubic (example)

[3.6] bcuint two-dimensional bicubic interpolation (example)

[3.6] splie2 construct two-dimensional spline (example)

[3.6] splin2 two-dimensional spline interpolation (example)

Chapter 4

[4.2] trapzd trapezoidal rule (example)

[4.2] qtrap integrate using trapezoidal rule (example)

[4.2] qsimp integrate using Simpson's rule (example)

[4.3] qromb integrate using Romberg adaptive method (example)

[4.4] midpnt extended midpoint rule (example)

[4.4] qromo integrate using open Romberg adaptive method (example)

[4.4] midinf integrate a function on a semi-infinite interval

[4.4] midsql integrate a function with lower square-root singularity

[4.4] midsqu integrate a function with upper square-root singularity

[4.4] midexp integrate a function that decreases exponentially

[4.5] qgaus integrate a function by Gaussian quadratures (example)

[4.5] gauleg Gauss-Legendre weights and abscissas (example)

[4.5] gaulag Gauss-Laguerre weights and abscissas (example)

[4.5] gauher Gauss-Hermite weights and abscissas (example)

[4.5] gaujac Gauss-Jacobi weights and abscissas (example)

[4.5] gaucof quadrature weights from orthogonal polynomials (example)

[4.5] orthog construct nonclassical orthogonal polynomials (example)

[4.6] quad3d integrate a function over a three-dimensional space (example)

Chapter 5

[5.1] eulsum sum a series by Euler--van Wijngaarden algorithm (example)

[5.3] ddpoly evaluate a polynomial and its derivatives (example)

[5.3] poldiv divide one polynomial by another (example)

[5.3] ratval evaluate a rational function

[5.7] dfridr numerical derivative by Ridders' method (example)

[5.8] chebft fit a Chebyshev polynomial to a function (example)

[5.8] chebev Chebyshev polynomial evaluation (example)

[5.9] chder derivative of a function already Chebyshev fitted (example)

[5.9] chint integrate a function already Chebyshev fitted (example)

[5.10] chebpc polynomial coefficients from a Chebyshev fit (example)

[5.10] pcshft polynomial coefficients of a shifted polynomial (example)

[5.11] pccheb inverse of chebpc; use to economize power series (example)

[5.12] pade Pade approximant from power series coefficients (example)

[5.13] ratlsq rational fit by least-squares method (example)

Chapter 6

[6.1] gammln logarithm of gamma function (example)

[6.1] factrl factorial function (example)

[6.1] bico binomial coefficients function (example)

[6.1] factln logarithm of factorial function (example)

[6.1] beta beta function (example)

[6.2] gammp incomplete gamma function (example)

[6.2] gammq complement of incomplete gamma function (example)

[6.2] gser series used by gammp and gammq (example)

[6.2] gcf continued fraction used by gammp and gammq (example)

[6.2] erf error function

[6.2] erfc complementary error function

[6.2] erfcc complementary error function, concise routine (example)

[6.3] expint exponential integral E_n (example)

[6.3] ei exponential integral Ei (example)

[6.4] betai incomplete beta function (example)

[6.4] betacf continued fraction used by betai

[6.5] bessj0 Bessel function J_0 (example)

[6.5] bessy0 Bessel function Y_0 (example)

[6.5] bessj1 Bessel function J_1 (example)

[6.5] bessy1 Bessel function Y_1 (example)

[6.5] bessy Bessel function Y of general integer order (example)

[6.5] bessj Bessel function J of general integer order (example)

[6.6] bessi0 modified Bessel function I_0 (example)

[6.6] bessk0 modified Bessel function K_0 (example)

[6.6] bessi1 modified Bessel function I_1 (example)

[6.6] bessk1 modified Bessel function K_1 (example)

[6.6] bessk modified Bessel function K of integer order (example)

[6.6] bessi modified Bessel function I of integer order (example)

[6.7] bessjy Bessel functions of fractional order (example)

[6.7] beschb Chebyshev expansion used by bessjy (example)

[6.7] bessik modified Bessel functions of fractional order (example)

[6.7] airy Airy functions

[6.7] sphbes spherical Bessel functions j_n and y_n (example)

[6.8] plgndr Legendre polynomials, associated (spherical harmonics) (example)

[6.9] frenel Fresnel integrals S(x) and C(x) (example)

[6.9] cisi cosine and sine integrals Ci and Si

[6.10] dawson Dawson's integral (example)

[6.11] rf Carlson's elliptic integral of the first kind (example)

[6.11] rd Carlson's elliptic integral of the second kind (example)

[6.11] rj Carlson's elliptic integral of the third kind (example)

[6.11] rc Carlson's degenerate elliptic integral (example)

[6.11] ellf Legendre elliptic integral of the first kind (example)

[6.11] elle Legendre elliptic integral of the second kind (example)

[6.11] ellpi Legendre elliptic integral of the third kind (example)

[6.11] sncndn Jacobian elliptic functions (example)

[6.12] hypgeo complex hypergeometric function (example)

[6.12] hypser complex hypergeometric function, series evaluation

[6.12] hypdrv complex hypergeometric function, derivative of

Chapter 7

[7.1] ran0 random deviate by Park and Miller minimal standard

[7.1] ran1 random deviate, minimal standard plus shuffle

[7.1] ran2 random deviate by L'Ecuyer long period plus shuffle

[7.1] ran3 random deviate by Knuth subtractive method

[7.2] expdev exponential random deviates (example)

[7.2] gasdev normally distributed random deviates (example)

[7.3] gamdev gamma-law distribution random deviates (example)

[7.3] poidev Poisson distributed random deviates (example)

[7.3] bnldev binomial distributed random deviates (example)

[7.4] irbit1 random bit sequence (example)

[7.4] irbit2 random bit sequence (example)

[7.5] psdes ``pseudo-DES'' hashing of 64 bits (example)

[7.5] ran4 random deviates from DES-like hashing (example)

[7.7] sobseq Sobol's quasi-random sequence (example)

[7.8] vegas adaptive multidimensional Monte Carlo integration (example)

[7.8] rebin sample rebinning used by vegas

[7.8] miser recursive multidimensional Monte Carlo integration (example)

[7.8] ranpt get random point, used by miser

Chapter 8

[8.1] piksrt sort an array by straight insertion (example)

[8.1] piksr2 sort two arrays by straight insertion (example)

[8.1] shell sort an array by Shell's method (example)

[8.2] sort sort an array by quicksort method (example)

[8.2] sort2 sort two arrays by quicksort method (example)

[8.3] hpsort sort an array by heapsort method (example)

[8.4] indexx construct an index for an array (example)

[8.4] sort3 sort, use an index to sort 3 or more arrays (example)

[8.4] rank construct a rank table for an array (example)

[8.5] select find the Nth largest in an array (example)

[8.5] selip find the Nth largest, without altering an array (example)

[8.5] hpsel find M largest values, without altering an array (example)

[8.6] eclass determine equivalence classes from list (example)

[8.6] eclazz determine equivalence classes from procedure (example)

Chapter 9

[9.0] scrsho graph a function to search for roots (example)

[9.1] zbrac outward search for brackets on roots (example)

[9.1] zbrak inward search for brackets on roots (example)

[9.1] rtbis find root of a function by bisection (example)

[9.2] rtflsp find root of a function by false-position (example)

[9.2] rtsec find root of a function by secant method (example)

[9.2] zriddr find root of a function by Ridders' method (example)

[9.3] zbrent find root of a function by Brent's method (example)

[9.4] rtnewt find root of a function by Newton-Raphson (example)

[9.4] rtsafe find root of a function by Newton-Raphson and bisection (example)

[9.5] laguer find a root of a polynomial by Laguerre's method (example)

[9.5] zroots roots of a polynomial by Laguerre's method with deflation (example)

[9.5] zrhqr roots of a polynomial by eigenvalue methods (example)

[9.5] qroot complex or double root of a polynomial, Bairstow (example)

[9.6] mnewt Newton's method for systems of equations (example)

[9.7] lnsrch search along a line, used by newt

[9.7] newt globally convergent multi-dimensional Newton's method (example)

[9.7] fdjac finite-difference Jacobian, used by newt

[9.7] fmin norm of a vector function, used by newt

[9.7] broydn secant method for systems of equations (example)

Chapter 10

[10.1] mnbrak bracket the minimum of a function (example)

[10.1] golden find minimum of a function by golden section search (example)

[10.2] brent find minimum of a function by Brent's method (example)

[10.3] dbrent find minimum of a function using derivative information (example)

[10.4] amoeba minimize in N-dimensions by downhill simplex method (example)

[10.4] amotry evaluate a trial point, used by amoeba

[10.5] powell minimize in N-dimensions by Powell's method (example)

[10.5] linmin minimum of a function along a ray in N-dimensions (example)

[10.5] f1dim function used by LINMIN (example)

[10.6] frprmn minimize in N-dimensions by conjugate gradient (example)

[10.6] df1dim alternative function used by LINMIN (example)

[10.7] dfpmin minimize in N-dimensions by variable metric method (example)

[10.8] simplx linear programming maximization of a linear function (example)

[10.8] simp1 linear programming, used by SIMPLX

[10.8] simp2 linear programming, used by SIMPLX

[10.8] simp3 linear programming, used by SIMPLX

[10.9] anneal traveling salesman problem by simulated annealing (example)

[10.9] revcst cost of a reversal, used by anneal

[10.9] reverse do a reversal, used by anneal

[10.9] trncst cost of a transposition, used by anneal

[10.9] trnspt do a transposition, used by anneal

[10.9] metrop Metropolis algorithm, used by anneal

[10.9] amebsa simulated annealing in continuous spaces (example)

[10.9] amotsa evaluate a trial point, used by amebsa

Chapter 11

[11.1] jacobi eigenvalues and eigenvectors of a symmetric matrix (example)

[11.1] eigsrt eigenvectors, sorts into order by eigenvalue (example)

[11.2] tred2 Householder reduction of a real, symmetric matrix (example)

[11.3] tqli eigensolution of a symmetric tridiagonal matrix (example)

[11.5] balanc balance a nonsymmetric matrix (example)

[11.5] elmhes reduce a general matrix to Hessenberg form (example)

[11.6] hqr eigenvalues of a Hessenberg matrix (example)

Chapter 12

[12.2] four1 fast Fourier transform (FFT) in one dimension (example)

[12.3] twofft fast Fourier transform of two real functions (example)

[12.3] realft fast Fourier transform of a single real function (example)

[12.3] sinft fast sine transform (example)

[12.3] cosft1 fast cosine transform with endpoints (example)

[12.3] cosft2 ``staggered'' fast cosine transform (example)

[12.4] fourn fast Fourier transform in multidimensions (example)

[12.5] rlft3 FFT of real data in two or three dimensions (example)

[12.6] fourfs FFT for huge data sets on external media (example)

[12.6] fourew rewind and permute files, used by fourfs

Chapter 13

[13.1] convlv convolution or deconvolution of data using FFT (example)

[13.2] correl correlation or autocorrelation of data using FFT (example)

[13.4] spctrm power spectrum estimation using FFT (example)

[13.6] memcof evaluate maximum entropy (MEM) coefficients (example)

[13.6] fixrts reflect roots of a polynomial into unit circle (example)

[13.6] predic linear prediction using MEM coefficients (example)

[13.7] evlmem power spectral estimation from MEM coefficients (example)

[13.8] period power spectrum of unevenly sampled data (example)

[13.8] fasper power spectrum of unevenly sampled larger data sets (example)

[13.8] spread extirpolate value into array, used by fasper

[13.9] dftcor compute endpoint corrections for Fourier integrals

[13.9] dftint high-accuracy Fourier integrals (example)

[13.10] wt1 one-dimensional discrete wavelet transform

[13.10] daub4 Daubechies 4-coefficient wavelet filter

[13.10] pwtset initialize coefficients for pwt

[13.10] pwt partial wavelet transform

[13.10] wtn multidimensional discrete wavelet transform

Chapter 14

[14.1] moment calculate moments of a data set (example)

[14.2] ttest Student's t-test for difference of means (example)

[14.2] avevar calculate mean and variance of a data set (example)

[14.2] tutest Student's t-test for means, case of unequal variances (example)

[14.2] tptest Student's t-test for means, case of paired data (example)

[14.2] ftest F-test for difference of variances (example)

[14.3] chsone chi-square test for difference between data and model (example)

[14.3] chstwo chi-square test for difference between two data sets (example)

[14.3] ksone Kolmogorov-Smirnov test of data against model (example)

[14.3] kstwo Kolmogorov-Smirnov test between two data sets (example)

[14.3] probks Kolmogorov-Smirnov probability function (example)

[14.4] cntab1 contingency table analysis using chi-square (example)

[14.4] cntab2 contingency table analysis using entropy measure (example)

[14.5] pearsn Pearson's correlation between two data sets (example)

[14.6] spear Spearman's rank correlation between two data sets (example)

[14.6] crank replaces array elements by their rank (example)

[14.6] kendl1 correlation between two data sets, Kendall's tau (example)

[14.6] kendl2 contingency table analysis using Kendall's tau (example)

[14.7] ks2d1s K--S test in two dimensions, data vs. model (example)

[14.7] quadct count points by quadrants, used by ks2d1s

[14.7] quadvl quadrant probabilities, used by ks2d1s

[14.7] ks2d2s K--S test in two dimensions, data vs. data (example)

[14.8] savgol Savitzky-Golay smoothing coefficients (example)

Chapter 15

[15.2] fit least-squares fit data to a straight line (example)

[15.3] fitexy fit data to a straight line, errors in both x and y (example)

[15.3] chixy used by fitexy to calculate a chi^2

[15.4] lfit general linear least-squares fit by normal equations (example)

[15.4] covsrt rearrange covariance matrix, used by LFIT (example)

[15.4] svdfit linear least-squares fit by singular value decomposition (example)

[15.4] svdvar variances from singular value decomposition (example)

[15.4] fpoly fit a polynomial using LFIT or SVDFIT (example)

[15.4] fleg fit a Legendre polynomial using LFIT or SVDFIT (example)

[15.5] mrqmin nonlinear least-squares fit, Marquardt's method (example)

[15.5] mrqcof used by MRQMIN to evaluate coefficients (example)

[15.5] fgauss fit a sum of Gaussians using MRQMIN (example)

[15.7] medfit fit data to a straight line robustly, least absolute deviation (example)

[15.7] rofunc fit data robustly, used by MEDFIT (example)

Chapter 16

[16.1] rk4 integrate one step of ODEs, fourth-order Runge-Kutta (example)

[16.1] rkdumb integrate ODEs by fourth-order Runge-Kutta (example)

[16.2] rkqs integrate one step of ODEs with accuracy monitoring (example)

[16.2] rkck Cash-Karp-Runge-Kutta step used by rkqs

[16.2] odeint integrate ODEs with accuracy monitoring (example)

[16.3] mmid integrate ODEs by modified midpoint method (example)

[16.4] bsstep integrate ODEs, Bulirsch-Stoer step (example)

[16.4] pzextr polynomial extrapolation, used by BSSTEP (example)

[16.4] rzextr rational function extrapolation, used by BSSTEP (example)

[16.5] stoerm integrate conservative second-order ODEs (example)

[16.6] stiff integrate stiff ODEs by fourth-order Rosenbrock (example)

[16.6] jacobn sample Jacobian routine for stiff

[16.6] derivs sample derivatives routine for stiff

[16.6] simpr integrate stiff ODEs by semi-implicit midpoint rule (example)

[16.6] stifbs integrate stiff ODEs, Bulirsch-Stoer step (example)

Chapter 17

[17.1] shoot solve two point boundary value problem by shooting

[17.2] shootf ditto, by shooting to a fitting point

[17.3] solvde two point boundary value problem, solve by relaxation

[17.3] bksub backsubstitution, used by SOLVDE

[17.3] pinvs diagonalize a sub-block, used by SOLVDE

[17.3] red reduce columns of a matrix, used by SOLVDE

[17.4] sfroid spheroidal functions by method of SOLVDE

[17.4] difeq spheroidal matrix coefficients, used by SFROID

[17.4] sphoot spheroidal functions by method of shoot

[17.4] sphfpt spheroidal functions by method of shootf (example)

Chapter 18

[18.1] fred2 solve linear Fredholm equations of the second kind (example)

[18.1] fredin interpolate solutions obtained with fred2 (example)

[18.2] voltra linear Volterra equations of the second kind (example)

[18.3] wwghts quadrature weights for an arbitrarily singular kernel

[18.3] kermom sample routine for moments of a singular kernel

[18.3] quadmx sample routine for a quadrature matrix

[18.3] fredex example of solving a singular Fredholm equation

Chapter 19

[19.5] sor elliptic PDE solved by successive overrelaxation method (example)

[19.6] mglin linear elliptic PDE solved by multigrid method (example)

[19.6] rstrct half-weighting restriction, used by mglin, mgfas

[19.6] interp bilinear prolongation, used by mglin, mgfas

[19.6] addint interpolate and add, used by mglin

[19.6] slvsml solve on coarsest grid, used by mglin

[19.6] relax Gauss-Seidel relaxation, used by mglin

[19.6] resid calculate residual, used by mglin

[19.6] copy utility used by mglin, mgfas

[19.6] fill0 utility used by mglin

[19.6] mgfas nonlinear elliptic PDE solved by multigrid method (example)

[19.6] relax2 Gauss-Seidel relaxation, used by mgfas

[19.6] slvsm2 solve on coarsest grid, used by mgfas

[19.6] lop applies nonlinear operator, used by mgfas

[19.6] matadd utility used by mgfas

[19.6] matsub utility used by mgfas

[19.6] anorm2 utility used by mgfas

Chapter 20

[20.1] machar diagnose computer's floating arithmetic (example)

[20.2] igray Gray code and its inverse (example)

[20.3] icrc1 cyclic redundancy checksum, used by icrc

[20.3] icrc cyclic redundancy checksum (example)

[20.3] decchk decimal check digit calculation or verification (example)

[20.4] hufmak construct a Huffman code

[20.4] hufapp append bits to a Huffman code, used by hufmak

[20.4] hufenc use Huffman code to encode and compress a character

[20.4] hufdec use Huffman code to decode and decompress a character

[20.5] arcmak construct an arithmetic code

[20.5] arcode encode or decode a character using arithmetic coding (example)

[20.5] arcsum add integer to byte string, used by arcode

[20.6] mpops multiple precision arithmetic, simpler operations

[20.6] mpmul multiple precision multiply, using FFT methods

[20.6] mpinv multiple precision reciprocal

[20.6] mpdiv multiple precision divide and remainder

[20.6] mpsqrt multiple precision square root

[20.6] mp2dfr multiple precision conversion to decimal base

[20.6] mppi multiple precision example, compute many digits of pi (example)

Numerical Recipes Second Edition (Legacy) C Routines by Chapter and Section