This is documentation for Mathematica 5, which was
based on an earlier version of the Wolfram Language.
View current documentation (Version 11.2)

Documentation / Mathematica / The Mathematica Book / Mathematica Reference Guide / Some Notes on Internal Implementation /

A.9.4 Numerical and Related Functions

Number representation and numerical evaluation

FilledSmallCircle Large integers and high-precision approximate numbers are stored as arrays of base or digits, depending on the lengths of machine integers.

FilledSmallCircle Precision is internally maintained as a floating-point number.

FilledSmallCircleIntegerDigits, RealDigits and related base conversion functions use recursive divide-and-conquer algorithms. Similar algorithms are used for number input and output.

FilledSmallCircleN uses an adaptive procedure to increase its internal working precision in order to achieve whatever overall precision is requested.

FilledSmallCircleFloor, Ceiling and related functions use an adaptive procedure similar to N to generate exact results from exact input.

Basic arithmetic

FilledSmallCircle Multiplication of large integers and high-precision approximate numbers is done using interleaved schoolbook, Karatsuba, three-way Toom-Cook and number-theoretic transform algorithms.

FilledSmallCircle Machine-code optimization for specific architectures is achieved by using GMP.

FilledSmallCircle Integer powers are found by a left-right binary decomposition algorithm.

FilledSmallCircle Reciprocals and rational powers of approximate numbers use Newton's method.

FilledSmallCircle Exact roots start from numerical estimates.

FilledSmallCircle Significance arithmetic is used for all arithmetic with approximate numbers beyond machine precision.

Pseudorandom numbers

FilledSmallCircleRandom uses the Wolfram rule 30 cellular automaton generator for integers.

FilledSmallCircle It uses a Marsaglia-Zaman subtract-with-borrow generator for real numbers.

Number-theoretical functions

FilledSmallCircleGCD interleaves the HGCD algorithm, the Jebelean-Sorenson-Weber accelerated GCD algorithm, and a combination of Euclid's algorithm and an algorithm based on iterative removal of powers of 2.

FilledSmallCirclePrimeQ first tests for divisibility using small primes, then uses the Miller-Rabin strong pseudoprime test base 2 and base 3, and then uses a Lucas test.

FilledSmallCircle As of 1997, this procedure is known to be correct only for , and it is conceivable that for larger it could claim a composite number to be prime.

FilledSmallCircle The package NumberTheory`PrimeQ` contains a much slower algorithm which has been proved correct for all . It can return an explicit certificate of primality.

FilledSmallCircleFactorInteger switches between trial division, Pollard , Pollard rho and quadratic sieve algorithms.

FilledSmallCircle The package NumberTheory`FactorIntegerECM` contains an elliptic curve algorithm suitable for factoring some very large integers.

FilledSmallCirclePrime and PrimePi use sparse caching and sieving. For large , the Lagarias-Miller-Odlyzko algorithm for PrimePi is used, based on asymptotic estimates of the density of primes, and is inverted to give Prime.

FilledSmallCircleLatticeReduce uses the Lenstra-Lenstra-Lovasz lattice reduction algorithm.

FilledSmallCircle To find a requested number of terms ContinuedFraction uses a modification of Lehmer's indirect method, with a self-restarting divide-and-conquer algorithm to reduce the numerical precision required at each step.

FilledSmallCircleContinuedFraction uses recurrence relations to find periodic continued fractions for quadratic irrationals.

FilledSmallCircleFromContinuedFraction uses iterated matrix multiplication optimized by a divide-and-conquer method.

Combinatorial functions

FilledSmallCircle Most combinatorial functions use sparse caching and recursion.

FilledSmallCircleFactorial, Binomial and related functions use a divide-and-conquer algorithm to balance the number of digits in subproducts.

FilledSmallCircleFibonacci[n] uses an iterative method based on the binary digit sequence of n.

FilledSmallCirclePartitionsP[n] uses Euler's pentagonal formula for small n, and the non-recursive Hardy-Ramanujan-Rademacher method for larger n.

FilledSmallCircleClebschGordan and related functions use generalized hypergeometric series.

Elementary transcendental functions

FilledSmallCircle Exponential and trigonometric functions use Taylor series, stable recursion by argument doubling, and functional relations.

