@c DO NOT EDIT!  Generated automatically by munge-texi.pl.

@c Copyright (C) 1996-2019 John W. Eaton
@c
@c This file is part of Octave.
@c
@c Octave is free software: you can redistribute it and/or modify it
@c under the terms of the GNU General Public License as published by
@c the Free Software Foundation, either version 3 of the License, or
@c (at your option) any later version.
@c
@c Octave is distributed in the hope that it will be useful, but
@c WITHOUT ANY WARRANTY; without even the implied warranty of
@c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
@c GNU General Public License for more details.
@c
@c You should have received a copy of the GNU General Public License
@c along with Octave; see the file COPYING.  If not, see
@c <https://www.gnu.org/licenses/>.

@node Arithmetic
@chapter 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 @dfn{mapping functions} apply the given operation individually to
each element when given a matrix argument.  For example:

@example
@group
sin ([1, 2; 3, 4])
     @result{}  0.84147   0.90930
         0.14112  -0.75680
@end group
@end example

@menu
* Exponents and Logarithms::
* Complex Arithmetic::
* Trigonometry::
* Sums and Products::
* Utility Functions::
* Special Functions::
* Rational Approximations::
* Coordinate Transformations::
* Mathematical Constants::
@end menu

@node Exponents and Logarithms
@section Exponents and Logarithms

@c exp libinterp/corefcn/mappers.cc
@anchor{XREFexp}
@deftypefn {} {} exp (@var{x})
Compute
@tex
$e^{x}$
@end tex
@ifnottex
@code{e^x}
@end ifnottex
for each element of @var{x}.

To compute the matrix exponential, see @ref{Linear Algebra}.
@xseealso{@ref{XREFlog,,log}}
@end deftypefn


@c expm1 libinterp/corefcn/mappers.cc
@anchor{XREFexpm1}
@deftypefn {} {} expm1 (@var{x})
Compute
@tex
$ e^{x} - 1 $
@end tex
@ifnottex
@code{exp (@var{x}) - 1}
@end ifnottex
accurately in the neighborhood of zero.
@xseealso{@ref{XREFexp,,exp}}
@end deftypefn


@c log libinterp/corefcn/mappers.cc
@anchor{XREFlog}
@deftypefn {} {} log (@var{x})
Compute the natural logarithm,
@tex
$\ln{(x)},$
@end tex
@ifnottex
@code{ln (@var{x})},
@end ifnottex
for each element of @var{x}.

To compute the matrix logarithm, see @ref{Linear Algebra}.
@xseealso{@ref{XREFexp,,exp}, @ref{XREFlog1p,,log1p}, @ref{XREFlog2,,log2}, @ref{XREFlog10,,log10}, @ref{XREFlogspace,,logspace}}
@end deftypefn


@c reallog scripts/specfun/reallog.m
@anchor{XREFreallog}
@deftypefn {} {} reallog (@var{x})
Return the real-valued natural logarithm of each element of @var{x}.

If any element results in a complex return value @code{reallog} aborts and
issues an error.
@xseealso{@ref{XREFlog,,log}, @ref{XREFrealpow,,realpow}, @ref{XREFrealsqrt,,realsqrt}}
@end deftypefn


@c log1p libinterp/corefcn/mappers.cc
@anchor{XREFlog1p}
@deftypefn {} {} log1p (@var{x})
Compute
@tex
$\ln{(1 + x)}$
@end tex
@ifnottex
@code{log (1 + @var{x})}
@end ifnottex
accurately in the neighborhood of zero.
@xseealso{@ref{XREFlog,,log}, @ref{XREFexp,,exp}, @ref{XREFexpm1,,expm1}}
@end deftypefn


@c log10 libinterp/corefcn/mappers.cc
@anchor{XREFlog10}
@deftypefn {} {} log10 (@var{x})
Compute the base-10 logarithm of each element of @var{x}.
@xseealso{@ref{XREFlog,,log}, @ref{XREFlog2,,log2}, @ref{XREFlogspace,,logspace}, @ref{XREFexp,,exp}}
@end deftypefn


@c log2 libinterp/corefcn/data.cc
@anchor{XREFlog2}
@deftypefn  {} {} log2 (@var{x})
@deftypefnx {} {[@var{f}, @var{e}] =} log2 (@var{x})
Compute the base-2 logarithm of each element of @var{x}.

If called with two output arguments, split @var{x} into
binary mantissa and exponent so that
@tex
${1 \over 2} \le \left| f \right| < 1$
@end tex
@ifnottex
@w{@code{1/2 <= abs(f) < 1}}
@end ifnottex
and @var{e} is an integer.  If
@tex
$x = 0$, $f = e = 0$.
@end tex
@ifnottex
@w{@code{x = 0}}, @w{@code{f = e = 0}}.
@end ifnottex
@xseealso{@ref{XREFpow2,,pow2}, @ref{XREFlog,,log}, @ref{XREFlog10,,log10}, @ref{XREFexp,,exp}}
@end deftypefn


@c pow2 scripts/specfun/pow2.m
@anchor{XREFpow2}
@deftypefn  {} {} pow2 (@var{x})
@deftypefnx {} {} pow2 (@var{f}, @var{e})
With one input argument, compute
@tex
$2^x$
@end tex
@ifnottex
2 .^ x
@end ifnottex
for each element of @var{x}.

With two input arguments, return
@tex
$f \cdot 2^e$.
@end tex
@ifnottex
f .* (2 .^ e).
@end ifnottex
@xseealso{@ref{XREFlog2,,log2}, @ref{XREFnextpow2,,nextpow2}, @ref{XREFpower,,power}}
@end deftypefn


@c nextpow2 scripts/general/nextpow2.m
@anchor{XREFnextpow2}
@deftypefn {} {@var{n} =} nextpow2 (@var{x})
Compute the exponent for the smallest power of two larger than the input.

For each element in the input array @var{x}, return the first integer
@var{n} such that
@tex
$2^n \ge |x|$.
@end tex
@ifnottex
2^n @geq{} abs (x).
@end ifnottex

@xseealso{@ref{XREFpow2,,pow2}, @ref{XREFlog2,,log2}}
@end deftypefn


@c realpow scripts/specfun/realpow.m
@anchor{XREFrealpow}
@deftypefn {} {} realpow (@var{x}, @var{y})
Compute the real-valued, element-by-element power operator.

This is equivalent to @w{@code{@var{x} .^ @var{y}}}, except that
@code{realpow} reports an error if any return value is complex.
@xseealso{@ref{XREFpower,,power}, @ref{XREFreallog,,reallog}, @ref{XREFrealsqrt,,realsqrt}}
@end deftypefn


@c sqrt libinterp/corefcn/mappers.cc
@anchor{XREFsqrt}
@deftypefn {} {} sqrt (@var{x})
Compute the square root of each element of @var{x}.

If @var{x} is negative, a complex result is returned.

To compute the matrix square root, see @ref{Linear Algebra}.
@xseealso{@ref{XREFrealsqrt,,realsqrt}, @ref{XREFnthroot,,nthroot}}
@end deftypefn


@c realsqrt scripts/specfun/realsqrt.m
@anchor{XREFrealsqrt}
@deftypefn {} {} realsqrt (@var{x})
Return the real-valued square root of each element of @var{x}.

If any element results in a complex return value @code{realsqrt} aborts and
issues an error.
@xseealso{@ref{XREFsqrt,,sqrt}, @ref{XREFrealpow,,realpow}, @ref{XREFreallog,,reallog}}
@end deftypefn


@c cbrt libinterp/corefcn/mappers.cc
@anchor{XREFcbrt}
@deftypefn {} {} cbrt (@var{x})
Compute the real cube root of each element of @var{x}.

Unlike @code{@var{x}^(1/3)}, the result will be negative if @var{x} is
negative.
@xseealso{@ref{XREFnthroot,,nthroot}}
@end deftypefn


@c nthroot scripts/specfun/nthroot.m
@anchor{XREFnthroot}
@deftypefn {} {} nthroot (@var{x}, @var{n})

Compute the real (non-complex) @var{n}-th root of @var{x}.

@var{x} must have all real entries and @var{n} must be a scalar.
If @var{n} is an even integer and @var{x} has negative entries then
@code{nthroot} aborts and issues an error.

Example:

@example
@group
nthroot (-1, 3)
@result{} -1
(-1) ^ (1 / 3)
@result{} 0.50000 - 0.86603i
@end group
@end example
@xseealso{@ref{XREFrealsqrt,,realsqrt}, @ref{XREFsqrt,,sqrt}, @ref{XREFcbrt,,cbrt}}
@end deftypefn


@node Complex Arithmetic
@section Complex Arithmetic

In the descriptions of the following functions,
@tex
$z$ is the complex number $x + iy$, where $i$ is defined as
$\sqrt{-1}$.
@end tex
@ifnottex
@var{z} is the complex number @var{x} + @var{i}@var{y}, where @var{i} is
defined as @code{sqrt (-1)}.
@end ifnottex

@c abs libinterp/corefcn/mappers.cc
@anchor{XREFabs}
@deftypefn {} {} abs (@var{z})
Compute the magnitude of @var{z}.

The magnitude is defined as
@tex
$|z| = \sqrt{x^2 + y^2}$.
@end tex
@ifnottex
|@var{z}| = @code{sqrt (x^2 + y^2)}.
@end ifnottex

For example:

@example
@group
abs (3 + 4i)
     @result{} 5
@end group
@end example
@xseealso{@ref{XREFarg,,arg}}
@end deftypefn


@c arg libinterp/corefcn/mappers.cc
@anchor{XREFarg}
@deftypefn  {} {} arg (@var{z})
@deftypefnx {} {} angle (@var{z})
Compute the argument, i.e., angle of @var{z}.

This is defined as,
@tex
$\theta = atan2 (y, x),$
@end tex
@ifnottex
@var{theta} = @code{atan2 (@var{y}, @var{x})},
@end ifnottex
in radians.

For example:

@example
@group
arg (3 + 4i)
     @result{} 0.92730
@end group
@end example
@xseealso{@ref{XREFabs,,abs}}
@end deftypefn


@c conj libinterp/corefcn/mappers.cc
@anchor{XREFconj}
@deftypefn {} {} conj (@var{z})
Return the complex conjugate of @var{z}.

The complex conjugate is defined as
@tex
$\bar{z} = x - iy$.
@end tex
@ifnottex
@code{conj (@var{z})} = @var{x} - @var{i}@var{y}.
@end ifnottex
@xseealso{@ref{XREFreal,,real}, @ref{XREFimag,,imag}}
@end deftypefn


@c cplxpair scripts/general/cplxpair.m
@anchor{XREFcplxpair}
@deftypefn  {} {} cplxpair (@var{z})
@deftypefnx {} {} cplxpair (@var{z}, @var{tol})
@deftypefnx {} {} cplxpair (@var{z}, @var{tol}, @var{dim})
Sort the numbers @var{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
@code{abs (imag (@var{z})) / abs (@var{z}) < @var{tol}}) are placed after
the complex pairs.

@var{tol} is a weighting factor in the range [0, 1) which determines the
tolerance of the matching.  The default value is @code{100 * eps} and the
resulting tolerance for a given complex pair is
@code{@var{tol} * abs (@var{z}(i)))}.

By default the complex pairs are sorted along the first non-singleton
dimension of @var{z}.  If @var{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 @var{tol}).
Note that there is no defined order for pairs with identical real parts but
differing imaginary parts.
@c Set example in small font to prevent overfull line

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


@c imag libinterp/corefcn/mappers.cc
@anchor{XREFimag}
@deftypefn {} {} imag (@var{z})
Return the imaginary part of @var{z} as a real number.
@xseealso{@ref{XREFreal,,real}, @ref{XREFconj,,conj}}
@end deftypefn


@c real libinterp/corefcn/mappers.cc
@anchor{XREFreal}
@deftypefn {} {} real (@var{z})
Return the real part of @var{z}.
@xseealso{@ref{XREFimag,,imag}, @ref{XREFconj,,conj}}
@end deftypefn


@node Trigonometry
@section Trigonometry

