Arithmetic
**********

Unless otherwise noted, all of the functions described in this chapter will work for real and complex scalar, vector, or matrix arguments.  Functions described as mapping functions apply the given operation individually to each element when given a matrix argument.  For example:

     sin ([1, 2; 3, 4])
          =>  0.84147   0.90930
              0.14112  -0.75680

Exponents and Logarithms
========================

 -- exp: (X)
     Compute e^x for each element of X.

     To compute the matrix exponential, see Linear Algebra. log

 -- expm1: (X)
     Compute exp (X) - 1 accurately in the neighborhood of zero. exp

 -- log: (X)
     Compute the natural logarithm, ln (X), for each element of X.

     To compute the matrix logarithm, see Linear Algebra. exp, log1p, log2, log10, logspace

 -- reallog: (X)
     Return the real-valued natural logarithm of each element of X.

     If any element results in a complex return value reallog aborts and issues an error. log, realpow, realsqrt

 -- log1p: (X)
     Compute log (1 + X) accurately in the neighborhood of zero. log, exp, expm1

 -- log10: (X)
     Compute the base-10 logarithm of each element of X. log, log2, logspace, exp

 -- Y =: log2 (X)
 -- [F, E] =: log2 (X)
     Compute the base-2 logarithm of each element of X.

     If called with one output, compute the base-2 logarithm such that 2^Y = X.

     If called with two output arguments, split X into binary mantissa (F) and exponent (E) such that X = F * 2^E where 1/2 <= abs (F) < 1 and E is an integer.  If x = 0, f = e = 0. pow2, log, log10, exp

 -- pow2: (X)
 -- pow2: (F, E)
     With one input argument, compute 2 .^ x for each element of X.

     With two input arguments, return f .* (2 .^ e). log2, nextpow2, power

 -- N =: nextpow2 (X)
     Compute the exponent for the smallest power of two larger than the input.

     For each element in the input array X, return the first integer N such that 2^n ≥ abs (x).

     pow2, log2

 -- realpow: (X, Y)
     Compute the real-valued, element-by-element power operator.

     This is equivalent to X .^ Y, except that realpow reports an error if any return value is complex. power, reallog, realsqrt

 -- sqrt: (X)
     Compute the square root of each element of X.

     If X is negative, a complex result is returned.

     To compute the matrix square root, see Linear Algebra. realsqrt, nthroot

 -- realsqrt: (X)
     Return the real-valued square root of each element of X.

     If any element results in a complex return value realsqrt aborts and issues an error. sqrt, realpow, reallog

 -- cbrt: (X)
     Compute the real-valued cube root of each element of X.

     Unlike X^(1/3), the result will be negative if X is negative.

     If any element of X is complex, cbrt aborts with an error. nthroot

 -- nthroot: (X, N)

     Compute the real (non-complex) N-th root of X.

     X must have all real entries and N must be a scalar. If N is an even integer and X has negative entries then nthroot aborts and issues an error.

     Example:

          nthroot (-1, 3)
          => -1
          (-1) ^ (1 / 3)
          => 0.50000 - 0.86603i

     realsqrt, sqrt, cbrt

Complex Arithmetic
==================