FilledSmallCircle Log and inverse trigonometric functions use Taylor series and functional relations.

Mathematical constants

FilledSmallCircle Values of constants are cached once computed.

FilledSmallCircle Binary splitting is used to subdivide computations of constants.

FilledSmallCirclePi is computed using the Chudnovsky formula.

FilledSmallCircleE is computed from its series expansion.

FilledSmallCircleEulerGamma uses the Brent-McMillan algorithm.

FilledSmallCircleCatalan is computed from a linearly convergent Ramanujan sum.

Special functions

FilledSmallCircle For machine precision most special functions use Mathematica-derived rational minimax approximations. The notes that follow apply mainly to arbitrary precision.

FilledSmallCircle Orthogonal polynomials use stable recursion formulas for polynomial cases and hypergeometric functions in general.

FilledSmallCircleGamma uses recursion, functional equations and the Binet asymptotic formula.

FilledSmallCircle Incomplete gamma and beta functions use hypergeometric series and continued fractions.

FilledSmallCirclePolyGamma uses Euler-Maclaurin summation, functional equations and recursion.

FilledSmallCirclePolyLog uses Euler-Maclaurin summation, expansions in terms of incomplete gamma functions and numerical quadrature.

FilledSmallCircleZeta and related functions use Euler-Maclaurin summation and functional equations. Near the critical strip they also use the Riemann-Siegel formula.

FilledSmallCircleStieltjesGamma uses Keiper's algorithm based on numerical quadrature of an integral representation of the zeta function.

FilledSmallCircle The error function and functions related to exponential integrals are all evaluated using incomplete gamma functions.

FilledSmallCircle The inverse error functions use binomial search and a high-order generalized Newton's method.

FilledSmallCircle Bessel functions use series and asymptotic expansions. For integer orders, some also use stable forward recursion.

FilledSmallCircle The hypergeometric functions use functional equations, stable recurrence relations, series expansions and asymptotic series. Methods from NSum and NIntegrate are also sometimes used.

FilledSmallCircleProductLog uses high-order Newton's method starting from rational approximations and asymptotic expansions.

FilledSmallCircle Elliptic integrals are evaluated using the descending Gauss transformation.

FilledSmallCircle Elliptic theta functions use series summation with recursive evaluation of series terms.

FilledSmallCircle Other elliptic functions mostly use arithmetic-geometric mean methods.

FilledSmallCircle Mathieu functions use Fourier series. The Mathieu characteristic functions use generalizations of Blanch's Newton method.

Numerical integration

FilledSmallCircle With Method->Automatic, NIntegrate uses GaussKronrod in one dimension, and MultiDimensional otherwise.

FilledSmallCircle If an explicit setting for MaxPoints is given, NIntegrate by default uses Method->QuasiMonteCarlo.

FilledSmallCircleGaussKronrod: adaptive Gaussian quadrature with error estimation based on evaluation at Kronrod points.

FilledSmallCircleDoubleExponential: non-adaptive double-exponential quadrature.

FilledSmallCircleTrapezoidal: elementary trapezoidal method.

FilledSmallCircleOscillatory: transformation to handle certain integrals containing trigonometric and Bessel functions.

FilledSmallCircleMultiDimensional: adaptive Genz-Malik algorithm.

FilledSmallCircleMonteCarlo: non-adaptive Monte Carlo.

FilledSmallCircleQuasiMonteCarlo: non-adaptive Halton-Hammersley-Wozniakowski algorithm.

Numerical sums and products

FilledSmallCircle If the ratio test does not give 1, the Wynn epsilon algorithm is applied to a sequence of partial sums or products.

FilledSmallCircle Otherwise Euler-Maclaurin summation is used with Integrate or NIntegrate.

Numerical differential equations

FilledSmallCircle For ordinary differential equations, NDSolve by default uses an LSODA approach, switching between a non-stiff Adams method and a stiff Gear backward differentiation formula method.

FilledSmallCircle For linear boundary value problems the Gel'fand-Lokutsiyevskii chasing method is used.

FilledSmallCircle Differential-algebraic equations use IDA, based on repeated BDF and Newton iteration methods.

FilledSmallCircle For -dimensional PDEs the method of lines is used.

FilledSmallCircleNDSolve supports explicit Method settings that cover most known methods from the literature.