Octave provides the following trigonometric functions where angles are
specified in radians.  To convert from degrees to radians multiply by
@tex
$\pi/180$
@end tex
@ifnottex
@code{pi/180}
@end ifnottex
or use the @code{deg2rad} function.  For example, @code{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 @samp{d} suffix.  As an example, @code{sin} expects an angle in radians while
@code{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: @url{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
@tex
$atan2(\pm0, -0)$ returns $\pm \pi$.
@end tex
@ifnottex
@code{atan2(+- 0, 0)} returns @code{+- pi}.
@end ifnottex

It should be noted that @sc{matlab} uses different definitions which apparently
do not distinguish -0.

@c deg2rad scripts/general/deg2rad.m
@anchor{XREFdeg2rad}
@deftypefn {} {@var{rad} =} deg2rad (@var{deg})

Convert degrees to radians.

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

The output @var{rad} is the same size and shape as @var{deg} with degrees
converted to radians using the conversion constant @code{pi/180}.

Example:

@example
@group
deg2rad ([0, 90, 180, 270, 360])
  @result{}  0.00000   1.57080   3.14159   4.71239   6.28319
@end group
@end example
@xseealso{@ref{XREFrad2deg,,rad2deg}}
@end deftypefn

@c rad2deg scripts/general/rad2deg.m
@anchor{XREFrad2deg}
@deftypefn {} {@var{deg} =} rad2deg (@var{rad})

Convert radians to degrees.

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

The output @var{deg} is the same size and shape as @var{rad} with radians
converted to degrees using the conversion constant @code{180/pi}.

Example:

@example
@group
rad2deg ([0, pi/2, pi, 3/2*pi, 2*pi])
  @result{}  0    90   180   270   360
@end group
@end example
@xseealso{@ref{XREFdeg2rad,,deg2rad}}
@end deftypefn


@c sin libinterp/corefcn/mappers.cc
@anchor{XREFsin}
@deftypefn {} {} sin (@var{x})
Compute the sine for each element of @var{x} in radians.
@xseealso{@ref{XREFasin,,asin}, @ref{XREFsind,,sind}, @ref{XREFsinh,,sinh}}
@end deftypefn

@c cos libinterp/corefcn/mappers.cc
@anchor{XREFcos}
@deftypefn {} {} cos (@var{x})
Compute the cosine for each element of @var{x} in radians.
@xseealso{@ref{XREFacos,,acos}, @ref{XREFcosd,,cosd}, @ref{XREFcosh,,cosh}}
@end deftypefn

@c tan libinterp/corefcn/mappers.cc
@anchor{XREFtan}
@deftypefn {} {} tan (@var{z})
Compute the tangent for each element of @var{x} in radians.
@xseealso{@ref{XREFatan,,atan}, @ref{XREFtand,,tand}, @ref{XREFtanh,,tanh}}
@end deftypefn

@c sec scripts/elfun/sec.m
@anchor{XREFsec}
@deftypefn {} {} sec (@var{x})
Compute the secant for each element of @var{x} in radians.
@xseealso{@ref{XREFasec,,asec}, @ref{XREFsecd,,secd}, @ref{XREFsech,,sech}}
@end deftypefn

@c csc scripts/elfun/csc.m
@anchor{XREFcsc}
@deftypefn {} {} csc (@var{x})
Compute the cosecant for each element of @var{x} in radians.
@xseealso{@ref{XREFacsc,,acsc}, @ref{XREFcscd,,cscd}, @ref{XREFcsch,,csch}}
@end deftypefn

@c cot scripts/elfun/cot.m
@anchor{XREFcot}
@deftypefn {} {} cot (@var{x})
Compute the cotangent for each element of @var{x} in radians.
@xseealso{@ref{XREFacot,,acot}, @ref{XREFcotd,,cotd}, @ref{XREFcoth,,coth}}
@end deftypefn


@c asin libinterp/corefcn/mappers.cc
@anchor{XREFasin}
@deftypefn {} {} asin (@var{x})
Compute the inverse sine in radians for each element of @var{x}.
@xseealso{@ref{XREFsin,,sin}, @ref{XREFasind,,asind}}
@end deftypefn

@c acos libinterp/corefcn/mappers.cc
@anchor{XREFacos}
@deftypefn {} {} acos (@var{x})
Compute the inverse cosine in radians for each element of @var{x}.
@xseealso{@ref{XREFcos,,cos}, @ref{XREFacosd,,acosd}}
@end deftypefn

@c atan libinterp/corefcn/mappers.cc
@anchor{XREFatan}
@deftypefn {} {} atan (@var{x})
Compute the inverse tangent in radians for each element of @var{x}.
@xseealso{@ref{XREFtan,,tan}, @ref{XREFatand,,atand}}
@end deftypefn

@c asec scripts/elfun/asec.m
@anchor{XREFasec}
@deftypefn {} {} asec (@var{x})
Compute the inverse secant in radians for each element of @var{x}.
@xseealso{@ref{XREFsec,,sec}, @ref{XREFasecd,,asecd}}
@end deftypefn

@c acsc scripts/elfun/acsc.m
@anchor{XREFacsc}
@deftypefn {} {} acsc (@var{x})
Compute the inverse cosecant in radians for each element of @var{x}.
@xseealso{@ref{XREFcsc,,csc}, @ref{XREFacscd,,acscd}}
@end deftypefn

@c acot scripts/elfun/acot.m
@anchor{XREFacot}
@deftypefn {} {} acot (@var{x})
Compute the inverse cotangent in radians for each element of @var{x}.
@xseealso{@ref{XREFcot,,cot}, @ref{XREFacotd,,acotd}}
@end deftypefn


@c sinh libinterp/corefcn/mappers.cc
@anchor{XREFsinh}
@deftypefn {} {} sinh (@var{x})
Compute the hyperbolic sine for each element of @var{x}.
@xseealso{@ref{XREFasinh,,asinh}, @ref{XREFcosh,,cosh}, @ref{XREFtanh,,tanh}}
@end deftypefn

@c cosh libinterp/corefcn/mappers.cc
@anchor{XREFcosh}
@deftypefn {} {} cosh (@var{x})
Compute the hyperbolic cosine for each element of @var{x}.
@xseealso{@ref{XREFacosh,,acosh}, @ref{XREFsinh,,sinh}, @ref{XREFtanh,,tanh}}
@end deftypefn

@c tanh libinterp/corefcn/mappers.cc
@anchor{XREFtanh}
@deftypefn {} {} tanh (@var{x})
Compute hyperbolic tangent for each element of @var{x}.
@xseealso{@ref{XREFatanh,,atanh}, @ref{XREFsinh,,sinh}, @ref{XREFcosh,,cosh}}
@end deftypefn

@c sech scripts/elfun/sech.m
@anchor{XREFsech}
@deftypefn {} {} sech (@var{x})
Compute the hyperbolic secant of each element of @var{x}.
@xseealso{@ref{XREFasech,,asech}}
@end deftypefn

@c csch scripts/elfun/csch.m
@anchor{XREFcsch}
@deftypefn {} {} csch (@var{x})
Compute the hyperbolic cosecant of each element of @var{x}.
@xseealso{@ref{XREFacsch,,acsch}}
@end deftypefn

@c coth scripts/elfun/coth.m
@anchor{XREFcoth}
@deftypefn {} {} coth (@var{x})
Compute the hyperbolic cotangent of each element of @var{x}.
@xseealso{@ref{XREFacoth,,acoth}}
@end deftypefn


@c asinh libinterp/corefcn/mappers.cc
@anchor{XREFasinh}
@deftypefn {} {} asinh (@var{x})
Compute the inverse hyperbolic sine for each element of @var{x}.
@xseealso{@ref{XREFsinh,,sinh}}
@end deftypefn

@c acosh libinterp/corefcn/mappers.cc
@anchor{XREFacosh}
@deftypefn {} {} acosh (@var{x})
Compute the inverse hyperbolic cosine for each element of @var{x}.
@xseealso{@ref{XREFcosh,,cosh}}
@end deftypefn

@c atanh libinterp/corefcn/mappers.cc
@anchor{XREFatanh}
@deftypefn {} {} atanh (@var{x})
Compute the inverse hyperbolic tangent for each element of @var{x}.
@xseealso{@ref{XREFtanh,,tanh}}
@end deftypefn

@c asech scripts/elfun/asech.m
@anchor{XREFasech}
@deftypefn {} {} asech (@var{x})
Compute the inverse hyperbolic secant of each element of @var{x}.
@xseealso{@ref{XREFsech,,sech}}
@end deftypefn

@c acsch scripts/elfun/acsch.m
@anchor{XREFacsch}
@deftypefn {} {} acsch (@var{x})
Compute the inverse hyperbolic cosecant of each element of @var{x}.
@xseealso{@ref{XREFcsch,,csch}}
@end deftypefn

@c acoth scripts/elfun/acoth.m
@anchor{XREFacoth}
@deftypefn {} {} acoth (@var{x})
Compute the inverse hyperbolic cotangent of each element of @var{x}.
@xseealso{@ref{XREFcoth,,coth}}
@end deftypefn


@c atan2 libinterp/corefcn/data.cc
@anchor{XREFatan2}
@deftypefn {} {} atan2 (@var{y}, @var{x})
Compute atan (@var{y} / @var{x}) for corresponding elements of @var{y}
and @var{x}.

@var{y} and @var{x} must match in size and orientation.  The signs of
elements of @var{y} and @var{x} are used to determine the quadrants of each
resulting value.

This function is equivalent to @code{arg (complex (@var{x}, @var{y}))}.
@xseealso{@ref{XREFtan,,tan}, @ref{XREFtand,,tand}, @ref{XREFtanh,,tanh}, @ref{XREFatanh,,atanh}}
@end deftypefn


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:

@example
@group
cosd (90)
     @result{} 0
cos (pi/2)
     @result{} 6.1230e-17
@end group
@end example

@c sind scripts/elfun/sind.m
@anchor{XREFsind}
@deftypefn {} {} sind (@var{x})
Compute the sine for each element of @var{x} in degrees.

Returns zero for elements where @code{@var{x}/180} is an integer.
@xseealso{@ref{XREFasind,,asind}, @ref{XREFsin,,sin}}
@end deftypefn

@c cosd scripts/elfun/cosd.m
@anchor{XREFcosd}
@deftypefn {} {} cosd (@var{x})
Compute the cosine for each element of @var{x} in degrees.

Returns zero for elements where @code{(@var{x}-90)/180} is an integer.
@xseealso{@ref{XREFacosd,,acosd}, @ref{XREFcos,,cos}}
@end deftypefn

@c tand scripts/elfun/tand.m
@anchor{XREFtand}
@deftypefn {} {} tand (@var{x})
Compute the tangent for each element of @var{x} in degrees.

Returns zero for elements where @code{@var{x}/180} is an integer and
@code{Inf} for elements where @code{(@var{x}-90)/180} is an integer.
@xseealso{@ref{XREFatand,,atand}, @ref{XREFtan,,tan}}
@end deftypefn

@c secd scripts/elfun/secd.m
@anchor{XREFsecd}
@deftypefn {} {} secd (@var{x})
Compute the secant for each element of @var{x} in degrees.
@xseealso{@ref{XREFasecd,,asecd}, @ref{XREFsec,,sec}}
@end deftypefn

@c cscd scripts/elfun/cscd.m
@anchor{XREFcscd}
@deftypefn {} {} cscd (@var{x})
Compute the cosecant for each element of @var{x} in degrees.
@xseealso{@ref{XREFacscd,,acscd}, @ref{XREFcsc,,csc}}
@end deftypefn

@c cotd scripts/elfun/cotd.m
@anchor{XREFcotd}
@deftypefn {} {} cotd (@var{x})
Compute the cotangent for each element of @var{x} in degrees.
@xseealso{@ref{XREFacotd,,acotd}, @ref{XREFcot,,cot}}
@end deftypefn


@c asind scripts/elfun/asind.m
@anchor{XREFasind}
@deftypefn {} {} asind (@var{x})
Compute the inverse sine in degrees for each element of @var{x}.
@xseealso{@ref{XREFsind,,sind}, @ref{XREFasin,,asin}}
@end deftypefn

@c acosd scripts/elfun/acosd.m
@anchor{XREFacosd}
@deftypefn {} {} acosd (@var{x})
Compute the inverse cosine in degrees for each element of @var{x}.
@xseealso{@ref{XREFcosd,,cosd}, @ref{XREFacos,,acos}}
@end deftypefn

@c atand scripts/elfun/atand.m
@anchor{XREFatand}
@deftypefn {} {} atand (@var{x})
Compute the inverse tangent in degrees for each element of @var{x}.
@xseealso{@ref{XREFtand,,tand}, @ref{XREFatan,,atan}}
@end deftypefn

@c atan2d scripts/elfun/atan2d.m
@anchor{XREFatan2d}
@deftypefn {} {} atan2d (@var{y}, @var{x})
Compute atan (@var{y} / @var{x}) in degrees for corresponding elements
from @var{y} and @var{x}.
@xseealso{@ref{XREFtand,,tand}, @ref{XREFatan2,,atan2}}
@end deftypefn

@c asecd scripts/elfun/asecd.m
@anchor{XREFasecd}
@deftypefn {} {} asecd (@var{x})
Compute the inverse secant in degrees for each element of @var{x}.
@xseealso{@ref{XREFsecd,,secd}, @ref{XREFasec,,asec}}
@end deftypefn

@c acscd scripts/elfun/acscd.m
@anchor{XREFacscd}
@deftypefn {} {} acscd (@var{x})
Compute the inverse cosecant in degrees for each element of @var{x}.
@xseealso{@ref{XREFcscd,,cscd}, @ref{XREFacsc,,acsc}}
@end deftypefn

@c acotd scripts/elfun/acotd.m
@anchor{XREFacotd}
@deftypefn {} {} acotd (@var{x})
Compute the inverse cotangent in degrees for each element of @var{x}.
@xseealso{@ref{XREFcotd,,cotd}, @ref{XREFacot,,acot}}
@end deftypefn


@node Sums and Products
@section Sums and Products

@c sum libinterp/corefcn/data.cc
@anchor{XREFsum}
@deftypefn  {} {} sum (@var{x})
@deftypefnx {} {} sum (@var{x}, @var{dim})
@deftypefnx {} {} sum (@dots{}, "native")
@deftypefnx {} {} sum (@dots{}, "double")
@deftypefnx {} {} sum (@dots{}, "extra")
Sum of elements along dimension @var{dim}.

If @var{dim} is omitted, it defaults to the first non-singleton dimension.

The optional @qcode{"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 @qcode{"native"} is given, then the operation is performed in the same
type as the original argument.

For example:

@example
@group
sum ([true, true])
   @result{} 2
sum ([true, true], "native")
   @result{} true
@end group
@end example

If @qcode{"double"} is given the sum is performed in double precision even for
single precision inputs.

For double precision inputs, the @qcode{"extra"} option will use a more
accurate algorithm than straightforward summation.  For single precision
inputs, @qcode{"extra"} is the same as @qcode{"double"}.  For all other data
type @qcode{"extra"} has no effect.
@xseealso{@ref{XREFcumsum,,cumsum}, @ref{XREFsumsq,,sumsq}, @ref{XREFprod,,prod}}
@end deftypefn


@c prod libinterp/corefcn/data.cc
@anchor{XREFprod}
@deftypefn  {} {} prod (@var{x})
@deftypefnx {} {} prod (@var{x}, @var{dim})
@deftypefnx {} {} prod (@dots{}, "native")
@deftypefnx {} {} prod (@dots{}, "double")
Product of elements along dimension @var{dim}.

If @var{dim} is omitted, it defaults to the first non-singleton dimension.

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

For example:

@example
@group
prod ([true, true])
   @result{} 1
prod ([true, true], "native")
   @result{} true
@end group
@end example

On the contrary, if @qcode{"double"} is given, the operation is performed
in double precision even for single precision inputs.
@xseealso{@ref{XREFcumprod,,cumprod}, @ref{XREFsum,,sum}}
@end deftypefn


@c cumsum libinterp/corefcn/data.cc
@anchor{XREFcumsum}
@deftypefn  {} {} cumsum (@var{x})
@deftypefnx {} {} cumsum (@var{x}, @var{dim})
@deftypefnx {} {} cumsum (@dots{}, "native")
@deftypefnx {} {} cumsum (@dots{}, "double")
Cumulative sum of elements along dimension @var{dim}.

If @var{dim} is omitted, it defaults to the first non-singleton dimension.
For example:

@example
@group
cumsum ([1, 2; 3, 4; 5, 6])
   @result{}  1   2
       4   6
       9  12
@end group
@end example

See @code{sum} for an explanation of the optional parameters @qcode{"native"}
and @qcode{"double"}.
@xseealso{@ref{XREFsum,,sum}, @ref{XREFcumprod,,cumprod}}
@end deftypefn


@c cumprod libinterp/corefcn/data.cc
@anchor{XREFcumprod}
@deftypefn  {} {} cumprod (@var{x})
@deftypefnx {} {} cumprod (@var{x}, @var{dim})
Cumulative product of elements along dimension @var{dim}.

If @var{dim} is omitted, it defaults to the first non-singleton dimension.
For example:

@example
@group
cumprod ([1, 2; 3, 4; 5, 6])
   @result{}  1   2
       3   8
      15  48
@end group
@end example
@xseealso{@ref{XREFprod,,prod}, @ref{XREFcumsum,,cumsum}}
@end deftypefn


@c sumsq libinterp/corefcn/data.cc
@anchor{XREFsumsq}
@deftypefn  {} {} sumsq (@var{x})
@deftypefnx {} {} sumsq (@var{x}, @var{dim})
Sum of squares of elements along dimension @var{dim}.

If @var{dim} is omitted, it defaults to the first non-singleton dimension.

This function is conceptually equivalent to computing

@example
sum (x .* conj (x), dim)
@end example

@noindent
but it uses less memory and avoids calling @code{conj} if @var{x} is real.
@xseealso{@ref{XREFsum,,sum}, @ref{XREFprod,,prod}}
@end deftypefn


@node Utility Functions
@section Utility Functions

@c ceil libinterp/corefcn/mappers.cc
@anchor{XREFceil}
@deftypefn {} {} ceil (@var{x})
Return the smallest integer not less than @var{x}.

This is equivalent to rounding towards positive infinity.

If @var{x} is complex, return
@code{ceil (real (@var{x})) + ceil (imag (@var{x})) * I}.

@example
@group
ceil ([-2.7, 2.7])
    @result{} -2    3
@end group
@end example
@xseealso{@ref{XREFfloor,,floor}, @ref{XREFround,,round}, @ref{XREFfix,,fix}}
@end deftypefn


@c fix libinterp/corefcn/mappers.cc
@anchor{XREFfix}
@deftypefn {} {} fix (@var{x})
Truncate fractional portion of @var{x} and return the integer portion.

This is equivalent to rounding towards zero.  If @var{x} is complex, return
@code{fix (real (@var{x})) + fix (imag (@var{x})) * I}.

@example
@group
fix ([-2.7, 2.7])
   @result{} -2    2
@end group
@end example
@xseealso{@ref{XREFceil,,ceil}, @ref{XREFfloor,,floor}, @ref{XREFround,,round}}
@end deftypefn


@c floor libinterp/corefcn/mappers.cc
@anchor{XREFfloor}
@deftypefn {} {} floor (@var{x})
Return the largest integer not greater than @var{x}.

This is equivalent to rounding towards negative infinity.  If @var{x} is
complex, return @code{floor (real (@var{x})) + floor (imag (@var{x})) * I}.

@example
@group
floor ([-2.7, 2.7])
     @result{} -3    2
@end group
@end example
@xseealso{@ref{XREFceil,,ceil}, @ref{XREFround,,round}, @ref{XREFfix,,fix}}
@end deftypefn


@c round libinterp/corefcn/mappers.cc
@anchor{XREFround}
@deftypefn {} {} round (@var{x})
Return the integer nearest to @var{x}.

If @var{x} is complex, return
@code{round (real (@var{x})) + round (imag (@var{x})) * I}.  If there
are two nearest integers, return the one further away from zero.

@example
@group
round ([-2.7, 2.7])
     @result{} -3    3
@end group
@end example
@xseealso{@ref{XREFceil,,ceil}, @ref{XREFfloor,,floor}, @ref{XREFfix,,fix}, @ref{XREFroundb,,roundb}}
@end deftypefn


@c roundb libinterp/corefcn/mappers.cc
@anchor{XREFroundb}
@deftypefn {} {} roundb (@var{x})
Return the integer nearest to @var{x}.  If there are two nearest
integers, return the even one (banker's rounding).

If @var{x} is complex,
return @code{roundb (real (@var{x})) + roundb (imag (@var{x})) * I}.
@xseealso{@ref{XREFround,,round}}
@end deftypefn


@c max libinterp/corefcn/max.cc
@anchor{XREFmax}
@deftypefn  {} {} max (@var{x})
@deftypefnx {} {} max (@var{x}, [], @var{dim})
@deftypefnx {} {[@var{w}, @var{iw}] =} max (@var{x})
@deftypefnx {} {} max (@var{x}, @var{y})
Find maximum values in the array @var{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, @code{max} operates along the first non-singleton
dimension.

If the optional third argument @var{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 (@var{x} and @var{y}), return the pairwise maximum according to
the rules for @ref{Broadcasting}.

Thus,

@example
max (max (@var{x}))
@end example

@noindent
returns the largest element of the 2-D matrix @var{x}, and

@example
@group
max (2:5, pi)
    @result{}  3.1416  3.1416  4.0000  5.0000
@end group
@end example

@noindent
compares each element of the range @code{2:5} with @code{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,

@example
@group
max ([-1 i 1 -i])
    @result{} -1
@end group
@end example

@noindent
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, @code{max} also returns
the first index of the maximum value(s).  Thus,

@example
@group
[x, ix] = max ([1, 3, 5, 2, 5])
    @result{}  x = 5
        ix = 3
@end group
@end example
@xseealso{@ref{XREFmin,,min}, @ref{XREFcummax,,cummax}, @ref{XREFcummin,,cummin}}
@end deftypefn


@c min libinterp/corefcn/max.cc
@anchor{XREFmin}
@deftypefn  {} {} min (@var{x})
@deftypefnx {} {} min (@var{x}, [], @var{dim})
@deftypefnx {} {[@var{w}, @var{iw}] =} min (@var{x})
@deftypefnx {} {} min (@var{x}, @var{y})
Find minimum values in the array @var{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, @code{min} operates along the first non-singleton
dimension.

If the optional third argument @var{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 (@var{x} and @var{y}), return the pairwise minimum according to
the rules for @ref{Broadcasting}.

Thus,

@example
min (min (@var{x}))
@end example

@noindent
returns the smallest element of the 2-D matrix @var{x}, and

@example
@group
min (2:5, pi)
    @result{}  2.0000  3.0000  3.1416  3.1416
@end group
@end example

@noindent
compares each element of the range @code{2:5} with @code{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,

@example
@group
min ([-1 i 1 -i])
    @result{} -i
@end group
@end example

@noindent
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, @code{min} also returns
the first index of the minimum value(s).  Thus,

@example
@group
[x, ix] = min ([1, 3, 0, 2, 0])
    @result{}  x = 0
        ix = 3
@end group
@end example
@xseealso{@ref{XREFmax,,max}, @ref{XREFcummin,,cummin}, @ref{XREFcummax,,cummax}}
@end deftypefn


@c cummax libinterp/corefcn/max.cc
@anchor{XREFcummax}
@deftypefn  {} {} cummax (@var{x})
@deftypefnx {} {} cummax (@var{x}, @var{dim})
@deftypefnx {} {[@var{w}, @var{iw}] =} cummax (@dots{})
Return the cumulative maximum values along dimension @var{dim}.

If @var{dim} is unspecified it defaults to column-wise operation.  For
example:

@example
@group
cummax ([1 3 2 6 4 5])
   @result{}  1  3  3  6  6  6
@end group
@end example

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

@example
@group
[w, iw] = cummax ([1 3 2 6 4 5])
@result{}
w =  1  3  3  6  6  6
iw = 1  2  2  4  4  4
@end group
@end example

@xseealso{@ref{XREFcummin,,cummin}, @ref{XREFmax,,max}, @ref{XREFmin,,min}}
@end deftypefn


@c cummin libinterp/corefcn/max.cc
@anchor{XREFcummin}
@deftypefn  {} {} cummin (@var{x})
@deftypefnx {} {} cummin (@var{x}, @var{dim})
@deftypefnx {} {[@var{w}, @var{iw}] =} cummin (@var{x})
Return the cumulative minimum values along dimension @var{dim}.

If @var{dim} is unspecified it defaults to column-wise operation.  For
example:

@example
@group
cummin ([5 4 6 2 3 1])
   @result{}  5  4  4  2  2  1
@end group
@end example

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

@example
@group
[w, iw] = cummin ([5 4 6 2 3 1])
@result{}
w =  5  4  4  2  2  1
iw = 1  2  2  4  4  6
@end group
@end example

@xseealso{@ref{XREFcummax,,cummax}, @ref{XREFmin,,min}, @ref{XREFmax,,max}}
@end deftypefn


@c hypot libinterp/corefcn/data.cc
@anchor{XREFhypot}
@deftypefn  {} {} hypot (@var{x}, @var{y})
@deftypefnx {} {} hypot (@var{x}, @var{y}, @var{z}, @dots{})
Compute the element-by-element square root of the sum of the squares of
@var{x} and @var{y}.

This is equivalent to
@code{sqrt (@var{x}.^2 + @var{y}.^2)}, but is calculated in a manner that
avoids overflows for large values of @var{x} or @var{y}.

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

@example
@group
hypot (hypot (@var{x}, @var{y}), @var{z})
hypot (hypot (hypot (@var{x}, @var{y}), @var{z}), @var{w}), etc.
@end group
@end example
@end deftypefn


@c gradient scripts/general/gradient.m
@anchor{XREFgradient}
@deftypefn  {} {@var{dx} =} gradient (@var{m})
@deftypefnx {} {[@var{dx}, @var{dy}, @var{dz}, @dots{}] =} gradient (@var{m})
@deftypefnx {} {[@dots{}] =} gradient (@var{m}, @var{s})
@deftypefnx {} {[@dots{}] =} gradient (@var{m}, @var{x}, @var{y}, @var{z}, @dots{})
@deftypefnx {} {[@dots{}] =} gradient (@var{f}, @var{x0})
@deftypefnx {} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{s})
@deftypefnx {} {[@dots{}] =} gradient (@var{f}, @var{x0}, @var{x}, @var{y}, @dots{})

Calculate the gradient of sampled data or a function.

If @var{m} is a vector, calculate the one-dimensional gradient of @var{m}.
If @var{m} is a matrix the gradient is calculated for each dimension.

@code{[@var{dx}, @var{dy}] = gradient (@var{m})} calculates the
one-dimensional gradient for @var{x} and @var{y} direction if @var{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 @var{s}
parameter.  If @var{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 @var{x}, @dots{} arguments.  Scalar values specify an equidistant
spacing.  Vector values for the @var{x}, @dots{} arguments specify the
coordinate for that dimension.  The length must match their respective
dimension of @var{m}.

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

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

If the first argument @var{f} is a function handle, the gradient of the
function at the points in @var{x0} is approximated using central difference.
For example, @code{gradient (@@cos, 0)} approximates the gradient of the
cosine function in the point @math{x0 = 0}.  As with sampled data, the
spacing values between the points from which the gradient is estimated can
be set via the @var{s} or @var{dx}, @var{dy}, @dots{} arguments.  By default
a spacing of 1 is used.
@xseealso{@ref{XREFdiff,,diff}, @ref{XREFdel2,,del2}}
@end deftypefn


@c dot libinterp/corefcn/dot.cc
@anchor{XREFdot}
@deftypefn {} {} dot (@var{x}, @var{y}, @var{dim})
Compute the dot product of two vectors.

If @var{x} and @var{y} are matrices, calculate the dot products along the
first non-singleton dimension.

If the optional argument @var{dim} is given, calculate the dot products
along this dimension.

Implementation Note: This is equivalent to
@code{sum (conj (@var{X}) .* @var{Y}, @var{dim})}, but avoids forming a
temporary array and is faster.  When @var{X} and @var{Y} are column vectors,
the result is equivalent to @code{@var{X}' * @var{Y}}.  Although, @code{dot}
is defined for integer arrays, the output may differ from the expected result
due to the limited range of integer objects.
@xseealso{@ref{XREFcross,,cross}, @ref{XREFdivergence,,divergence}}
@end deftypefn


@c cross scripts/linear-algebra/cross.m
@anchor{XREFcross}
@deftypefn  {} {} cross (@var{x}, @var{y})
@deftypefnx {} {} cross (@var{x}, @var{y}, @var{dim})
Compute the vector cross product of two 3-dimensional vectors @var{x} and
@var{y}.

If @var{x} and @var{y} are matrices, the cross product is applied along the
first dimension with three elements.

The optional argument  @var{dim} forces the cross product to be calculated
along the specified dimension.

Example Code:

@example
@group
cross ([1, 1, 0], [0, 1, 1])
  @result{}
       1   -1   1
@end group
@end example

@xseealso{@ref{XREFdot,,dot}, @ref{XREFcurl,,curl}, @ref{XREFdivergence,,divergence}}
@end deftypefn


@c divergence scripts/general/divergence.m
@anchor{XREFdivergence}
@deftypefn  {} {@var{div} =} divergence (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz})
@deftypefnx {} {@var{div} =} divergence (@var{fx}, @var{fy}, @var{fz})
@deftypefnx {} {@var{div} =} divergence (@var{x}, @var{y}, @var{fx}, @var{fy})
@deftypefnx {} {@var{div} =} divergence (@var{fx}, @var{fy})
Calculate divergence of a vector field given by the arrays @var{fx},
@var{fy}, and @var{fz} or @var{fx}, @var{fy} respectively.

@tex
$$
div F(x,y,z) = \partial_x{F} + \partial_y{F} + \partial_z{F}
$$
@end tex
@ifnottex

@example
@group
                  d               d               d
div F(x,y,z)  =   -- F(x,y,z)  +  -- F(x,y,z)  +  -- F(x,y,z)
                  dx              dy              dz
@end group
@end example

@end ifnottex
The coordinates of the vector field can be given by the arguments @var{x},
@var{y}, @var{z} or @var{x}, @var{y} respectively.

@xseealso{@ref{XREFcurl,,curl}, @ref{XREFgradient,,gradient}, @ref{XREFdel2,,del2}, @ref{XREFdot,,dot}}
@end deftypefn


@c curl scripts/general/curl.m
@anchor{XREFcurl}
@deftypefn  {} {[@var{cx}, @var{cy}, @var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{z}, @var{fx}, @var{fy}, @var{fz})
@deftypefnx {} {[@var{cz}, @var{v}] =} curl (@var{x}, @var{y}, @var{fx}, @var{fy})
@deftypefnx {} {[@dots{}] =} curl (@var{fx}, @var{fy}, @var{fz})
@deftypefnx {} {[@dots{}] =} curl (@var{fx}, @var{fy})
@deftypefnx {} {@var{v} =} curl (@dots{})
Calculate curl of vector field given by the arrays @var{fx}, @var{fy}, and
@var{fz} or @var{fx}, @var{fy} respectively.
@tex
$$ curl F(x,y,z) = \left( {\partial{d} \over \partial{y}} F_z - {\partial{d} \over \partial{z}} F_y, {\partial{d} \over \partial{z}} F_x - {\partial{d} \over \partial{x}} F_z, {\partial{d} \over \partial{x}} F_y - {\partial{d} \over \partial{y}} F_x \right)$$
@end tex
@ifnottex

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

@end ifnottex
The coordinates of the vector field can be given by the arguments @var{x},
@var{y}, @var{z} or @var{x}, @var{y} respectively.  @var{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.
@xseealso{@ref{XREFdivergence,,divergence}, @ref{XREFgradient,,gradient}, @ref{XREFdel2,,del2}, @ref{XREFcross,,cross}}
@end deftypefn


@c del2 scripts/general/del2.m
@anchor{XREFdel2}
@deftypefn  {} {@var{L} =} del2 (@var{M})
@deftypefnx {} {@var{L} =} del2 (@var{M}, @var{h})
@deftypefnx {} {@var{L} =} del2 (@var{M}, @var{dx}, @var{dy}, @dots{})

Calculate the discrete Laplace
@tex
operator $( \nabla^2 )$.
@end tex
@ifnottex
operator.
@end ifnottex

For a 2-dimensional matrix @var{M} this is defined as
@tex
$$L = {1 \over 4} \left( {d^2 \over dx^2} M(x,y) + {d^2 \over dy^2} M(x,y) \right)$$
@end tex
@ifnottex

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

@end ifnottex
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 @var{h}, which is a
scalar defining the equidistant spacing in all dimensions.  Alternatively,
the spacing in each dimension may be defined separately by @var{dx},
@var{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
@var{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

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

@xseealso{@ref{XREFgradient,,gradient}, @ref{XREFdiff,,diff}}
@end deftypefn


@c factorial scripts/specfun/factorial.m
@anchor{XREFfactorial}
@deftypefn {} {} factorial (@var{n})
Return the factorial of @var{n} where @var{n} is a real non-negative
integer.

If @var{n} is a scalar, this is equivalent to @code{prod (1:@var{n})}.  For
vector or matrix arguments, return the factorial of each element in the
array.

For non-integers see the generalized factorial function @code{gamma}.
Note that the factorial function grows large quite quickly, and even
with double precision values overflow will occur if @var{n} > 171.  For
such cases consider @code{gammaln}.
@xseealso{@ref{XREFprod,,prod}, @ref{XREFgamma,,gamma}, @ref{XREFgammaln,,gammaln}}
@end deftypefn


@c factor scripts/specfun/factor.m
@anchor{XREFfactor}
@deftypefn  {} {@var{pf} =} factor (@var{q})
@deftypefnx {} {[@var{pf}, @var{n}] =} factor (@var{q})
Return the prime factorization of @var{q}.

The prime factorization is defined as @code{prod (@var{pf}) == @var{q}}
where every element of @var{pf} is a prime number.  If @code{@var{q} == 1},
return 1.  The output @var{pf} is of the same numeric class as the input.

With two output arguments, return the unique prime factors @var{pf} and
their multiplicities.  That is,
@code{prod (@var{pf} .^ @var{n}) == @var{q}}.

Implementation Note: The input @var{q} must be less than @code{flintmax}
(9.0072e+15) in order to factor correctly.
@xseealso{@ref{XREFgcd,,gcd}, @ref{XREFlcm,,lcm}, @ref{XREFisprime,,isprime}, @ref{XREFprimes,,primes}}
@end deftypefn


@c gcd libinterp/corefcn/gcd.cc
@anchor{XREFgcd}
@deftypefn  {} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})
@deftypefnx {} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @var{a2}, @dots{})
Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}.

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 @var{v1}, @dots{}, contain integer vectors such
that,

@tex
$g = v_1 a_1 + v_2 a_2 + \cdots$
@end tex
@ifnottex

@example
@var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}
@end example

@end ifnottex

Example code:

@example
@group
gcd ([15, 9], [20, 18])
   @result{}  5  9
@end group
@end example

@xseealso{@ref{XREFlcm,,lcm}, @ref{XREFfactor,,factor}, @ref{XREFisprime,,isprime}}
@end deftypefn


@c lcm scripts/specfun/lcm.m
@anchor{XREFlcm}
@deftypefn  {} {} lcm (@var{x}, @var{y})
@deftypefnx {} {} lcm (@var{x}, @var{y}, @dots{})
Compute the least common multiple of @var{x} and @var{y}, or of the list of
all arguments.

All elements must be numeric and of the same size or scalar.
@xseealso{@ref{XREFfactor,,factor}, @ref{XREFgcd,,gcd}, @ref{XREFisprime,,isprime}}
@end deftypefn


@c rem libinterp/corefcn/data.cc
@anchor{XREFrem}
@deftypefn {} {} rem (@var{x}, @var{y})
Return the remainder of the division @code{@var{x} / @var{y}}.

The remainder is computed using the expression

@example
x - y .* fix (x ./ y)
@end example

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 @sc{matlab}.  Any floating
point integers greater than @code{flintmax} (2^53 for double) will not compute
correctly.  For larger integer values convert the input to @code{uint64} before
calling this function.

By convention,

@example
@group
rem (@var{x}, 0) = NaN  if @var{x} is a floating point variable
rem (@var{x}, 0) = 0    if @var{x} is an integer variable
rem (@var{x}, @var{y})  returns a value with the signbit from @var{x}
@end group
@end example

For the opposite conventions see the @code{mod} function.  In general,
@code{rem} is best when computing the remainder after division of two
@emph{positive} numbers.  For negative numbers, or when the values are
periodic, @code{mod} is a better choice.
@xseealso{@ref{XREFmod,,mod}}
@end deftypefn


@c mod libinterp/corefcn/data.cc
@anchor{XREFmod}
@deftypefn {} {} mod (@var{x}, @var{y})
Compute the modulo of @var{x} and @var{y}.

Conceptually this is given by

@example
x - y .* floor (x ./ y)
@end example

@noindent
and is written such that the correct modulus is returned for integer types.
This function handles negative values correctly.  That is,
@w{@code{mod (-1, 3)}} is 2, not -1, as @w{@code{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 @sc{matlab}.  Any floating
point integers greater than @code{flintmax} (2^53 for double) will not compute
correctly.  For larger integer values convert the input to @code{uint64} before
calling this function.

By convention,

@example
@group
mod (@var{x}, 0) = @var{x}
mod (@var{x}, @var{y})      returns a value with the signbit from @var{y}
@end group
@end example

For the opposite conventions see the @code{rem} function.  In general,
@code{mod} is a better choice than @code{rem} when any of the inputs are
negative numbers or when the values are periodic.
@xseealso{@ref{XREFrem,,rem}}
@end deftypefn


@c primes scripts/specfun/primes.m
@anchor{XREFprimes}
@deftypefn {} {@var{p} =} primes (@var{n})
Return all primes up to @var{n}.

The output data class (double, single, uint32, etc.@:) is the same as the
input class of @var{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
@math{k} primes less than
@tex
$k \log (5 k)$.
@end tex
@ifnottex
k*log (5*k).
@end ifnottex

See also @code{list_primes} if you need a specific number @var{n} of primes.
@xseealso{@ref{XREFlist_primes,,list_primes}, @ref{XREFisprime,,isprime}}
@end deftypefn


@c list_primes scripts/miscellaneous/list_primes.m
@anchor{XREFlist_primes}
@deftypefn  {} {} list_primes ()
@deftypefnx {} {} list_primes (@var{n})
List the first @var{n} primes.

If @var{n} is unspecified, the first 25 primes are listed.
@xseealso{@ref{XREFprimes,,primes}, @ref{XREFisprime,,isprime}}
@end deftypefn


@c sign libinterp/corefcn/mappers.cc
@anchor{XREFsign}
@deftypefn {} {} sign (@var{x})
Compute the @dfn{signum} function.

This is defined as
@tex
$$
{\rm sign} (@var{x}) = \cases{1,&$x>0$;\cr 0,&$x=0$;\cr -1,&$x<0$.\cr}
$$
@end tex
@ifnottex

@example
@group
           -1, x < 0;
sign (x) =  0, x = 0;
            1, x > 0.
@end group
@end example

@end ifnottex

For complex arguments, @code{sign} returns @code{x ./ abs (@var{x})}.

Note that @code{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 @code{signbit} function.
@xseealso{@ref{XREFsignbit,,signbit}}
@end deftypefn


@c signbit libinterp/corefcn/mappers.cc
@anchor{XREFsignbit}
@deftypefn {} {} signbit (@var{x})
Return logical true if the value of @var{x} has its sign bit set and false
otherwise.

This behavior is consistent with the other logical functions.
See @ref{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 @code{x < 0.0}, because IEEE 754 floating point
allows zero to be signed.  The comparison @code{-0.0 < 0.0} is false,
but @code{signbit (-0.0)} will return a nonzero value.
@xseealso{@ref{XREFsign,,sign}}
@end deftypefn


@node Special Functions
@section Special Functions

@c airy libinterp/corefcn/besselj.cc
@anchor{XREFairy}
@deftypefn {} {[@var{a}, @var{ierr}] =} airy (@var{k}, @var{z}, @var{opt})
Compute Airy functions of the first and second kind, and their derivatives.

@example
@group
 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))))
@end group
@end example

The function call @code{airy (@var{z})} is equivalent to
@code{airy (0, @var{z})}.

The result is the same size as @var{z}.

If requested, @var{ierr} contains the following status information and
is the same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

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

@item
Error---no computation, algorithm termination condition not met,
return @code{NaN}.
@end enumerate
@end deftypefn


@c besselj libinterp/corefcn/besselj.cc
@anchor{XREFbesselj}
@deftypefn  {} {@var{J} =} besselj (@var{alpha}, @var{x})
@deftypefnx {} {@var{J} =} besselj (@var{alpha}, @var{x}, @var{opt})
@deftypefnx {} {[@var{J}, @var{ierr}] =} besselj (@dots{})
Compute Bessel functions of the first kind.

The order of the Bessel function @var{alpha} must be real.  The points for
evaluation @var{x} may be complex.

If the optional argument @var{opt} is 1 or true, the result @var{J} is
multiplied by @w{@code{exp (-abs (imag (@var{x})))}}.

If @var{alpha} is a scalar, the result is the same size as @var{x}.  If @var{x}
is a scalar, the result is the same size as @var{alpha}.  If @var{alpha} is a
row vector and @var{x} is a column vector, the result is a matrix with
@code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
size.

If requested, @var{ierr} contains the following status information and is the
same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

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

@item
Error---no computation, algorithm termination condition not met, return
@code{NaN}.
@end enumerate

@xseealso{@ref{XREFbessely,,bessely}, @ref{XREFbesseli,,besseli}, @ref{XREFbesselk,,besselk}, @ref{XREFbesselh,,besselh}}
@end deftypefn


@c bessely libinterp/corefcn/besselj.cc
@anchor{XREFbessely}
@deftypefn  {} {@var{Y} =} bessely (@var{alpha}, @var{x})
@deftypefnx {} {@var{Y} =} bessely (@var{alpha}, @var{x}, @var{opt})
@deftypefnx {} {[@var{Y}, @var{ierr}] =} bessely (@dots{})
Compute Bessel functions of the second kind.

The order of the Bessel function @var{alpha} must be real.  The points for
evaluation @var{x} may be complex.

If the optional argument @var{opt} is 1 or true, the result @var{Y} is
multiplied by @w{@code{exp (-abs (imag (@var{x})))}}.

If @var{alpha} is a scalar, the result is the same size as @var{x}.  If @var{x}
is a scalar, the result is the same size as @var{alpha}.  If @var{alpha} is a
row vector and @var{x} is a column vector, the result is a matrix with
@code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
size.

If requested, @var{ierr} contains the following status information and is the
same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

@item
Complete loss of significance by argument reduction, return @code{NaN}.

@item
Error---no computation, algorithm termination condition not met, return
@code{NaN}.
@end enumerate

@xseealso{@ref{XREFbesselj,,besselj}, @ref{XREFbesseli,,besseli}, @ref{XREFbesselk,,besselk}, @ref{XREFbesselh,,besselh}}
@end deftypefn


@c besseli libinterp/corefcn/besselj.cc
@anchor{XREFbesseli}
@deftypefn  {} {@var{I} =} besseli (@var{alpha}, @var{x})
@deftypefnx {} {@var{I} =} besseli (@var{alpha}, @var{x}, @var{opt})
@deftypefnx {} {[@var{I}, @var{ierr}] =} besseli (@dots{})
Compute modified Bessel functions of the first kind.

The order of the Bessel function @var{alpha} must be real.  The points for
evaluation @var{x} may be complex.

If the optional argument @var{opt} is 1 or true, the result @var{I} is
multiplied by @w{@code{exp (-abs (real (@var{x})))}}.

If @var{alpha} is a scalar, the result is the same size as @var{x}.  If @var{x}
is a scalar, the result is the same size as @var{alpha}.  If @var{alpha} is a
row vector and @var{x} is a column vector, the result is a matrix with
@code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
size.

If requested, @var{ierr} contains the following status information and is the
same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

@item
Complete loss of significance by argument reduction, return @code{NaN}.

@item
Error---no computation, algorithm termination condition not met, return
@code{NaN}.
@end enumerate

@xseealso{@ref{XREFbesselk,,besselk}, @ref{XREFbesselj,,besselj}, @ref{XREFbessely,,bessely}, @ref{XREFbesselh,,besselh}}
@end deftypefn


@c besselk libinterp/corefcn/besselj.cc
@anchor{XREFbesselk}
@deftypefn  {} {@var{K} =} besselk (@var{alpha}, @var{x})
@deftypefnx {} {@var{K} =} besselk (@var{alpha}, @var{x}, @var{opt})
@deftypefnx {} {[@var{K}, @var{ierr}] =} besselk (@dots{})

Compute modified Bessel functions of the second kind.

The order of the Bessel function @var{alpha} must be real.  The points for
evaluation @var{x} may be complex.

If the optional argument @var{opt} is 1 or true, the result @var{K} is
multiplied by @w{@code{exp (@var{x})}}.

If @var{alpha} is a scalar, the result is the same size as @var{x}.  If @var{x}
is a scalar, the result is the same size as @var{alpha}.  If @var{alpha} is a
row vector and @var{x} is a column vector, the result is a matrix with
@code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
size.

If requested, @var{ierr} contains the following status information and is the
same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

@item
Complete loss of significance by argument reduction, return @code{NaN}.

@item
Error---no computation, algorithm termination condition not met, return
@code{NaN}.
@end enumerate

@xseealso{@ref{XREFbesseli,,besseli}, @ref{XREFbesselj,,besselj}, @ref{XREFbessely,,bessely}, @ref{XREFbesselh,,besselh}}
@end deftypefn


@c besselh libinterp/corefcn/besselj.cc
@anchor{XREFbesselh}
@deftypefn  {} {@var{H} =} besselh (@var{alpha}, @var{x})
@deftypefnx {} {@var{H} =} besselh (@var{alpha}, @var{k}, @var{x})
@deftypefnx {} {@var{H} =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})
@deftypefnx {} {[@var{H}, @var{ierr}] =} besselh (@dots{})
Compute Bessel functions of the third kind (Hankel functions).

The order of the Bessel function @var{alpha} must be real.  The kind of Hankel
function is specified by @var{k} and may be either first (@var{k} = 1) or
second (@var{k} = 2).  The default is Hankel functions of the first kind.  The
points for evaluation @var{x} may be complex.

If the optional argument @var{opt} is 1 or true, the result is multiplied
by @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
@var{k} = 2.

If @var{alpha} is a scalar, the result is the same size as @var{x}.  If @var{x}
is a scalar, the result is the same size as @var{alpha}.  If @var{alpha} is a
row vector and @var{x} is a column vector, the result is a matrix with
@code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
size.

If requested, @var{ierr} contains the following status information and is the
same size as the result.

@enumerate 0
@item
Normal return.

@item
Input error, return @code{NaN}.

@item
Overflow, return @code{Inf}.

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

@item
Complete loss of significance by argument reduction, return @code{NaN}.

@item
Error---no computation, algorithm termination condition not met, return
@code{NaN}.
@end enumerate

@xseealso{@ref{XREFbesselj,,besselj}, @ref{XREFbessely,,bessely}, @ref{XREFbesseli,,besseli}, @ref{XREFbesselk,,besselk}}
@end deftypefn


@c beta scripts/specfun/beta.m
@anchor{XREFbeta}
@deftypefn {} {} beta (@var{a}, @var{b})
Compute the Beta function for real inputs @var{a} and @var{b}.

The Beta function definition is
@tex
$$
 B (a, b) = {\Gamma (a) \Gamma (b) \over \Gamma (a + b)}.
$$
@end tex
@ifnottex

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

@end ifnottex

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.
@xref{XREFbetaln,,betaln}, for computing the logarithm of the Beta function
in an efficient manner.
@xseealso{@ref{XREFbetaln,,betaln}, @ref{XREFbetainc,,betainc}, @ref{XREFbetaincinv,,betaincinv}}
@end deftypefn


@c betainc scripts/specfun/betainc.m
@anchor{XREFbetainc}
@deftypefn  {} {} betainc (@var{x}, @var{a}, @var{b})
@deftypefnx {} {} betainc (@var{x}, @var{a}, @var{b}, @var{tail})
Compute the incomplete beta function.

This is defined as
@tex
$$
I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
$$
@end tex
@ifnottex

@example
@group
               x
              /
             |
I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt
             |
             /
            0
@end group
@end example

@end ifnottex

with real @var{x} in the range [0,1].  The inputs @var{a} and @var{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, @var{tail} is @qcode{"lower"} and the incomplete beta function
integrated from 0 to @var{x} is computed.  If @var{tail} is @qcode{"upper"}
then the complementary function integrated from @var{x} to 1 is calculated.
The two choices are related by

betainc (@var{x}, @var{a}, @var{b}, @qcode{"upper"}) =
1 - betainc (@var{x}, @var{a}, @var{b}, @qcode{"lower"}).

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

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

@xseealso{@ref{XREFbeta,,beta}, @ref{XREFbetaincinv,,betaincinv}, @ref{XREFbetaln,,betaln}}
@end deftypefn


@c betaincinv scripts/specfun/betaincinv.m
@anchor{XREFbetaincinv}
@deftypefn  {} {} betaincinv (@var{y}, @var{a}, @var{b})
@deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "lower")
@deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "upper")
Compute the inverse of the normalized incomplete beta function.

The normalized incomplete beta function is defined as
@tex
$$
 I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
$$
@end tex
@ifnottex

@example
@group
               x
              /
             |
I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt
             |
             /
            0
@end group
@end example

@end ifnottex

If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})}
is returned for each of the other inputs.

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

The variable @var{y} must be in the interval [0,1], while @var{a} and
@var{b} must be real and strictly positive.

By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete
beta function integrated from 0 to @var{x} is computed.  If @var{tail} is
@qcode{"upper"} then the complementary function integrated from @var{x} to 1
is inverted.

The function is computed by standard Newton's method, by solving
@tex
$$
 y - I_x (a, b) = 0
$$
@end tex
@ifnottex

@example
@var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0
@end example

@end ifnottex

@xseealso{@ref{XREFbetainc,,betainc}, @ref{XREFbeta,,beta}, @ref{XREFbetaln,,betaln}}
@end deftypefn


@c betaln scripts/specfun/betaln.m
@anchor{XREFbetaln}
@deftypefn {} {} betaln (@var{a}, @var{b})
Compute the natural logarithm of the Beta function for real inputs @var{a}
and @var{b}.

@code{betaln} is defined as
@tex
$$
 {\rm betaln} (a, b) = \ln (B (a,b)) \equiv \ln ({\Gamma (a) \Gamma (b) \over \Gamma (a + b)}).
$$
@end tex
@ifnottex

@example
betaln (a, b) = log (beta (a, b))
@end example

@end ifnottex
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.
@xseealso{@ref{XREFbeta,,beta}, @ref{XREFbetainc,,betainc}, @ref{XREFbetaincinv,,betaincinv}, @ref{XREFgammaln,,gammaln}}
@end deftypefn


@c bincoeff scripts/general/bincoeff.m
@anchor{XREFbincoeff}
@deftypefn {} {} bincoeff (@var{n}, @var{k})
Return the binomial coefficient of @var{n} and @var{k}.

The binomial coefficient is defined as
@tex
$$
 {n \choose k} = {n (n-1) (n-2) \cdots (n-k+1) \over k!}
$$
@end tex
@ifnottex

@example
@group
 /   \
 | n |    n (n-1) (n-2) @dots{} (n-k+1)
 |   |  = -------------------------
 | k |               k!
 \   /
@end group
@end example

@end ifnottex
For example:

@example
@group
bincoeff (5, 2)
   @result{} 10
@end group
@end example

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

@xseealso{@ref{XREFnchoosek,,nchoosek}}
@end deftypefn


@c commutation_matrix scripts/linear-algebra/commutation_matrix.m
@anchor{XREFcommutation_matrix}
@deftypefn {} {} commutation_matrix (@var{m}, @var{n})
Return the commutation matrix
@tex
 $K_{m,n}$
@end tex
@ifnottex
K(m,n)
@end ifnottex
which is the unique
@tex
 $m n \times m n$
@end tex
@ifnottex
@var{m}*@var{n} by @var{m}*@var{n}
@end ifnottex
matrix such that
@tex
 $K_{m,n} \cdot {\rm vec} (A) = {\rm vec} (A^T)$
@end tex
@ifnottex
@math{K(m,n) * vec(A) = vec(A')}
@end ifnottex
for all
@tex
 $m\times n$
@end tex
@ifnottex
@math{m} by @math{n}
@end ifnottex
matrices
@tex
 $A$.
@end tex
@ifnottex
@math{A}.
@end ifnottex

If only one argument @var{m} is given,
@tex
 $K_{m,m}$
@end tex
@ifnottex
@math{K(m,m)}
@end ifnottex
is returned.

See @nospell{Magnus and Neudecker} (1988), @cite{Matrix Differential
Calculus with Applications in Statistics and Econometrics}.
@end deftypefn


@c cosint scripts/specfun/cosint.m
@anchor{XREFcosint}
@deftypefn {} {} cosint (@var{x})
Compute the cosine integral function:
@tex
$$
{\rm Ci} (x) = - \int_x^\infty {{\cos (t)} \over t} dt
$$
@end tex
@ifnottex

@example
@group
            +oo
           /
Ci (x) = - | (cos (t)) / t dt
           /
          x
@end group
@end example

@end ifnottex
An equivalent definition is
@tex
$$
{\rm Ci} (x) = \gamma + \log (x) + \int_0^x {{\cos (t) - 1} \over t} dt
$$
@end tex
@ifnottex

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

@end ifnottex
Reference:

@nospell{M. Abramowitz and I.A. Stegun},
@cite{Handbook of Mathematical Functions}, 1964.

@xseealso{@ref{XREFsinint,,sinint}, @ref{XREFexpint,,expint}, @ref{XREFcos,,cos}}

@end deftypefn


@c duplication_matrix scripts/linear-algebra/duplication_matrix.m
@anchor{XREFduplication_matrix}
@deftypefn {} {} duplication_matrix (@var{n})
Return the duplication matrix
@tex
 $D_n$
@end tex
@ifnottex
@nospell{@math{Dn}}
@end ifnottex
which is the unique
@tex
 $n^2 \times n(n+1)/2$
@end tex
@ifnottex
@math{n^2} by @math{n*(n+1)/2}
@end ifnottex
matrix such that
@tex
 $D_n * {\rm vech} (A) = {\rm vec} (A)$
@end tex
@ifnottex
@nospell{@math{Dn vech (A) = vec (A)}}
@end ifnottex
for all symmetric
@tex
 $n \times n$
@end tex
@ifnottex
@math{n} by @math{n}
@end ifnottex
matrices
@tex
 $A$.
@end tex
@ifnottex
@math{A}.
@end ifnottex

See @nospell{Magnus and Neudecker} (1988), @cite{Matrix Differential
Calculus with Applications in Statistics and Econometrics}.
@end deftypefn


@c dawson libinterp/corefcn/mappers.cc
@anchor{XREFdawson}
@deftypefn {} {} dawson (@var{z})
Compute the Dawson (scaled imaginary error) function.

The Dawson function is defined as
@tex
$$
 {\sqrt{\pi} \over 2} e^{-z^2} {\rm erfi} (z) \equiv -i {\sqrt{\pi} \over 2} e^{-z^2} {\rm erf} (iz)
$$
@end tex
@ifnottex

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

@end ifnottex
@xseealso{@ref{XREFerfc,,erfc}, @ref{XREFerf,,erf}, @ref{XREFerfcx,,erfcx}, @ref{XREFerfi,,erfi}, @ref{XREFerfinv,,erfinv}, @ref{XREFerfcinv,,erfcinv}}
@end deftypefn


@c ellipj libinterp/corefcn/ellipj.cc
@anchor{XREFellipj}
@deftypefn  {} {[@var{sn}, @var{cn}, @var{dn}, @var{err}] =} ellipj (@var{u}, @var{m})
@deftypefnx {} {[@var{sn}, @var{cn}, @var{dn}, @var{err}] =} ellipj (@var{u}, @var{m}, @var{tol})
Compute the Jacobi elliptic functions @var{sn}, @var{cn}, and @var{dn}
of complex argument @var{u} and real parameter @var{m}.

If @var{m} is a scalar, the results are the same size as @var{u}.
If @var{u} is a scalar, the results are the same size as @var{m}.
If @var{u} is a column vector and @var{m} is a row vector, the
results are matrices with @code{length (@var{u})} rows and
@code{length (@var{m})} columns.  Otherwise, @var{u} and
@var{m} must conform in size and the results will be the same size as the
inputs.

The value of @var{u} may be complex.
The value of @var{m} must be 0 @leq{} @var{m} @leq{} 1.

The optional input @var{tol} is currently ignored (@sc{matlab} uses this to
allow faster, less accurate approximation).

If requested, @var{err} contains the following status information
and is the same size as the result.

@enumerate 0
@item
Normal return.

@item
Error---no computation, algorithm termination condition not met,
return @code{NaN}.
@end enumerate

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

@xseealso{@ref{XREFellipke,,ellipke}}
@end deftypefn


@c ellipke scripts/specfun/ellipke.m
@anchor{XREFellipke}
@deftypefn  {} {@var{k} =} ellipke (@var{m})
@deftypefnx {} {@var{k} =} ellipke (@var{m}, @var{tol})
@deftypefnx {} {[@var{k}, @var{e}] =} ellipke (@dots{})
Compute complete elliptic integrals of the first K(@var{m}) and second
E(@var{m}) kind.

@var{m} must be a scalar or real array with -Inf @leq{} @var{m} @leq{} 1.

The optional input @var{tol} controls the stopping tolerance of the
algorithm and defaults to @code{eps (class (@var{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

@tex
$$
{\rm K} (m) = \int_0^1 {dt \over \sqrt{(1 - t^2) (1 - m t^2)}}
$$
@end tex
@ifnottex

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

@end ifnottex

Elliptic integrals of the second kind are defined as

@tex
$$
{\rm E} (m) = \int_0^1 {\sqrt{1 - m t^2} \over \sqrt{1 - t^2}} dt
$$
@end tex
@ifnottex

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

@end ifnottex

Reference: Milton @nospell{Abramowitz} and Irene A. @nospell{Stegun},
@cite{Handbook of Mathematical Functions}, Chapter 17, Dover, 1965.
@xseealso{@ref{XREFellipj,,ellipj}}
@end deftypefn


@c erf libinterp/corefcn/mappers.cc
@anchor{XREFerf}
@deftypefn {} {} erf (@var{z})
Compute the error function.

The error function is defined as
@tex
$$
 {\rm erf} (z) = {2 \over \sqrt{\pi}}\int_0^z e^{-t^2} dt
$$
@end tex
@ifnottex

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

@end ifnottex
@xseealso{@ref{XREFerfc,,erfc}, @ref{XREFerfcx,,erfcx}, @ref{XREFerfi,,erfi}, @ref{XREFdawson,,dawson}, @ref{XREFerfinv,,erfinv}, @ref{XREFerfcinv,,erfcinv}}
@end deftypefn


@c erfc libinterp/corefcn/mappers.cc
@anchor{XREFerfc}
@deftypefn {} {} erfc (@var{z})
Compute the complementary error function.

The complementary error function is defined as
@tex
$1 - {\rm erf} (z)$.
@end tex
@ifnottex
@w{@code{1 - erf (@var{z})}}.
@end ifnottex
@xseealso{@ref{XREFerfcinv,,erfcinv}, @ref{XREFerfcx,,erfcx}, @ref{XREFerfi,,erfi}, @ref{XREFdawson,,dawson}, @ref{XREFerf,,erf}, @ref{XREFerfinv,,erfinv}}
@end deftypefn


@c erfcx libinterp/corefcn/mappers.cc
@anchor{XREFerfcx}
@deftypefn {} {} erfcx (@var{z})
Compute the scaled complementary error function.

The scaled complementary error function is defined as
@tex
$$
 e^{z^2} {\rm erfc} (z) \equiv e^{z^2} (1 - {\rm erf} (z))
$$
@end tex
@ifnottex

@example
exp (z^2) * erfc (z)
@end example

@end ifnottex
@xseealso{@ref{XREFerfc,,erfc}, @ref{XREFerf,,erf}, @ref{XREFerfi,,erfi}, @ref{XREFdawson,,dawson}, @ref{XREFerfinv,,erfinv}, @ref{XREFerfcinv,,erfcinv}}
@end deftypefn


@c erfi libinterp/corefcn/mappers.cc
@anchor{XREFerfi}
@deftypefn {} {} erfi (@var{z})
Compute the imaginary error function.

The imaginary error function is defined as
@tex
$$
 -i {\rm erf} (iz)
$$
@end tex
@ifnottex

@example
-i * erf (i*z)
@end example

@end ifnottex
@xseealso{@ref{XREFerfc,,erfc}, @ref{XREFerf,,erf}, @ref{XREFerfcx,,erfcx}, @ref{XREFdawson,,dawson}, @ref{XREFerfinv,,erfinv}, @ref{XREFerfcinv,,erfcinv}}
@end deftypefn


@c erfinv libinterp/corefcn/mappers.cc
@anchor{XREFerfinv}
@deftypefn {} {} erfinv (@var{x})
Compute the inverse error function.

The inverse error function is defined such that

@example
erf (@var{y}) == @var{x}
@end example
@xseealso{@ref{XREFerf,,erf}, @ref{XREFerfc,,erfc}, @ref{XREFerfcx,,erfcx}, @ref{XREFerfi,,erfi}, @ref{XREFdawson,,dawson}, @ref{XREFerfcinv,,erfcinv}}
@end deftypefn


@c erfcinv libinterp/corefcn/mappers.cc
@anchor{XREFerfcinv}
@deftypefn {} {} erfcinv (@var{x})
Compute the inverse complementary error function.

The inverse complementary error function is defined such that

@example
erfc (@var{y}) == @var{x}
@end example
@xseealso{@ref{XREFerfc,,erfc}, @ref{XREFerf,,erf}, @ref{XREFerfcx,,erfcx}, @ref{XREFerfi,,erfi}, @ref{XREFdawson,,dawson}, @ref{XREFerfinv,,erfinv}}
@end deftypefn


@c expint scripts/specfun/expint.m
@anchor{XREFexpint}
@deftypefn {} {} expint (@var{x})
Compute the exponential integral.

The exponential integral is defined as:

@tex
$$
{\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt
$$
@end tex
@ifnottex

@example
@group
           +oo
          /
          | exp (-t)
E_1 (x) = | -------- dt
          |    t
          /
         x
@end group
@end example

@end ifnottex

Note: For compatibility, this function uses the @sc{matlab} definition
of the exponential integral.  Most other sources refer to this particular
value as @math{E_1 (x)}, and the exponential integral as
@tex
$$
{\rm Ei} (x) = - \int_{-x}^\infty {e^{-t} \over t} dt.
$$
@end tex
@ifnottex

@example
@group
            +oo
           /
           | exp (-t)
Ei (x) = - | -------- dt
           |    t
           /
         -x
@end group
@end example

@end ifnottex
The two definitions are related, for positive real values of @var{x}, by
@tex
$
E_1 (-x) = -{\rm Ei} (x) - i\pi.
$
@end tex
@ifnottex
@w{@code{E_1 (-x) = -Ei (x) - i*pi}}.
@end ifnottex

References:

@nospell{M. Abramowitz and I.A. Stegun},
@cite{Handbook of Mathematical Functions}, 1964.

@nospell{N. Bleistein and R.A. Handelsman},
@cite{Asymptotic expansions of integrals}, 1986.

@xseealso{@ref{XREFcosint,,cosint}, @ref{XREFsinint,,sinint}, @ref{XREFexp,,exp}}
@end deftypefn


@c gamma libinterp/corefcn/mappers.cc
@anchor{XREFgamma}
@deftypefn {} {} gamma (@var{z})
Compute the Gamma function.

The Gamma function is defined as
@tex
$$
 \Gamma (z) = \int_0^\infty t^{z-1} e^{-t} dt.
$$
@end tex
@ifnottex

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

@end ifnottex

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 (@code{gammaln}) in calculations to minimize
loss of precision.  The final result is then
@code{exp (@var{result_using_gammaln}).}
@xseealso{@ref{XREFgammainc,,gammainc}, @ref{XREFgammaln,,gammaln}, @ref{XREFfactorial,,factorial}}
@end deftypefn


@c gammainc scripts/specfun/gammainc.m
@anchor{XREFgammainc}
@deftypefn  {} {} gammainc (@var{x}, @var{a})
@deftypefnx {} {} gammainc (@var{x}, @var{a}, @var{tail})
Compute the normalized incomplete gamma function.

This is defined as
@tex
$$
 \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
$$
@end tex
@ifnottex

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

@end ifnottex
with the limiting value of 1 as @var{x} approaches infinity.
The standard notation is @math{P(a,x)}, e.g., @nospell{Abramowitz} and
@nospell{Stegun} (6.5.1).

If @var{a} is scalar, then @code{gammainc (@var{x}, @var{a})} is returned
for each element of @var{x} and vice versa.

If neither @var{x} nor @var{a} is scalar then the sizes of @var{x} and
@var{a} must agree, and @code{gammainc} is applied element-by-element.
The elements of @var{a} must be non-negative.

By default, @var{tail} is @qcode{"lower"} and the incomplete gamma function
integrated from 0 to @var{x} is computed.  If @var{tail} is @qcode{"upper"}
then the complementary function integrated from @var{x} to infinity is
calculated.

If @var{tail} is @qcode{"scaledlower"}, then the lower incomplete gamma
function is multiplied by
@tex
$\Gamma(a+1)\exp(x)x^{-a}$.
@end tex
@ifnottex
@math{gamma(a+1)*exp(x)/(x^a)}.
@end ifnottex
If @var{tail} is @qcode{"scaledupper"}, then the upper incomplete gamma
function is multiplied by the same quantity.

References:

@nospell{M. Abramowitz and I.A. Stegun},
@cite{Handbook of mathematical functions},
@nospell{Dover publications, Inc.}, 1972.

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

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

@xseealso{@ref{XREFgamma,,gamma}, @ref{XREFgammaincinv,,gammaincinv}, @ref{XREFgammaln,,gammaln}}
@end deftypefn


@c gammaincinv scripts/specfun/gammaincinv.m
@anchor{XREFgammaincinv}
@deftypefn  {} {} gammaincinv (@var{y}, @var{a})
@deftypefnx {} {} gammaincinv (@var{y}, @var{a}, @var{tail})
Compute the inverse of the normalized incomplete gamma function.

The normalized incomplete gamma function is defined as
@tex
$$
 \gamma (x, a) = {1 \over {\Gamma (a)}}\displaystyle{\int_0^x t^{a-1} e^{-t} dt}
$$
@end tex
@ifnottex

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

@end ifnottex

and @code{gammaincinv (gammainc (@var{x}, @var{a}), @var{a}) = @var{x}}
for each non-negative value of @var{x}.  If @var{a} is scalar then
@code{gammaincinv (@var{y}, @var{a})} is returned for each element of
@var{y} and vice versa.

If neither @var{y} nor @var{a} is scalar then the sizes of @var{y} and
@var{a} must agree, and @code{gammaincinv} is applied element-by-element.
The variable @var{y} must be in the interval @math{[0,1]} while @var{a} must
be real and positive.

By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete
gamma function integrated from 0 to @var{x} is computed.  If @var{tail} is
@qcode{"upper"}, then the complementary function integrated from @var{x} to
infinity is inverted.

The function is computed with Newton's method by solving
@tex
$$
 y - \gamma (x, a) = 0
$$
@end tex
@ifnottex

@example
@var{y} - gammainc (@var{x}, @var{a}) = 0
@end example

@end ifnottex

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

@xseealso{@ref{XREFgammainc,,gammainc}, @ref{XREFgamma,,gamma}, @ref{XREFgammaln,,gammaln}}
@end deftypefn


@c legendre scripts/specfun/legendre.m
@anchor{XREFlegendre}
@deftypefn  {} {@var{l} =} legendre (@var{n}, @var{x})
@deftypefnx {} {@var{l} =} legendre (@var{n}, @var{x}, @var{normalization})
Compute the associated Legendre function of degree @var{n} and order
@var{m} = 0 @dots{} @var{n}.

The value @var{n} must be a real non-negative integer.

@var{x} is a vector with real-valued elements in the range [-1, 1].

The optional argument @var{normalization} may be one of @qcode{"unnorm"},
@qcode{"sch"}, or @qcode{"norm"}.  The default if no normalization is given
is @qcode{"unnorm"}.

When the optional argument @var{normalization} is @qcode{"unnorm"}, compute
the associated Legendre function of degree @var{n} and order @var{m} and
return all values for @var{m} = 0 @dots{} @var{n}.  The return value has one
dimension more than @var{x}.

The associated Legendre function of degree @var{n} and order @var{m}:

@tex
$$
P^m_n(x) = (-1)^m (1-x^2)^{m/2}{d^m\over {dx^m}}P_n (x)
$$
@end tex
@ifnottex

@example
@group
 m         m      2  m/2   d^m
P(x) = (-1) * (1-x  )    * ----  P(x)
 n                         dx^m   n
@end group
@end example

@end ifnottex

@noindent
with Legendre polynomial of degree @var{n}:

@tex
$$
P(x) = {1\over{2^n n!}}\biggl({d^n\over{dx^n}}(x^2 - 1)^n\biggr)
$$
@end tex
@ifnottex

@example
@group
          1    d^n   2    n
P(x) = ------ [----(x - 1) ]
 n     2^n n!  dx^n
@end group
@end example

@end ifnottex

@noindent
@code{legendre (3, [-1.0, -0.9, -0.8])} returns the matrix:

@example
@group
 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
@end group
@end example

When the optional argument @var{normalization} is @qcode{"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 @var{n} and order 0:

@tex
$$
SP^0_n (x) = P^0_n (x)
$$
@end tex
@ifnottex

@example
@group
  0      0
SP(x) = P(x)
  n      n
@end group
@end example

@end ifnottex

For Legendre functions of degree n and order m:

@tex
$$
SP^m_n (x) = P^m_n (x)(-1)^m\biggl({2(n-m)!\over{(n+m)!}}\biggl)^{0.5}
$$
@end tex
@ifnottex

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

@end ifnottex

When the optional argument @var{normalization} is @qcode{"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 @var{n} and order @var{m}

@tex
$$
NP^m_n (x) = P^m_n (x)(-1)^m\biggl({(n+0.5)(n-m)!\over{(n+m)!}}\biggl)^{0.5}
$$
@end tex
@ifnottex

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

@end ifnottex
@end deftypefn


@anchor{XREFgammaln}
@c lgamma libinterp/corefcn/mappers.cc
@anchor{XREFlgamma}
@deftypefn  {} {} gammaln (@var{x})
@deftypefnx {} {} lgamma (@var{x})
Return the natural logarithm of the gamma function of @var{x}.
@xseealso{@ref{XREFgamma,,gamma}, @ref{XREFgammainc,,gammainc}}
@end deftypefn


@c psi libinterp/corefcn/psi.cc
@anchor{XREFpsi}
@deftypefn  {} {} psi (@var{z})
@deftypefnx {} {} psi (@var{k}, @var{z})
Compute the psi (polygamma) function.

The polygamma functions are the @var{k}th derivative of the logarithm
of the gamma function.  If unspecified, @var{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:

@tex
$$
\Psi (z) = {d (log (\Gamma (z))) \over dx}
$$
@end tex
@ifnottex

@example
@group
psi (z) = d (log (gamma (z))) / dx
@end group
@end example

@end ifnottex

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

@xseealso{@ref{XREFgamma,,gamma}, @ref{XREFgammainc,,gammainc}, @ref{XREFgammaln,,gammaln}}
@end deftypefn


@c sinint scripts/specfun/sinint.m
@anchor{XREFsinint}
@deftypefn {} {} sinint (@var{x})
Compute the sine integral function:
@tex
$$
{\rm Si} (x) = \int_0^x {\sin (t) \over t} dt
$$
@end tex
@ifnottex

@example
@group
           x
          /
Si (x) =  | sin (t) / t dt
          /
         0
@end group
@end example

@end ifnottex

Reference:
@nospell{M. Abramowitz and I.A. Stegun},
@cite{Handbook of Mathematical Functions}, 1964.

@xseealso{@ref{XREFcosint,,cosint}, @ref{XREFexpint,,expint}, @ref{XREFsin,,sin}}
@end deftypefn


@node Rational Approximations
@section Rational Approximations

@c rat scripts/general/rat.m
@anchor{XREFrat}
@deftypefn  {} {@var{s} =} rat (@var{x})
@deftypefnx {} {@var{s} =} rat (@var{x}, @var{tol})
@deftypefnx {} {[@var{n}, @var{d}] =} rat (@dots{})

Find a rational approximation of @var{x} to within the tolerance defined by
@var{tol}.

If unspecified, the default tolerance is @code{1e-6 * norm (@var{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 @var{x} such
that @code{@var{x} = @var{n} ./ @var{d}}.

For example:

@example
@group
s = rat (pi)
@result{} s = 3 + 1/(7 + 1/16)

[n, d] = rat (pi)
@result{} n =  355
@result{} d =  113

n / d - pi
@result{} 0.00000026676
@end group
@end example

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

@xseealso{@ref{XREFrats,,rats}, @ref{XREFformat,,format}}
@end deftypefn


@c rats libinterp/corefcn/pr-output.cc
@anchor{XREFrats}
@deftypefn  {} {@var{s} =} rats (@var{x})
@deftypefnx {} {@var{s} =} rats (@var{x}, @var{len})
Convert @var{x} into a rational approximation represented as a string.

A rational approximation to a floating point number is a simple fraction
with numerator @var{N} and denominator @var{D} such that
@code{@var{x} = @var{N}/@var{D}}.

The optional second argument defines the maximum length of the string
representing the elements of @var{x}.  By default, @var{len} is 9.

If the length of the smallest possible rational approximation exceeds
@var{len}, an asterisk (*) padded with spaces will be returned instead.

Example conversion from matrix to string, and back again.

@example
@group
r = rats (hilb (4));
x = str2num (r)
@end group
@end example

@xseealso{@ref{XREFrat,,rat}, @ref{XREFformat,,format}}
@end deftypefn


@node Coordinate Transformations
@section Coordinate Transformations

@c cart2pol scripts/general/cart2pol.m
@anchor{XREFcart2pol}
@deftypefn  {} {[@var{theta}, @var{r}] =} cart2pol (@var{x}, @var{y})
@deftypefnx {} {[@var{theta}, @var{r}, @var{z}] =} cart2pol (@var{x}, @var{y}, @var{z})
@deftypefnx {} {[@var{theta}, @var{r}] =} cart2pol (@var{C})
@deftypefnx {} {[@var{theta}, @var{r}, @var{z}] =} cart2pol (@var{C})
@deftypefnx {} {@var{P} =} cart2pol (@dots{})

Transform Cartesian coordinates to polar or cylindrical coordinates.

The inputs @var{x}, @var{y} (, and @var{z}) must be the same shape, or
scalar.  If called with a single matrix argument then each row of @var{C}
represents the Cartesian coordinate (@var{x}, @var{y} (, @var{z})).

@var{theta} describes the angle relative to the positive x-axis.

@var{r} is the distance to the z-axis @w{(0, 0, z)}.

If only a single return argument is requested then return a matrix @var{P}
where each row represents one polar/(cylindrical) coordinate
(@var{theta}, @var{phi} (, @var{z})).
@xseealso{@ref{XREFpol2cart,,pol2cart}, @ref{XREFcart2sph,,cart2sph}, @ref{XREFsph2cart,,sph2cart}}
@end deftypefn


@c pol2cart scripts/general/pol2cart.m
@anchor{XREFpol2cart}
@deftypefn  {} {[@var{x}, @var{y}] =} pol2cart (@var{theta}, @var{r})
@deftypefnx {} {[@var{x}, @var{y}, @var{z}] =} pol2cart (@var{theta}, @var{r}, @var{z})
@deftypefnx {} {[@var{x}, @var{y}] =} pol2cart (@var{P})
@deftypefnx {} {[@var{x}, @var{y}, @var{z}] =} pol2cart (@var{P})
@deftypefnx {} {@var{C} =} pol2cart (@dots{})
Transform polar or cylindrical coordinates to Cartesian coordinates.

The inputs @var{theta}, @var{r}, (and @var{z}) must be the same shape, or
scalar.  If called with a single matrix argument then each row of @var{P}
represents the polar/(cylindrical) coordinate (@var{theta}, @var{r}
(, @var{z})).

@var{theta} describes the angle relative to the positive x-axis.

@var{r} is the distance to the z-axis (0, 0, z).

If only a single return argument is requested then return a matrix @var{C}
where each row represents one Cartesian coordinate
(@var{x}, @var{y} (, @var{z})).
@xseealso{@ref{XREFcart2pol,,cart2pol}, @ref{XREFsph2cart,,sph2cart}, @ref{XREFcart2sph,,cart2sph}}
@end deftypefn


@c cart2sph scripts/general/cart2sph.m
@anchor{XREFcart2sph}
@deftypefn  {} {[@var{theta}, @var{phi}, @var{r}] =} cart2sph (@var{x}, @var{y}, @var{z})
@deftypefnx {} {[@var{theta}, @var{phi}, @var{r}] =} cart2sph (@var{C})
@deftypefnx {} {@var{S} =} cart2sph (@dots{})
Transform Cartesian coordinates to spherical coordinates.

The inputs @var{x}, @var{y}, and @var{z} must be the same shape, or scalar.
If called with a single matrix argument then each row of @var{C} represents
the Cartesian coordinate (@var{x}, @var{y}, @var{z}).

@var{theta} describes the angle relative to the positive x-axis.

@var{phi} is the angle relative to the xy-plane.

@var{r} is the distance to the origin @w{(0, 0, 0)}.

If only a single return argument is requested then return a matrix @var{S}
where each row represents one spherical coordinate
(@var{theta}, @var{phi}, @var{r}).
@xseealso{@ref{XREFsph2cart,,sph2cart}, @ref{XREFcart2pol,,cart2pol}, @ref{XREFpol2cart,,pol2cart}}
@end deftypefn


@c sph2cart scripts/general/sph2cart.m
@anchor{XREFsph2cart}
@deftypefn  {} {[@var{x}, @var{y}, @var{z}] =} sph2cart (@var{theta}, @var{phi}, @var{r})
@deftypefnx {} {[@var{x}, @var{y}, @var{z}] =} sph2cart (@var{S})
@deftypefnx {} {@var{C} =} sph2cart (@dots{})
Transform spherical coordinates to Cartesian coordinates.

The inputs @var{theta}, @var{phi}, and @var{r} must be the same shape, or
scalar.  If called with a single matrix argument then each row of @var{S}
represents the spherical coordinate (@var{theta}, @var{phi}, @var{r}).

@var{theta} describes the angle relative to the positive x-axis.

@var{phi} is the angle relative to the xy-plane.

@var{r} is the distance to the origin @w{(0, 0, 0)}.

If only a single return argument is requested then return a matrix @var{C}
where each row represents one Cartesian coordinate
(@var{x}, @var{y}, @var{z}).
@xseealso{@ref{XREFcart2sph,,cart2sph}, @ref{XREFpol2cart,,pol2cart}, @ref{XREFcart2pol,,cart2pol}}
@end deftypefn


@node Mathematical Constants
@section Mathematical Constants

@c e libinterp/corefcn/data.cc
@anchor{XREFe}
@deftypefn  {} {} e
@deftypefnx {} {} e (@var{n})
@deftypefnx {} {} e (@var{n}, @var{m})
@deftypefnx {} {} e (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} e (@dots{}, @var{class})
Return a scalar, matrix, or N-dimensional array whose elements are all equal
to the base of natural logarithms.

The constant
@tex
$e$ satisfies the equation $\log (e) = 1$.
@end tex
@ifnottex
@samp{e} satisfies the equation @code{log} (e) = 1.
@end ifnottex

When called with no arguments, return a scalar with the value @math{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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFlog,,log}, @ref{XREFexp,,exp}, @ref{XREFpi,,pi}, @ref{XREFI,,I}}
@end deftypefn


@c pi libinterp/corefcn/data.cc
@anchor{XREFpi}
@deftypefn  {} {} pi
@deftypefnx {} {} pi (@var{n})
@deftypefnx {} {} pi (@var{n}, @var{m})
@deftypefnx {} {} pi (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} pi (@dots{}, @var{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
@tex
diameter($\pi$).
@end tex
@ifnottex
diameter.
@end ifnottex

Internally, @code{pi} is computed as @samp{4.0 * atan (1.0)}.

When called with no arguments, return a scalar with the value of
@tex
$\pi$.
@end tex
@ifnottex
pi.
@end ifnottex

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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFe,,e}, @ref{XREFI,,I}}
@end deftypefn


@c I libinterp/corefcn/data.cc
@anchor{XREFI}
@c List other forms of function in documentation index
@findex i
@findex j
@findex J

@deftypefn  {} {} I
@deftypefnx {} {} I (@var{n})
@deftypefnx {} {} I (@var{n}, @var{m})
@deftypefnx {} {} I (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} I (@dots{}, @var{class})
Return a scalar, matrix, or N-dimensional array whose elements are all equal
to the pure imaginary unit, defined as
@tex
$\sqrt{-1}$.
@end tex
@ifnottex
@w{@code{sqrt (-1)}}.
@end ifnottex

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 @math{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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFe,,e}, @ref{XREFpi,,pi}, @ref{XREFlog,,log}, @ref{XREFexp,,exp}}
@end deftypefn


@c Inf libinterp/corefcn/data.cc
@anchor{XREFInf}
@c List other form of function in documentation index
@findex inf

@deftypefn  {} {} Inf
@deftypefnx {} {} Inf (@var{n})
@deftypefnx {} {} Inf (@var{n}, @var{m})
@deftypefnx {} {} Inf (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} Inf (@dots{}, @var{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.

@example
@group
[ 1/0 e^800 ]
@result{} Inf   Inf
@end group
@end example

When called with no arguments, return a scalar with the value @samp{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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFisinf,,isinf}, @ref{XREFNaN,,NaN}}
@end deftypefn


@c NaN libinterp/corefcn/data.cc
@anchor{XREFNaN}
@c List other form of function in documentation index
@findex nan

@deftypefn  {} {} NaN
@deftypefnx {} {} NaN (@var{n})
@deftypefnx {} {} NaN (@var{n}, @var{m})
@deftypefnx {} {} NaN (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} NaN (@dots{}, @var{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
@tex
($\infty - \infty$), zero divided by zero ($0/0$),
@end tex
@ifnottex
(Inf - Inf), zero divided by zero (0/0),
@end ifnottex
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 @code{isnan} function.

When called with no arguments, return a scalar with the value @samp{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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFisnan,,isnan}, @ref{XREFInf,,Inf}}
@end deftypefn


@c eps libinterp/corefcn/data.cc
@anchor{XREFeps}
@deftypefn  {} {} eps
@deftypefnx {} {} eps (@var{x})
@deftypefnx {} {} eps (@var{n}, @var{m})
@deftypefnx {} {} eps (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} eps (@dots{}, @var{class})
Return a scalar, matrix or N-dimensional array whose elements are all eps,
the machine precision.

More precisely, @code{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,
@code{eps} is approximately
@tex
$2.2204\times10^{-16}$ for double precision and $1.1921\times10^{-7}$
@end tex
@ifnottex
2.2204e-16 for double precision and 1.1921e-07
@end ifnottex
for single precision.

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

Given a single argument @var{x}, return the distance between @var{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 @var{class} specifies the return
type and may be either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFrealmax,,realmax}, @ref{XREFrealmin,,realmin}, @ref{XREFintmax,,intmax}, @ref{XREFflintmax,,flintmax}}
@end deftypefn


@c realmax libinterp/corefcn/data.cc
@anchor{XREFrealmax}
@deftypefn  {} {} realmax
@deftypefnx {} {} realmax (@var{n})
@deftypefnx {} {} realmax (@var{n}, @var{m})
@deftypefnx {} {} realmax (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} realmax (@dots{}, @var{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, @code{realmax} is approximately
@tex
$1.7977\times10^{308}$ for double precision and $3.4028\times10^{38}$
@end tex
@ifnottex
1.7977e+308 for double precision and 3.4028e+38
@end ifnottex
for single precision.

When called with no arguments, return a scalar with the value
@code{realmax (@qcode{"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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFrealmin,,realmin}, @ref{XREFintmax,,intmax}, @ref{XREFflintmax,,flintmax}, @ref{XREFeps,,eps}}
@end deftypefn


@c realmin libinterp/corefcn/data.cc
@anchor{XREFrealmin}
@deftypefn  {} {} realmin
@deftypefnx {} {} realmin (@var{n})
@deftypefnx {} {} realmin (@var{n}, @var{m})
@deftypefnx {} {} realmin (@var{n}, @var{m}, @var{k}, @dots{})
@deftypefnx {} {} realmin (@dots{}, @var{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, @code{realmin} is approximately
@tex
$2.2251\times10^{-308}$ for double precision and $1.1755\times10^{-38}$
@end tex
@ifnottex
2.2251e-308 for double precision and 1.1755e-38
@end ifnottex
for single precision.

When called with no arguments, return a scalar with the value
@code{realmin (@qcode{"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 @var{class} specifies the return type and may be
either @qcode{"double"} or @qcode{"single"}.
@xseealso{@ref{XREFrealmax,,realmax}, @ref{XREFintmin,,intmin}, @ref{XREFeps,,eps}}
@end deftypefn