In the descriptions of the following functions, Z is the complex number X + IY, where I is defined as sqrt (-1).

 -- abs: (Z)
     Compute the magnitude of Z.

     The magnitude is defined as |Z| = sqrt (x^2 + y^2).

     For example:

          abs (3 + 4i)
               => 5

     arg

 -- arg: (Z)
 -- angle: (Z)
     Compute the argument, i.e., angle of Z.

     This is defined as, THETA = atan2 (Y, X), in radians.

     For example:

          arg (3 + 4i)
               => 0.92730

     abs

 -- conj: (Z)
     Return the complex conjugate of Z.

     The complex conjugate is defined as conj (Z) = X - IY. real, imag

 -- cplxpair: (Z)
 -- cplxpair: (Z, TOL)
 -- cplxpair: (Z, TOL, DIM)
     Sort the numbers Z into complex conjugate pairs ordered by increasing real part.

     The negative imaginary complex numbers are placed first within each pair. All real numbers (those with abs (imag (Z)) / abs (Z) < TOL) are placed after the complex pairs.

     TOL is a weighting factor in the range [0, 1) which determines the tolerance of the matching.  The default value is 100 * eps and the resulting tolerance for a given complex pair is TOL * abs (Z(i))).

     By default the complex pairs are sorted along the first non-singleton dimension of Z.  If DIM is specified, then the complex pairs are sorted along this dimension.

     Signal an error if some complex numbers could not be paired.  Signal an error if all complex numbers are not exact conjugates (to within TOL). Note that there is no defined order for pairs with identical real parts but differing imaginary parts.

          cplxpair (exp (2i*pi*[0:4]'/5)) == exp (2i*pi*[3; 2; 4; 1; 0]/5)

 -- imag: (Z)
     Return the imaginary part of Z as a real number. real, conj

 -- real: (Z)
     Return the real part of Z. imag, conj

Trigonometry
============

Octave provides the following trigonometric functions where angles are specified in radians.  To convert from degrees to radians multiply by pi/180 or use the deg2rad function.  For example, sin (30 * pi/180) returns the sine of 30 degrees.  As an alternative, Octave provides a number of trigonometric functions which work directly on an argument specified in degrees.  These functions are named after the base trigonometric function with a d suffix.  As an example, sin expects an angle in radians while sind expects an angle in degrees.

Octave uses the C library trigonometric functions.  It is expected that these functions are defined by the ISO/IEC 9899 Standard.  This Standard is available at: http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1124.pdf. Section F.9.1 deals with the trigonometric functions.  The behavior of most of the functions is relatively straightforward.  However, there are some exceptions to the standard behavior.  Many of the exceptions involve the behavior for -0.  The most complex case is atan2.  Octave exactly implements the behavior given in the Standard.  Including atan2(+- 0, 0) returns +- pi.

It should be noted that MATLAB uses different definitions which apparently do not distinguish -0.

 -- RAD =: deg2rad (DEG)

     Convert degrees to radians.

     The input DEG must be a scalar, vector, or N-dimensional array of double or single floating point values.  DEG may be complex in which case the real and imaginary components are converted separately.

     The output RAD is the same size and shape as DEG with degrees converted to radians using the conversion constant pi/180.

     Example:

          deg2rad ([0, 90, 180, 270, 360])
            =>  0.00000   1.57080   3.14159   4.71239   6.28319

     rad2deg

 -- DEG =: rad2deg (RAD)

     Convert radians to degrees.

     The input RAD must be a scalar, vector, or N-dimensional array of double or single floating point values.  RAD may be complex in which case the real and imaginary components are converted separately.

     The output DEG is the same size and shape as RAD with radians converted to degrees using the conversion constant 180/pi.

     Example:

          rad2deg ([0, pi/2, pi, 3/2*pi, 2*pi])
            =>  0    90   180   270   360

     deg2rad

 -- sin: (X)
     Compute the sine for each element of X in radians. asin, sind, sinh

 -- cos: (X)
     Compute the cosine for each element of X in radians. acos, cosd, cosh

 -- tan: (Z)
     Compute the tangent for each element of X in radians. atan, tand, tanh

 -- sec: (X)
     Compute the secant for each element of X in radians. asec, secd, sech

 -- csc: (X)
     Compute the cosecant for each element of X in radians. acsc, cscd, csch

 -- cot: (X)
     Compute the cotangent for each element of X in radians. acot, cotd, coth

 -- asin: (X)
     Compute the inverse sine in radians for each element of X. sin, asind

 -- acos: (X)
     Compute the inverse cosine in radians for each element of X. cos, acosd

 -- atan: (X)
     Compute the inverse tangent in radians for each element of X. tan, atand

 -- asec: (X)
     Compute the inverse secant in radians for each element of X. sec, asecd

 -- acsc: (X)
     Compute the inverse cosecant in radians for each element of X. csc, acscd

 -- acot: (X)
     Compute the inverse cotangent in radians for each element of X. cot, acotd

 -- sinh: (X)
     Compute the hyperbolic sine for each element of X. asinh, cosh, tanh

 -- cosh: (X)
     Compute the hyperbolic cosine for each element of X. acosh, sinh, tanh

 -- tanh: (X)
     Compute hyperbolic tangent for each element of X. atanh, sinh, cosh

 -- sech: (X)
     Compute the hyperbolic secant of each element of X. asech

 -- csch: (X)
     Compute the hyperbolic cosecant of each element of X. acsch

 -- coth: (X)
     Compute the hyperbolic cotangent of each element of X. acoth

 -- asinh: (X)
     Compute the inverse hyperbolic sine for each element of X. sinh

 -- acosh: (X)
     Compute the inverse hyperbolic cosine for each element of X. cosh

 -- atanh: (X)
     Compute the inverse hyperbolic tangent for each element of X. tanh

 -- asech: (X)
     Compute the inverse hyperbolic secant of each element of X. sech

 -- acsch: (X)
     Compute the inverse hyperbolic cosecant of each element of X. csch

 -- acoth: (X)
     Compute the inverse hyperbolic cotangent of each element of X. coth

 -- atan2: (Y, X)
     Compute atan (Y / X) for corresponding elements of Y and X.

     Y and X must match in size and orientation.  The signs of elements of Y and X are used to determine the quadrants of each resulting value.

     This function is equivalent to arg (complex (X, Y)). tan, tand, tanh, atanh

Octave provides the following trigonometric functions where angles are specified in degrees.  These functions produce true zeros at the appropriate intervals rather than the small round-off error that occurs when using radians.  For example:

     cosd (90)
          => 0
     cos (pi/2)
          => 6.1230e-17

 -- sind: (X)
     Compute the sine for each element of X in degrees.

     Returns zero for elements where X/180 is an integer. asind, sin

 -- cosd: (X)
     Compute the cosine for each element of X in degrees.

     Returns zero for elements where (X-90)/180 is an integer. acosd, cos

 -- tand: (X)
     Compute the tangent for each element of X in degrees.

     Returns zero for elements where X/180 is an integer and Inf for elements where (X-90)/180 is an integer. atand, tan

 -- secd: (X)
     Compute the secant for each element of X in degrees. asecd, sec

 -- cscd: (X)
     Compute the cosecant for each element of X in degrees. acscd, csc

 -- cotd: (X)
     Compute the cotangent for each element of X in degrees. acotd, cot

 -- asind: (X)
     Compute the inverse sine in degrees for each element of X. sind, asin

 -- acosd: (X)
     Compute the inverse cosine in degrees for each element of X. cosd, acos

 -- atand: (X)
     Compute the inverse tangent in degrees for each element of X. tand, atan

 -- atan2d: (Y, X)
     Compute atan (Y / X) in degrees for corresponding elements from Y and X. tand, atan2

 -- asecd: (X)
     Compute the inverse secant in degrees for each element of X. secd, asec

 -- acscd: (X)
     Compute the inverse cosecant in degrees for each element of X. cscd, acsc

 -- acotd: (X)
     Compute the inverse cotangent in degrees for each element of X. cotd, acot

Sums and Products
=================

 -- sum: (X)
 -- sum: (X, DIM)
 -- sum: (..., "native")
 -- sum: (..., "double")
 -- sum: (..., "extra")
     Sum of elements along dimension DIM.

     If DIM is omitted, it defaults to the first non-singleton dimension.

     The optional "type" input determines the class of the variable used for calculations.  By default, operations on floating point inputs (double or single) are performed in their native data type, while operations on integer, logical, and character data types are performed using doubles.  If the argument "native" is given, then the operation is performed in the same type as the original argument.

     For example:

          sum ([true, true])
             => 2
          sum ([true, true], "native")
             => true

     If "double" is given the sum is performed in double precision even for single precision inputs.

     For double precision inputs, the "extra" option will use a more accurate algorithm than straightforward summation.  For single precision inputs, "extra" is the same as "double".  For all other data type "extra" has no effect. cumsum, sumsq, prod

 -- prod: (X)
 -- prod: (X, DIM)
 -- prod: (..., "native")
 -- prod: (..., "double")
     Product of elements along dimension DIM.

     If DIM is omitted, it defaults to the first non-singleton dimension.

     The optional "type" input determines the class of the variable used for calculations.  If the argument "native" is given, then the operation is performed in the same type as the original argument, rather than the default double type.

     For example:

          prod ([true, true])
             => 1
          prod ([true, true], "native")
             => true

     On the contrary, if "double" is given, the operation is performed in double precision even for single precision inputs. cumprod, sum

 -- cumsum: (X)
 -- cumsum: (X, DIM)
 -- cumsum: (..., "native")
 -- cumsum: (..., "double")
     Cumulative sum of elements along dimension DIM.

     If DIM is omitted, it defaults to the first non-singleton dimension. For example:

          cumsum ([1, 2; 3, 4; 5, 6])
             =>  1   2
                 4   6
                 9  12

     See sum for an explanation of the optional parameters "native" and "double". sum, cumprod

 -- cumprod: (X)
 -- cumprod: (X, DIM)
     Cumulative product of elements along dimension DIM.

     If DIM is omitted, it defaults to the first non-singleton dimension. For example:

          cumprod ([1, 2; 3, 4; 5, 6])
             =>  1   2
                 3   8
                15  48

     prod, cumsum

 -- sumsq: (X)
 -- sumsq: (X, DIM)
     Sum of squares of elements along dimension DIM.

     If DIM is omitted, it defaults to the first non-singleton dimension.

     This function is conceptually equivalent to computing

          sum (x .* conj (x), dim)

     but it uses less memory and avoids calling conj if X is real. sum, prod

Utility Functions
=================

 -- ceil: (X)
     Return the smallest integer not less than X.

     This is equivalent to rounding towards positive infinity.

     If X is complex, return ceil (real (X)) + ceil (imag (X)) * I.

          ceil ([-2.7, 2.7])
              => -2    3

     floor, round, fix

 -- fix: (X)
     Truncate fractional portion of X and return the integer portion.

     This is equivalent to rounding towards zero.  If X is complex, return fix (real (X)) + fix (imag (X)) * I.

          fix ([-2.7, 2.7])
             => -2    2

     ceil, floor, round

 -- floor: (X)
     Return the largest integer not greater than X.

     This is equivalent to rounding towards negative infinity.  If X is complex, return floor (real (X)) + floor (imag (X)) * I.

          floor ([-2.7, 2.7])
               => -3    2

     ceil, round, fix

 -- round: (X)
     Return the integer nearest to X.

     If X is complex, return round (real (X)) + round (imag (X)) * I.  If there are two nearest integers, return the one further away from zero.

          round ([-2.7, 2.7])
               => -3    3

     ceil, floor, fix, roundb

 -- roundb: (X)
     Return the integer nearest to X.  If there are two nearest integers, return the even one (banker's rounding).

     If X is complex, return roundb (real (X)) + roundb (imag (X)) * I. round

 -- max: (X)
 -- max: (X, [], DIM)
 -- [W, IW] =: max (X)
 -- max: (X, Y)
     Find maximum values in the array X.

     For a vector argument, return the maximum value.  For a matrix argument, return a row vector with the maximum value of each column.  For a multi-dimensional array, max operates along the first non-singleton dimension.

     If the optional third argument DIM is present then operate along this dimension.  In this case the second argument is ignored and should be set to the empty matrix.

     For two inputs (X and Y), return the pairwise maximum according to the rules for Broadcasting.

     Thus,

          max (max (X))

     returns the largest element of the 2-D matrix X, and

          max (2:5, pi)
              =>  3.1416  3.1416  4.0000  5.0000

     compares each element of the range 2:5 with pi, and returns a row vector of the maximum values.

     For complex arguments, the magnitude of the elements are used for comparison.  If the magnitudes are identical, then the results are ordered by phase angle in the range (-pi, pi].  Hence,

          max ([-1 i 1 -i])
              => -1

     because all entries have magnitude 1, but -1 has the largest phase angle with value pi.

     If called with one input and two output arguments, max also returns the first index of the maximum value(s).  Thus,

          [x, ix] = max ([1, 3, 5, 2, 5])
              =>  x = 5
                  ix = 3

     min, cummax, cummin

 -- min: (X)
 -- min: (X, [], DIM)
 -- [W, IW] =: min (X)
 -- min: (X, Y)
     Find minimum values in the array X.

     For a vector argument, return the minimum value.  For a matrix argument, return a row vector with the minimum value of each column.  For a multi-dimensional array, min operates along the first non-singleton dimension.

     If the optional third argument DIM is present then operate along this dimension.  In this case the second argument is ignored and should be set to the empty matrix.

     For two inputs (X and Y), return the pairwise minimum according to the rules for Broadcasting.

     Thus,

          min (min (X))

     returns the smallest element of the 2-D matrix X, and

          min (2:5, pi)
              =>  2.0000  3.0000  3.1416  3.1416

     compares each element of the range 2:5 with pi, and returns a row vector of the minimum values.

     For complex arguments, the magnitude of the elements are used for comparison.  If the magnitudes are identical, then the results are ordered by phase angle in the range (-pi, pi].  Hence,

          min ([-1 i 1 -i])
              => -i

     because all entries have magnitude 1, but -i has the smallest phase angle with value -pi/2.

     If called with one input and two output arguments, min also returns the first index of the minimum value(s).  Thus,

          [x, ix] = min ([1, 3, 0, 2, 0])
              =>  x = 0
                  ix = 3

     max, cummin, cummax

 -- cummax: (X)
 -- cummax: (X, DIM)
 -- [W, IW] =: cummax (...)
     Return the cumulative maximum values along dimension DIM.

     If DIM is unspecified it defaults to column-wise operation.  For example:

          cummax ([1 3 2 6 4 5])
             =>  1  3  3  6  6  6

     If called with two output arguments the index of the maximum value is also returned.

          [w, iw] = cummax ([1 3 2 6 4 5])
          =>
          w =  1  3  3  6  6  6
          iw = 1  2  2  4  4  4

     cummin, max, min

 -- cummin: (X)
 -- cummin: (X, DIM)
 -- [W, IW] =: cummin (X)
     Return the cumulative minimum values along dimension DIM.

     If DIM is unspecified it defaults to column-wise operation.  For example:

          cummin ([5 4 6 2 3 1])
             =>  5  4  4  2  2  1

     If called with two output arguments the index of the minimum value is also returned.

          [w, iw] = cummin ([5 4 6 2 3 1])
          =>
          w =  5  4  4  2  2  1
          iw = 1  2  2  4  4  6

     cummax, min, max

 -- hypot: (X, Y)
 -- hypot: (X, Y, Z, ...)
     Compute the element-by-element square root of the sum of the squares of X and Y.

     This is equivalent to sqrt (X.^2 + Y.^2), but is calculated in a manner that avoids overflows for large values of X or Y.

     hypot can also be called with more than 2 arguments; in this case, the arguments are accumulated from left to right:

          hypot (hypot (X, Y), Z)
          hypot (hypot (hypot (X, Y), Z), W), etc.

 -- DX =: gradient (M)
 -- [DX, DY, DZ, ...] =: gradient (M)
 -- [...] =: gradient (M, S)
 -- [...] =: gradient (M, X, Y, Z, ...)
 -- [...] =: gradient (F, X0)
 -- [...] =: gradient (F, X0, S)
 -- [...] =: gradient (F, X0, X, Y, ...)

     Calculate the gradient of sampled data or a function.

     If M is a vector, calculate the one-dimensional gradient of M. If M is a matrix the gradient is calculated for each dimension.

     [DX, DY] = gradient (M) calculates the one-dimensional gradient for X and Y direction if M is a matrix.  Additional return arguments can be use for multi-dimensional matrices.

     A constant spacing between two points can be provided by the S parameter.  If S is a scalar, it is assumed to be the spacing for all dimensions.  Otherwise, separate values of the spacing can be supplied by the X, ... arguments.  Scalar values specify an equidistant spacing.  Vector values for the X, ... arguments specify the coordinate for that dimension.  The length must match their respective dimension of M.

     At boundary points a linear extrapolation is applied.  Interior points are calculated with the first approximation of the numerical gradient

          y'(i) = 1/(x(i+1)-x(i-1)) * (y(i-1)-y(i+1)).

     If the first argument F is a function handle, the gradient of the function at the points in X0 is approximated using central difference. For example, gradient (@cos, 0) approximates the gradient of the cosine function in the point x0 = 0.  As with sampled data, the spacing values between the points from which the gradient is estimated can be set via the S or DX, DY, ... arguments.  By default a spacing of 1 is used. diff, del2

 -- dot: (X, Y, DIM)
     Compute the dot product of two vectors.

     If X and Y are matrices, calculate the dot products along the first non-singleton dimension.

     If the optional argument DIM is given, calculate the dot products along this dimension.

     Implementation Note: This is equivalent to sum (conj (X) .* Y, DIM), but avoids forming a temporary array and is faster.  When X and Y are column vectors, the result is equivalent to X' * Y.  Although, dot is defined for integer arrays, the output may differ from the expected result due to the limited range of integer objects. cross, divergence

 -- cross: (X, Y)
 -- cross: (X, Y, DIM)
     Compute the vector cross product of two 3-dimensional vectors X and Y.

     If X and Y are matrices, the cross product is applied along the first dimension with three elements.

     The optional argument  DIM forces the cross product to be calculated along the specified dimension.

     Example Code:

          cross ([1, 1, 0], [0, 1, 1])
            =>
                 1   -1   1

     dot, curl, divergence

 -- DIV =: divergence (X, Y, Z, FX, FY, FZ)
 -- DIV =: divergence (FX, FY, FZ)
 -- DIV =: divergence (X, Y, FX, FY)
 -- DIV =: divergence (FX, FY)
     Calculate divergence of a vector field given by the arrays FX, FY, and FZ or FX, FY respectively.

                            d               d               d
          div F(x,y,z)  =   -- F(x,y,z)  +  -- F(x,y,z)  +  -- F(x,y,z)
                            dx              dy              dz

     The coordinates of the vector field can be given by the arguments X, Y, Z or X, Y respectively.

     curl, gradient, del2, dot

 -- [CX, CY, CZ, V] =: curl (X, Y, Z, FX, FY, FZ)
 -- [CZ, V] =: curl (X, Y, FX, FY)
 -- [...] =: curl (FX, FY, FZ)
 -- [...] =: curl (FX, FY)
 -- V =: curl (...)
     Calculate curl of vector field given by the arrays FX, FY, and FZ or FX, FY respectively.

                            / d         d       d         d       d         d     \
          curl F(x,y,z)  =  | -- Fz  -  -- Fy,  -- Fx  -  -- Fz,  -- Fy  -  -- Fx |
                            \ dy        dz      dz        dx      dx        dy    /

     The coordinates of the vector field can be given by the arguments X, Y, Z or X, Y respectively.  V calculates the scalar component of the angular velocity vector in direction of the z-axis for two-dimensional input.  For three-dimensional input the scalar rotation is calculated at each grid point in direction of the vector field at that point. divergence, gradient, del2, cross

 -- L =: del2 (M)
 -- L =: del2 (M, H)
 -- L =: del2 (M, DX, DY, ...)

     Calculate the discrete Laplace operator.

     For a 2-dimensional matrix M this is defined as

                1    / d^2            d^2         \
          L  = --- * | ---  M(x,y) +  ---  M(x,y) |
                4    \ dx^2           dy^2        /

     For N-dimensional arrays the sum in parentheses is expanded to include second derivatives over the additional higher dimensions.

     The spacing between evaluation points may be defined by H, which is a scalar defining the equidistant spacing in all dimensions.  Alternatively, the spacing in each dimension may be defined separately by DX, DY, etc.  A scalar spacing argument defines equidistant spacing, whereas a vector argument can be used to specify variable spacing.  The length of the spacing vectors must match the respective dimension of M.  The default spacing value is 1.

     Dimensions with fewer than 3 data points are skipped.  Boundary points are calculated from the linear extrapolation of interior points.

     Example: Second derivative of 2*x^3

          f = @(x) 2*x.^3;
          dd = @(x) 12*x;
          x = 1:6;
          L = 4*del2 (f(x));
          assert (L, dd (x));

     gradient, diff

 -- factorial: (N)
     Return the factorial of N where N is a real non-negative integer.

     If N is a scalar, this is equivalent to prod (1:N).  For vector or matrix arguments, return the factorial of each element in the array.

     For non-integers see the generalized factorial function gamma. Note that the factorial function grows large quite quickly, and even with double precision values overflow will occur if N > 171.  For such cases consider gammaln. prod, gamma, gammaln

 -- PF =: factor (Q)
 -- [PF, N] =: factor (Q)
     Return the prime factorization of Q.

     The prime factorization is defined as prod (PF) == Q where every element of PF is a prime number.  If Q == 1, return 1.  The output PF is of the same numeric class as the input.

     With two output arguments, return the unique prime factors PF and their multiplicities.  That is, prod (PF .^ N) == Q.

     Implementation Note: The input Q must be less than flintmax (9.0072e+15) in order to factor correctly. gcd, lcm, isprime, primes

 -- G =: gcd (A1, A2, ...)
 -- [G, V1, ...] =: gcd (A1, A2, ...)
     Compute the greatest common divisor of A1, A2, ....

     All arguments must be the same size or scalar.  For arrays, the greatest common divisor is calculated for each element individually.  All elements must be ordinary or Gaussian (complex) integers.  Note that for Gaussian integers, the gcd is only unique up to a phase factor (multiplication by 1, -1, i, or -i), so an arbitrary greatest common divisor among the four possible is returned.

     Optional return arguments V1, ..., contain integer vectors such that,

          G = V1 .* A1 + V2 .* A2 + ...

     Example code:

          gcd ([15, 9], [20, 18])
             =>  5  9

     lcm, factor, isprime

 -- lcm: (X, Y)
 -- lcm: (X, Y, ...)
     Compute the least common multiple of X and Y, or of the list of all arguments.

     All elements must be numeric and of the same size or scalar. factor, gcd, isprime

 -- rem: (X, Y)
     Return the remainder of the division X / Y.

     The remainder is computed using the expression

          x - y .* fix (x ./ y)

     An error message is printed if the dimensions of the arguments do not agree, or if either argument is complex.

     Programming Notes: When calculating with floating point numbers (double, single), values within a few eps of an integer will be rounded to that integer before computation for compatibility with MATLAB.  Any floating point integers greater than flintmax (2^53 for double) will not compute correctly.  For larger integer values convert the input to uint64 before calling this function.

     By convention,

          rem (X, 0) = NaN  if X is a floating point variable
          rem (X, 0) = 0    if X is an integer variable
          rem (X, Y)  returns a value with the signbit from X

     For the opposite conventions see the mod function.  In general, rem is best when computing the remainder after division of two positive numbers.  For negative numbers, or when the values are periodic, mod is a better choice. mod

 -- mod: (X, Y)
     Compute the modulo of X and Y.

     Conceptually this is given by

          x - y .* floor (x ./ y)

     and is written such that the correct modulus is returned for integer types. This function handles negative values correctly.  That is, mod (-1, 3) is 2, not -1, as rem (-1, 3) returns.

     An error results if the dimensions of the arguments do not agree, or if either of the arguments is complex.

     Programming Notes: When calculating with floating point numbers (double, single), values within a few eps of an integer will be rounded to that integer before computation for compatibility with MATLAB.  Any floating point integers greater than flintmax (2^53 for double) will not compute correctly.  For larger integer values convert the input to uint64 before calling this function.

     By convention,

          mod (X, 0) = X
          mod (X, Y)      returns a value with the signbit from Y

     For the opposite conventions see the rem function.  In general, mod is a better choice than rem when any of the inputs are negative numbers or when the values are periodic. rem

 -- P =: primes (N)
     Return all primes up to N.

     The output data class (double, single, uint32, etc.) is the same as the input class of N.  The algorithm used is the Sieve of Eratosthenes.

     Note: If you need a specific number of primes you can use the fact that the distance from one prime to the next is, on average, proportional to the logarithm of the prime.  Integrating, one finds that there are about k primes less than k*log (5*k).

     See also list_primes if you need a specific number N of primes. list_primes, isprime

 -- list_primes: ()
 -- list_primes: (N)
     List the first N primes.

     If N is unspecified, the first 25 primes are listed. primes, isprime

 -- sign: (X)
     Compute the signum function.

     This is defined as

                     -1, x < 0;
          sign (x) =  0, x = 0;
                      1, x > 0.

     For complex arguments, sign returns x ./ abs (X).

     Note that sign (-0.0) is 0.  Although IEEE 754 floating point allows zero to be signed, 0.0 and -0.0 compare equal.  If you must test whether zero is signed, use the signbit function. signbit

 -- signbit: (X)
     Return logical true if the value of X has its sign bit set and false otherwise.

     This behavior is consistent with the other logical functions. See Logical Values.  The behavior differs from the C language function which returns nonzero if the sign bit is set.

     This is not the same as x < 0.0, because IEEE 754 floating point allows zero to be signed.  The comparison -0.0 < 0.0 is false, but signbit (-0.0) will return a nonzero value. sign

Special Functions
=================

 -- [A, IERR] =: airy (K, Z, OPT)
     Compute Airy functions of the first and second kind, and their derivatives.

           K   Function   Scale factor (if "opt" is supplied)
          ---  --------   ---------------------------------------
           0   Ai (Z)     exp ((2/3) * Z * sqrt (Z))
           1   dAi(Z)/dZ  exp ((2/3) * Z * sqrt (Z))
           2   Bi (Z)     exp (-abs (real ((2/3) * Z * sqrt (Z))))
           3   dBi(Z)/dZ  exp (-abs (real ((2/3) * Z * sqrt (Z))))

     The function call airy (Z) is equivalent to airy (0, Z).

     The result is the same size as Z.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Loss of significance by argument reduction, output may be inaccurate.

      5. Error---no computation, algorithm termination condition not met, return NaN.

 -- J =: besselj (ALPHA, X)
 -- J =: besselj (ALPHA, X, OPT)
 -- [J, IERR] =: besselj (...)
     Compute Bessel functions of the first kind.

     The order of the Bessel function ALPHA must be real.  The points for evaluation X may be complex.

     If the optional argument OPT is 1 or true, the result J is multiplied by exp (-abs (imag (X))).

     If ALPHA is a scalar, the result is the same size as X.  If X is a scalar, the result is the same size as ALPHA.  If ALPHA is a row vector and X is a column vector, the result is a matrix with length (X) rows and length (ALPHA) columns. Otherwise, ALPHA and X must conform and the result will be the same size.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Loss of significance by argument reduction, output may be inaccurate.

      5. Error---no computation, algorithm termination condition not met, return NaN.

     bessely, besseli, besselk, besselh

 -- Y =: bessely (ALPHA, X)
 -- Y =: bessely (ALPHA, X, OPT)
 -- [Y, IERR] =: bessely (...)
     Compute Bessel functions of the second kind.

     The order of the Bessel function ALPHA must be real.  The points for evaluation X may be complex.

     If the optional argument OPT is 1 or true, the result Y is multiplied by exp (-abs (imag (X))).

     If ALPHA is a scalar, the result is the same size as X.  If X is a scalar, the result is the same size as ALPHA.  If ALPHA is a row vector and X is a column vector, the result is a matrix with length (X) rows and length (ALPHA) columns. Otherwise, ALPHA and X must conform and the result will be the same size.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Complete loss of significance by argument reduction, return NaN.

      5. Error---no computation, algorithm termination condition not met, return NaN.

     besselj, besseli, besselk, besselh

 -- I =: besseli (ALPHA, X)
 -- I =: besseli (ALPHA, X, OPT)
 -- [I, IERR] =: besseli (...)
     Compute modified Bessel functions of the first kind.

     The order of the Bessel function ALPHA must be real.  The points for evaluation X may be complex.

     If the optional argument OPT is 1 or true, the result I is multiplied by exp (-abs (real (X))).

     If ALPHA is a scalar, the result is the same size as X.  If X is a scalar, the result is the same size as ALPHA.  If ALPHA is a row vector and X is a column vector, the result is a matrix with length (X) rows and length (ALPHA) columns. Otherwise, ALPHA and X must conform and the result will be the same size.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Complete loss of significance by argument reduction, return NaN.

      5. Error---no computation, algorithm termination condition not met, return NaN.

     besselk, besselj, bessely, besselh

 -- K =: besselk (ALPHA, X)
 -- K =: besselk (ALPHA, X, OPT)
 -- [K, IERR] =: besselk (...)

     Compute modified Bessel functions of the second kind.

     The order of the Bessel function ALPHA must be real.  The points for evaluation X may be complex.

     If the optional argument OPT is 1 or true, the result K is multiplied by exp (X).

     If ALPHA is a scalar, the result is the same size as X.  If X is a scalar, the result is the same size as ALPHA.  If ALPHA is a row vector and X is a column vector, the result is a matrix with length (X) rows and length (ALPHA) columns. Otherwise, ALPHA and X must conform and the result will be the same size.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Complete loss of significance by argument reduction, return NaN.

      5. Error---no computation, algorithm termination condition not met, return NaN.

     besseli, besselj, bessely, besselh

 -- H =: besselh (ALPHA, X)
 -- H =: besselh (ALPHA, K, X)
 -- H =: besselh (ALPHA, K, X, OPT)
 -- [H, IERR] =: besselh (...)
     Compute Bessel functions of the third kind (Hankel functions).

     The order of the Bessel function ALPHA must be real.  The kind of Hankel function is specified by K and may be either first (K = 1) or second (K = 2).  The default is Hankel functions of the first kind.  The points for evaluation X may be complex.

     If the optional argument OPT is 1 or true, the result is multiplied by exp (-I*X) for K = 1 or exp (I*X) for K = 2.

     If ALPHA is a scalar, the result is the same size as X.  If X is a scalar, the result is the same size as ALPHA.  If ALPHA is a row vector and X is a column vector, the result is a matrix with length (X) rows and length (ALPHA) columns. Otherwise, ALPHA and X must conform and the result will be the same size.

     If requested, IERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Input error, return NaN.

      2. Overflow, return Inf.

      3. Loss of significance by argument reduction results in less than half of machine accuracy.

      4. Complete loss of significance by argument reduction, return NaN.

      5. Error---no computation, algorithm termination condition not met, return NaN.

     besselj, bessely, besseli, besselk

 -- beta: (A, B)
     Compute the Beta function for real inputs A and B.

     The Beta function definition is

          beta (a, b) = gamma (a) * gamma (b) / gamma (a + b).

     The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly. See betaln, for computing the logarithm of the Beta function in an efficient manner. betaln, betainc, betaincinv

 -- betainc: (X, A, B)
 -- betainc: (X, A, B, TAIL)
     Compute the incomplete beta function.

     This is defined as

                                    x
                                   /
                           1       |
          I_x (a, b) = ----------  | t^(a-1) (1-t)^(b-1) dt
                       beta (a,b)  |
                                   /
                                  0

     with real X in the range [0,1].  The inputs A and B must be real and strictly positive (> 0).  If one of the inputs is not a scalar then the other inputs must be scalar or of compatible dimensions.

     By default, TAIL is "lower" and the incomplete beta function integrated from 0 to X is computed.  If TAIL is "upper" then the complementary function integrated from X to 1 is calculated. The two choices are related by

     betainc (X, A, B, "upper") = 1 - betainc (X, A, B, "lower").

     betainc uses a more sophisticated algorithm than subtraction to get numerically accurate results when the "lower" value is small.

     Reference: A. Cuyt, V. Brevik Petersen, B. Verdonk, H. Waadeland, W.B. Jones, Handbook of Continued Fractions for Special Functions, ch. 18.

     beta, betaincinv, betaln

 -- betaincinv: (Y, A, B)
 -- betaincinv: (Y, A, B, "lower")
 -- betaincinv: (Y, A, B, "upper")
     Compute the inverse of the normalized incomplete beta function.

     The normalized incomplete beta function is defined as

                                    x
                                   /
                           1       |
          I_x (a, b) = ----------  | t^(a-1) (1-t)^(b-1) dt
                       beta (a,b)  |
                                   /
                                  0

     If two inputs are scalar, then betaincinv (Y, A, B) is returned for each of the other inputs.

     If two or more inputs are not scalar, the sizes of them must agree, and betaincinv is applied element-by-element.

     The variable Y must be in the interval [0,1], while A and B must be real and strictly positive.

     By default, TAIL is "lower" and the inverse of the incomplete beta function integrated from 0 to X is computed.  If TAIL is "upper" then the complementary function integrated from X to 1 is inverted.

     The function is computed by standard Newton's method, by solving

          Y - betainc (X, A, B) = 0

     betainc, beta, betaln

 -- betaln: (A, B)
     Compute the natural logarithm of the Beta function for real inputs A and B.

     betaln is defined as

          betaln (a, b) = log (beta (a, b))

     and is calculated in a way to reduce the occurrence of underflow.

     The Beta function can grow quite large and it is often more useful to work with the logarithm of the output rather than the function directly. beta, betainc, betaincinv, gammaln

 -- bincoeff: (N, K)
     Return the binomial coefficient of N and K.

     The binomial coefficient is defined as

           /   \
           | n |    n (n-1) (n-2) ... (n-k+1)
           |   |  = -------------------------
           | k |               k!
           \   /

     For example:

          bincoeff (5, 2)
             => 10

     In most cases, the nchoosek function is faster for small scalar integer arguments.  It also warns about loss of precision for big arguments.

     nchoosek

 -- commutation_matrix: (M, N)
     Return the commutation matrix K(m,n) which is the unique M*N by M*N matrix such that K(m,n) * vec(A) = vec(A') for all m by n matrices A.

     If only one argument M is given, K(m,m) is returned.

     See Magnus and Neudecker (1988), Matrix Differential Calculus with Applications in Statistics and Econometrics.

 -- cosint: (X)
     Compute the cosine integral function:

                      +oo
                     /
          Ci (x) = - | (cos (t)) / t dt
                     /
                    x

     An equivalent definition is

                                       x
                                      /
                                      |  cos (t) - 1
          Ci (x) = gamma + log (x) +  | -------------  dt
                                      |        t
                                      /
                                     0

     Reference:

     M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1964.

     sinint, expint, cos

 -- duplication_matrix: (N)
     Return the duplication matrix Dn which is the unique n^2 by n*(n+1)/2 matrix such that Dn vech (A) = vec (A) for all symmetric n by n matrices A.

     See Magnus and Neudecker (1988), Matrix Differential Calculus with Applications in Statistics and Econometrics.

 -- dawson: (Z)
     Compute the Dawson (scaled imaginary error) function.

     The Dawson function is defined as

          (sqrt (pi) / 2) * exp (-z^2) * erfi (z)

     erfc, erf, erfcx, erfi, erfinv, erfcinv

 -- [SN, CN, DN, ERR] =: ellipj (U, M)
 -- [SN, CN, DN, ERR] =: ellipj (U, M, TOL)
     Compute the Jacobi elliptic functions SN, CN, and DN of complex argument U and real parameter M.

     If M is a scalar, the results are the same size as U. If U is a scalar, the results are the same size as M. If U is a column vector and M is a row vector, the results are matrices with length (U) rows and length (M) columns.  Otherwise, U and M must conform in size and the results will be the same size as the inputs.

     The value of U may be complex. The value of M must be 0 ≤ M ≤ 1.

     The optional input TOL is currently ignored (MATLAB uses this to allow faster, less accurate approximation).

     If requested, ERR contains the following status information and is the same size as the result.

      0. Normal return.

      1. Error---no computation, algorithm termination condition not met, return NaN.

     Reference: Milton Abramowitz and Irene A Stegun, Handbook of Mathematical Functions, Chapter 16 (Sections 16.4, 16.13, and 16.15), Dover, 1965.

     ellipke

 -- K =: ellipke (M)
 -- K =: ellipke (M, TOL)
 -- [K, E] =: ellipke (...)
     Compute complete elliptic integrals of the first K(M) and second E(M) kind.

     M must be a scalar or real array with -Inf ≤ M ≤ 1.

     The optional input TOL controls the stopping tolerance of the algorithm and defaults to eps (class (M)).  The tolerance can be increased to compute a faster, less accurate approximation.

     When called with one output only elliptic integrals of the first kind are returned.

     Mathematical Note:

     Elliptic integrals of the first kind are defined as

                   1
                  /               dt
          K (m) = | ------------------------------
                  / sqrt ((1 - t^2)*(1 - m*t^2))
                 0

     Elliptic integrals of the second kind are defined as

                   1
                  /  sqrt (1 - m*t^2)
          E (m) = |  ------------------ dt
                  /  sqrt (1 - t^2)
                 0

     Reference: Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions, Chapter 17, Dover, 1965. ellipj

 -- erf: (Z)
     Compute the error function.

     The error function is defined as

                                  z
                        2        /
          erf (z) = --------- *  | e^(-t^2) dt
                    sqrt (pi)    /
                              t=0

     erfc, erfcx, erfi, dawson, erfinv, erfcinv

 -- erfc: (Z)
     Compute the complementary error function.

     The complementary error function is defined as 1 - erf (Z). erfcinv, erfcx, erfi, dawson, erf, erfinv

 -- erfcx: (Z)
     Compute the scaled complementary error function.

     The scaled complementary error function is defined as

          exp (z^2) * erfc (z)

     erfc, erf, erfi, dawson, erfinv, erfcinv

 -- erfi: (Z)
     Compute the imaginary error function.

     The imaginary error function is defined as

          -i * erf (i*z)

     erfc, erf, erfcx, dawson, erfinv, erfcinv

 -- erfinv: (X)
     Compute the inverse error function.

     The inverse error function is defined such that

          erf (Y) == X

     erf, erfc, erfcx, erfi, dawson, erfcinv

 -- erfcinv: (X)
     Compute the inverse complementary error function.

     The inverse complementary error function is defined such that

          erfc (Y) == X

     erfc, erf, erfcx, erfi, dawson, erfinv

 -- expint: (X)
     Compute the exponential integral.

     The exponential integral is defined as:

                     +oo
                    /
                    | exp (-t)
          E_1 (x) = | -------- dt
                    |    t
                    /
                   x

     Note: For compatibility, this function uses the MATLAB definition of the exponential integral.  Most other sources refer to this particular value as E_1 (x), and the exponential integral as

                      +oo
                     /
                     | exp (-t)
          Ei (x) = - | -------- dt
                     |    t
                     /
                   -x

     The two definitions are related, for positive real values of X, by E_1 (-x) = -Ei (x) - i*pi.

     References:

     M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1964.

     N. Bleistein and R.A. Handelsman, Asymptotic expansions of integrals, 1986.

     cosint, sinint, exp

 -- gamma: (Z)
     Compute the Gamma function.

     The Gamma function is defined as

                       infinity
                      /
          gamma (z) = | t^(z-1) exp (-t) dt.
                      /
                   t=0

     Programming Note: The gamma function can grow quite large even for small input values.  In many cases it may be preferable to use the natural logarithm of the gamma function (gammaln) in calculations to minimize loss of precision.  The final result is then exp (RESULT_USING_GAMMALN). gammainc, gammaln, factorial

 -- gammainc: (X, A)
 -- gammainc: (X, A, TAIL)
     Compute the normalized incomplete gamma function.

     This is defined as

                                          x
                                 1       /
          gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt
                            gamma (a)    /
                                      t=0

     with the limiting value of 1 as X approaches infinity. The standard notation is P(a,x), e.g., Abramowitz and Stegun (6.5.1).

     If A is scalar, then gammainc (X, A) is returned for each element of X and vice versa.

     If neither X nor A is scalar then the sizes of X and A must agree, and gammainc is applied element-by-element. The elements of A must be non-negative.

     By default, TAIL is "lower" and the incomplete gamma function integrated from 0 to X is computed.  If TAIL is "upper" then the complementary function integrated from X to infinity is calculated.

     If TAIL is "scaledlower", then the lower incomplete gamma function is multiplied by gamma(a+1)*exp(x)/(x^a). If TAIL is "scaledupper", then the upper incomplete gamma function is multiplied by the same quantity.

     References:

     M. Abramowitz and I.A. Stegun, Handbook of mathematical functions, Dover publications, Inc., 1972.

     W. Gautschi, A computational procedure for incomplete gamma functions, ACM Trans. Math Software, pp. 466--481, Vol 5, No. 4, 2012.

     W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran 77, ch. 6.2, Vol 1, 1992.

     gamma, gammaincinv, gammaln

 -- gammaincinv: (Y, A)
 -- gammaincinv: (Y, A, TAIL)
     Compute the inverse of the normalized incomplete gamma function.

     The normalized incomplete gamma function is defined as

                                          x
                                 1       /
          gammainc (x, a) = ---------    | exp (-t) t^(a-1) dt
                            gamma (a)    /
                                      t=0

     and gammaincinv (gammainc (X, A), A) = X for each non-negative value of X.  If A is scalar then gammaincinv (Y, A) is returned for each element of Y and vice versa.

     If neither Y nor A is scalar then the sizes of Y and A must agree, and gammaincinv is applied element-by-element. The variable Y must be in the interval [0,1] while A must be real and positive.

     By default, TAIL is "lower" and the inverse of the incomplete gamma function integrated from 0 to X is computed.  If TAIL is "upper", then the complementary function integrated from X to infinity is inverted.

     The function is computed with Newton's method by solving

          Y - gammainc (X, A) = 0

     Reference: A. Gil, J. Segura, and N. M. Temme, Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function ratios, SIAM J. Sci. Computing, pp. A2965--A2981, Vol 34, 2012.

     gammainc, gamma, gammaln

 -- L =: legendre (N, X)
 -- L =: legendre (N, X, NORMALIZATION)
     Compute the associated Legendre function of degree N and order M = 0 ... N.

     The value N must be a real non-negative integer.

     X is a vector with real-valued elements in the range [-1, 1].

     The optional argument NORMALIZATION may be one of "unnorm", "sch", or "norm".  The default if no normalization is given is "unnorm".

     When the optional argument NORMALIZATION is "unnorm", compute the associated Legendre function of degree N and order M and return all values for M = 0 ... N.  The return value has one dimension more than X.

     The associated Legendre function of degree N and order M:

           m         m      2  m/2   d^m
          P(x) = (-1) * (1-x  )    * ----  P(x)
           n                         dx^m   n

     with Legendre polynomial of degree N:

                    1    d^n   2    n
          P(x) = ------ [----(x - 1) ]
           n     2^n n!  dx^n

     legendre (3, [-1.0, -0.9, -0.8]) returns the matrix:

           x  |   -1.0   |   -0.9   |   -0.8
          ------------------------------------
          m=0 | -1.00000 | -0.47250 | -0.08000
          m=1 |  0.00000 | -1.99420 | -1.98000
          m=2 |  0.00000 | -2.56500 | -4.32000
          m=3 |  0.00000 | -1.24229 | -3.24000

     When the optional argument NORMALIZATION is "sch", compute the Schmidt semi-normalized associated Legendre function.  The Schmidt semi-normalized associated Legendre function is related to the unnormalized Legendre functions by the following:

     For Legendre functions of degree N and order 0:

            0      0
          SP(x) = P(x)
            n      n

     For Legendre functions of degree n and order m:

            m      m         m    2(n-m)! 0.5
          SP(x) = P(x) * (-1)  * [-------]
            n      n              (n+m)!

     When the optional argument NORMALIZATION is "norm", compute the fully normalized associated Legendre function.  The fully normalized associated Legendre function is related to the unnormalized associated Legendre functions by the following:

     For Legendre functions of degree N and order M

            m      m         m    (n+0.5)(n-m)! 0.5
          NP(x) = P(x) * (-1)  * [-------------]
            n      n                  (n+m)!

 -- gammaln: (X)
 -- lgamma: (X)
     Return the natural logarithm of the gamma function of X. gamma, gammainc

 -- psi: (Z)
 -- psi: (K, Z)
     Compute the psi (polygamma) function.

     The polygamma functions are the Kth derivative of the logarithm of the gamma function.  If unspecified, K defaults to zero.  A value of zero computes the digamma function, a value of 1, the trigamma function, and so on.

     The digamma function is defined:

          psi (z) = d (log (gamma (z))) / dx

     When computing the digamma function (when K equals zero), Z can have any value real or complex value.  However, for polygamma functions (K higher than 0), Z must be real and non-negative.

     gamma, gammainc, gammaln

 -- sinint: (X)
     Compute the sine integral function:

                     x
                    /
          Si (x) =  | sin (t) / t dt
                    /
                   0

     Reference: M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, 1964.

     cosint, expint, sin

Rational Approximations
=======================

 -- S =: rat (X)
 -- S =: rat (X, TOL)
 -- [N, D] =: rat (...)

     Find a rational approximation of X to within the tolerance defined by TOL.

     If unspecified, the default tolerance is 1e-6 * norm (X(:), 1).

     When called with one output argument, return a string containing a continued fraction expansion (multiple terms).

     When called with two output arguments, return numeric matrices for the numerator and denominator of a fractional representation of X such that X = N ./ D.

     For example:

          s = rat (pi)
          => s = 3 + 1/(7 + 1/16)

          [n, d] = rat (pi)
          => n =  355
          => d =  113

          n / d - pi
          => 0.00000026676

     Programming Note: With one output rat produces a string which is a continued fraction expansion.  To produce a string which is a simple fraction (one numerator, one denominator) use rats.

     rats, format

 -- S =: rats (X)
 -- S =: rats (X, LEN)
     Convert X into a rational approximation represented as a string.

     A rational approximation to a floating point number is a simple fraction with numerator N and denominator D such that X = N/D.

     The optional second argument defines the maximum length of the string representing the elements of X.  By default, LEN is 9.

     If the length of the smallest possible rational approximation exceeds LEN, an asterisk (*) padded with spaces will be returned instead.

     Example conversion from matrix to string, and back again.

          r = rats (hilb (4));
          x = str2num (r)

     rat, format

Coordinate Transformations
==========================

 -- [THETA, R] =: cart2pol (X, Y)
 -- [THETA, R, Z] =: cart2pol (X, Y, Z)
 -- [THETA, R] =: cart2pol (C)
 -- [THETA, R, Z] =: cart2pol (C)
 -- P =: cart2pol (...)

     Transform Cartesian coordinates to polar or cylindrical coordinates.

     The inputs X, Y (, and Z) must be the same shape, or scalar.  If called with a single matrix argument then each row of C represents the Cartesian coordinate (X, Y (, Z)).

     THETA describes the angle relative to the positive x-axis.

     R is the distance to the z-axis (0, 0, z).

     If only a single return argument is requested then return a matrix P where each row represents one polar/(cylindrical) coordinate (THETA, PHI (, Z)). pol2cart, cart2sph, sph2cart

 -- [X, Y] =: pol2cart (THETA, R)
 -- [X, Y, Z] =: pol2cart (THETA, R, Z)
 -- [X, Y] =: pol2cart (P)
 -- [X, Y, Z] =: pol2cart (P)
 -- C =: pol2cart (...)
     Transform polar or cylindrical coordinates to Cartesian coordinates.

     The inputs THETA, R, (and Z) must be the same shape, or scalar.  If called with a single matrix argument then each row of P represents the polar/(cylindrical) coordinate (THETA, R (, Z)).

     THETA describes the angle relative to the positive x-axis.

     R is the distance to the z-axis (0, 0, z).

     If only a single return argument is requested then return a matrix C where each row represents one Cartesian coordinate (X, Y (, Z)). cart2pol, sph2cart, cart2sph

 -- [THETA, PHI, R] =: cart2sph (X, Y, Z)
 -- [THETA, PHI, R] =: cart2sph (C)
 -- S =: cart2sph (...)
     Transform Cartesian coordinates to spherical coordinates.

     The inputs X, Y, and Z must be the same shape, or scalar. If called with a single matrix argument then each row of C represents the Cartesian coordinate (X, Y, Z).

     THETA describes the angle relative to the positive x-axis.

     PHI is the angle relative to the xy-plane.

     R is the distance to the origin (0, 0, 0).

     If only a single return argument is requested then return a matrix S where each row represents one spherical coordinate (THETA, PHI, R). sph2cart, cart2pol, pol2cart

 -- [X, Y, Z] =: sph2cart (THETA, PHI, R)
 -- [X, Y, Z] =: sph2cart (S)
 -- C =: sph2cart (...)
     Transform spherical coordinates to Cartesian coordinates.

     The inputs THETA, PHI, and R must be the same shape, or scalar.  If called with a single matrix argument then each row of S represents the spherical coordinate (THETA, PHI, R).

     THETA describes the angle relative to the positive x-axis.

     PHI is the angle relative to the xy-plane.

     R is the distance to the origin (0, 0, 0).

     If only a single return argument is requested then return a matrix C where each row represents one Cartesian coordinate (X, Y, Z). cart2sph, pol2cart, cart2pol

Mathematical Constants
======================

 -- e:
 -- e: (N)
 -- e: (N, M)
 -- e: (N, M, K, ...)
 -- e: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the base of natural logarithms.

     The constant e satisfies the equation log (e) = 1.

     When called with no arguments, return a scalar with the value e.

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". log, exp, pi, I

 -- pi:
 -- pi: (N)
 -- pi: (N, M)
 -- pi: (N, M, K, ...)
 -- pi: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the ratio of the circumference of a circle to its diameter.

     Internally, pi is computed as 4.0 * atan (1.0).

     When called with no arguments, return a scalar with the value of pi.

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". e, I

 -- I:
 -- I: (N)
 -- I: (N, M)
 -- I: (N, M, K, ...)
 -- I: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the pure imaginary unit, defined as sqrt (-1).

     I, and its equivalents i, j, and J, are functions so any of the names may be reused for other purposes (such as i for a counter variable).

     When called with no arguments, return a scalar with the value i.

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". e, pi, log, exp

 -- Inf:
 -- Inf: (N)
 -- Inf: (N, M)
 -- Inf: (N, M, K, ...)
 -- Inf: (..., CLASS)
     Return a scalar, matrix or N-dimensional array whose elements are all equal to the IEEE representation for positive infinity.

     Infinity is produced when results are too large to be represented using the IEEE floating point format for numbers.  Two common examples which produce infinity are division by zero and overflow.

          [ 1/0 e^800 ]
          => Inf   Inf

     When called with no arguments, return a scalar with the value Inf.

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". isinf, NaN

 -- NaN:
 -- NaN: (N)
 -- NaN: (N, M)
 -- NaN: (N, M, K, ...)
 -- NaN: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the IEEE symbol NaN (Not a Number).

     NaN is the result of operations which do not produce a well defined numerical result.  Common operations which produce a NaN are arithmetic with infinity (Inf - Inf), zero divided by zero (0/0), and any operation involving another NaN value (5 + NaN).

     Note that NaN always compares not equal to NaN (NaN != NaN).  This behavior is specified by the IEEE standard for floating point arithmetic.  To find NaN values, use the isnan function.

     When called with no arguments, return a scalar with the value NaN.

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". isnan, Inf

 -- eps:
 -- eps: (X)
 -- eps: (N, M)
 -- eps: (N, M, K, ...)
 -- eps: (..., CLASS)
     Return a scalar, matrix or N-dimensional array whose elements are all eps, the machine precision.

     More precisely, eps is the relative spacing between any two adjacent numbers in the machine's floating point system.  This number is obviously system dependent.  On machines that support IEEE floating point arithmetic, eps is approximately 2.2204e-16 for double precision and 1.1921e-07 for single precision.

     When called with no arguments, return a scalar with the value eps (1.0).

     Given a single argument X, return the distance between X and the next largest value.

     When called with more than one argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.  The optional argument CLASS specifies the return type and may be either "double" or "single". realmax, realmin, intmax, flintmax

 -- realmax:
 -- realmax: (N)
 -- realmax: (N, M)
 -- realmax: (N, M, K, ...)
 -- realmax: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the largest floating point number that is representable.

     The actual value is system dependent.  On machines that support IEEE floating point arithmetic, realmax is approximately 1.7977e+308 for double precision and 3.4028e+38 for single precision.

     When called with no arguments, return a scalar with the value realmax ("double").

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". realmin, intmax, flintmax, eps

 -- realmin:
 -- realmin: (N)
 -- realmin: (N, M)
 -- realmin: (N, M, K, ...)
 -- realmin: (..., CLASS)
     Return a scalar, matrix, or N-dimensional array whose elements are all equal to the smallest normalized floating point number that is representable.

     The actual value is system dependent.  On machines that support IEEE floating point arithmetic, realmin is approximately 2.2251e-308 for double precision and 1.1755e-38 for single precision.

     When called with no arguments, return a scalar with the value realmin ("double").

     When called with a single argument, return a square matrix with the dimension specified.

     When called with more than one scalar argument the first two arguments are taken as the number of rows and columns and any further arguments specify additional matrix dimensions.

     The optional argument CLASS specifies the return type and may be either "double" or "single". realmax, intmin, eps