FilledSmallCircle The code for NDSolve and related functions is about 1400 pages long.

Approximate equation solving and optimization

FilledSmallCircle Polynomial root finding is done based on the Jenkins-Traub algorithm.

FilledSmallCircle For sparse linear systems, Solve and NSolve use several efficient numerical methods, mostly based on Gauss factoring with Markowitz products (approximately 250 pages of code).

FilledSmallCircle For systems of algebraic equations, NSolve computes a numerical Gröbner basis using an efficient monomial ordering, then uses eigensystem methods to extract numerical roots.

FilledSmallCircleFindRoot uses a damped Newton's method, the secant method and Brent's method.

FilledSmallCircle With Method->Automatic and two starting values, FindMinimum uses Brent's principal axis method. With one starting value for each variable, FindMinimum uses BFGS quasi-Newton methods, with a limited memory variant for large systems.

FilledSmallCircle If the function to be minimized is a sum of squares, FindMinimum uses the Levenberg-Marquardt method (Method->LevenbergMarquardt).

FilledSmallCircleLinearProgramming uses simplex and revised simplex methods, and with Method->"InteriorPoint" uses primal-dual interior point methods.

FilledSmallCircle For linear cases, NMinimize and NMaximize use the same methods as LinearProgramming. For nonlinear cases, they use Nelder-Mead methods, supplemented by differential evolution, especially when integer variables are present.

Data manipulation

FilledSmallCircleFourier uses the FFT algorithm with decomposition of the length into prime factors. When the prime factors are large, fast convolution methods are used to maintain asymptotic complexity.

FilledSmallCircle For real input, Fourier uses a real transform method.

FilledSmallCircleListConvolve and ListCorrelate use FFT algorithms when possible. For exact integer inputs, enough digits are computed to deduce exact integer results.

FilledSmallCircleInterpolatingFunction uses divided differences to construct Lagrange or Hermite interpolating polynomials.

FilledSmallCircleFit works using singular value decomposition. FindFit uses the same method for the linear least-squares case, the Levenberg-Marquardt method for nonlinear least-squares, and general FindMinimum methods for other norms.

FilledSmallCircleCellularAutomaton uses bit-packed parallel operations with bit slicing. For elementary rules, absolutely optimal Boolean functions are used, while for totalistic rules, just-in-time-compiled bit-packed tables are used. In two dimensions, sparse bit-packed arrays are used when possible, with only active clusters updated.

Approximate numerical linear algebra

FilledSmallCircle Machine-precision matrices are typically converted to a special internal representation for processing.

FilledSmallCircleSparseArray with rules involving patterns uses cylindrical algebraic decomposition to find connected array components. Sparse arrays are stored internally using compressed sparse row formats, generalized for tensors of arbitrary rank.

FilledSmallCircle For dense arrays, LAPACK algorithms extended for arbitrary precision are used when appropriate.

FilledSmallCircle BLAS technology is used to optimize for particular machine architectures.

FilledSmallCircleLUDecomposition, Inverse, RowReduce and Det use Gaussian elimination with partial pivoting. LinearSolve uses the same methods, together with iterative improvement for high-precision numbers.

FilledSmallCircle For sparse arrays, LinearSolve uses UMFPACK multifrontal direct solver methods and with Method->"Krylov" uses Krylov iterative methods preconditioned by an incomplete LU factorization.Eigenvalues and Eigenvectors use ARPACK Arnoldi methods.

FilledSmallCircleSingularValueDecomposition uses the QR algorithm with Givens rotations. PseudoInverse, NullSpace and MatrixRank are based on SingularValueDecomposition.

FilledSmallCircleQRDecomposition uses Householder transformations.

FilledSmallCircleSchurDecomposition uses QR iteration.

FilledSmallCircleMatrixExp uses Schur decomposition.

Exact numerical linear algebra

FilledSmallCircleInverse and LinearSolve use efficient row reduction based on numerical approximation.

FilledSmallCircle With Modulus->n, modular Gaussian elimination is used.

FilledSmallCircleDet uses modular methods and row reduction, constructing a result using the Chinese Remainder Theorem.

FilledSmallCircleEigenvalues works by interpolating the characteristic polynomial.

FilledSmallCircleMatrixExp uses Putzer's method or Jordan decomposition.