   Copyright (C) 1996, 1997 John W. Eaton.

   Permission is granted to make and distribute verbatim copies of this
manual provided the copyright notice and this permission notice are
preserved on all copies.

   Permission is granted to copy and distribute modified versions of
this manual under the conditions for verbatim copying, provided that
the entire resulting derived work is distributed under the terms of a
permission notice identical to this one.

   Permission is granted to copy and distribute translations of this
manual into another language, under the above conditions for modified
versions.

Specialized Two-Dimensional Plots
=================================

 - Function File:  bar (X, Y)
     Given two vectors of x-y data, `bar' produces a bar graph.

     If only one argument is given, it is taken as a vector of y-values
     and the x coordinates are taken to be the indices of the elements.

     If two output arguments are specified, the data are generated but
     not plotted.  For example,

          bar (x, y);

     and

          [xb, yb] = bar (x, y);
          plot (xb, yb);

     are equivalent.

 - Function File:  contour (Z, N)
 - Function File:  contour (X, Y, Z, N)
     Make a contour plot of the three-dimensional surface described by
     Z.  Someone needs to improve `gnuplot''s contour routines before
     this will be very useful.

 - Function File:  hist (Y, X, NORM)
     Produce histogram counts or plots.

     With one vector input argument, plot a histogram of the values with
     10 bins.  The range of the histogram bins is determined by the
     range of the data.

     Given a second scalar argument, use that as the number of bins.

     Given a second vector argument, use that as the centers of the
     bins, with the width of the bins determined from the adjacent
     values in the vector.

     If third argument is provided, the histogram is normalised such
     that the sum of the bars is equal to NORM.

     Extreme values are lumped in the first and last bins.

     With two output arguments, produce the values NN and XX such that
     `bar (XX, NN)' will plot the histogram.

 - Function File:  loglog (ARGS)
     Make a two-dimensional plot using log scales for both axes.  See
     the description of `plot' for a description of the arguments that
     `loglog' will accept.

 - Function File:  polar (THETA, RHO, FMT)
     Make a two-dimensional plot given polar the coordinates THETA and
     RHO.

     The optional third argument specifies the line type.

 - Function File:  semilogx (ARGS)
     Make a two-dimensional plot using a log scale for the X axis.  See
     the description of `plot' for a description of the arguments that
     `semilogx' will accept.

 - Function File:  semilogy (ARGS)
     Make a two-dimensional plot using a log scale for the Y axis.  See
     the description of `plot' for a description of the arguments that
     `semilogy' will accept.

 - Function File:  stairs (X, Y)
     Given two vectors of x-y data, bar produces a `stairstep' plot.

     If only one argument is given, it is taken as a vector of y-values
     and the x coordinates are taken to be the indices of the elements.

     If two output arguments are specified, the data are generated but
     not plotted.  For example,

          stairs (x, y);

     and

          [xs, ys] = stairs (x, y);
          plot (xs, ys);

     are equivalent.

 - Function File:  errorbar (ARGS)
     This function produces two-dimensional plots with errorbars. Many
     different combinations of arguments are possible.  The simplest
     form is

          errorbar (Y, EY)

     where the first argument is taken as the set of Y coordinates and
     the second argument EY is taken as the errors of the Y values. X
     coordinates are taken to be the indices of the elements, starting
     with 1.

     If more than two arguments are given, they are interpreted as

          errorbar (X, Y, ..., FMT ...)

     where after X and Y there can be up to four error parameters such
     as EY, EX, LY, UY etc., depending on the plot type. Any number of
     argument sets may appear, as long as they are separated with a
     format string FMT.

     If Y is a matrix, X and error parameters must also be matrices
     having same dimensions. The columns of Y are plotted versus the
     corresponding columns of X and errorbars are drawn from the
     corresponding columns of error parameters.

     If FMT is missing, yerrorbars ("~") plot style is assumed.  If the
     FMT argument is supplied, it is interpreted as in normal plots
     (See __pltopt__). In addition the following plot styles are
     supported by errorbar:

    `~'
          Set yerrorbars plot style (default).

    `>'
          Set xerrorbars plot style.

    `~>'
          Set xyerrorbars plot style.

    `#'
          Set boxes plot style.

    `#~'
          Set boxerrorbars plot style.

    `#~>'
          Set boxxyerrorbars plot style.

     Examples:

          errorbar(X, Y, EX, ">")

     xerrorbar plot of Y versus X with X errorbars drawn from X-EX to
     X+EX.

          errorbar(X, Y1, EY, "~", X, Y2, LY, UY)

     Two yerrorbar plots with Y1 and Y2 versus X.  Errorbars for Y1 are
     drawn from Y1-EY to Y1+EY, errorbars for Y2 from Y2-LY to Y2+UY.

          errorbar(X, Y, LX, UX, LY, UY, "~>")

     xyerrorbar plot of Y versus X where X errorbars are drawn from
     X-LX to X+UX and Y errorbars from Y-LY to Y+UY.

 - Function File:  loglogerr (ARGS)
     This function produces two-dimensional plots on double logarithm
     axis with errorbars. Many different combinations of arguments are
     possible.  The most used form is

          loglogerr (X, Y, EY, FMT)

     which produces a double logarithm plot of Y versus X with errors
     in the Y-scale defined by EY and the plot format defined by FMT.
     See errorbar for available formats and additional information.

 - Function File:  semilogxerr (ARGS)
     This function produces two-dimensional plots on a semilogarithm
     axis with errorbars. Many different combinations of arguments are
     possible.  The most used form is

          semilogxerr (X, Y, EY, FMT)

     which produces a semi-logarithm plot of Y versus X with errors in
     the Y-scale defined by EY and the plot format defined by FMT. See
     errorbar for available formats and additional information.

 - Function File:  semilogyerr (ARGS)
     This function produces two-dimensional plots on a semilogarithm
     axis with errorbars. Many different combinations of arguments are
     possible.  The most used form is

          semilogyerr (X, Y, EY, FMT)

     which produces a semi-logarithm plot of Y versus X with errors in
     the Y-scale defined by EY and the plot format defined by FMT. See
     errorbar for available formats and additional information.

Three-Dimensional Plotting
==========================

 - Command: gsplot RANGES EXPRESSION USING TITLE STYLE
     Generate a 3-dimensional plot.

     The RANGES, USING, TITLE, and STYLE arguments are optional, and
     the USING, TITLE and STYLE qualifiers may appear in any order
     after the expression.  You may plot multiple expressions with a
     single command by separating them with commas.  Each expression
     may have its own set of qualifiers.

     The optional item RANGES has the syntax

          [ x_lo : x_up ] [ y_lo : y_up ] [ z_lo : z_up ]

     and may be used to specify the ranges for the axes of the plot,
     independent of the actual range of the data.  The range for the y
     and z axes and any of the individual limits may be omitted.  A
     range `[:]' indicates that the default limits should be used.  This
     normally means that a range just large enough to include all the
     data points will be used.

     The expression to be plotted must not contain any literal matrices
     (e.g.  `[ 1, 2; 3, 4 ]') since it is nearly impossible to
     distinguish a plot range from a matrix of data.

     See the help for `gnuplot' for a description of the syntax for the
     optional items.

     By default, the `gsplot' command plots each column of the
     expression as the z value, using the row index as the x value, and
     the column index as the y value.  The indices are counted from
     zero, not one.  For example,

          gsplot rand (5, 2)

     will plot a random surface, with the x and y values taken from the
     row and column indices of the matrix.

     If parametric plotting mode is set (using the command `gset
     parametric', then `gsplot' takes the columns of the matrix three
     at a time as the x, y and z values that define a line in three
     space.  Any extra columns are ignored, and the x and y values are
     expected to be sorted.  For example, with `parametric' set, it
     makes sense to plot a matrix like

          1 1 3 2 1 6 3 1 9
          1 2 2 2 2 5 3 2 8
          1 3 1 2 3 4 3 3 7

     but not `rand (5, 30)'.

   The MATLAB-style three-dimensional plotting commands are:

 - Function File:  mesh (X, Y, Z)
     Plot a mesh given matrices X, and Y from `meshdom' and a matrix Z
     corresponding to the X and Y coordinates of the mesh.  If X and Y
     are vectors, then a typical vertex is (X(j), Y(i), Z(i,j)).  Thus,
     columns of Z correspond to different X values and rows of Z
     correspond to different Y values.

 - Function File: [XX, YY] = meshgrid (X, Y)
 - Function File: [XX, YY] = meshgrid (X)
     Given vectors of X and Y coordinates, return two matrices
     corresponding to the X and Y coordinates of a mesh.  The rows of
     XX are copies of X, and the columns of YY are copies of Y.

 - Function File:  meshdom (X, Y)
     Given vectors of X and Y coordinates, return two matrices
     corresponding to the X and Y coordinates of the mesh.

     Note: this function is provided for compatibility with older
     versions of MATLAB.  You should use `meshgrid' instead.

Plot Annotations
================

 - Function File:  grid (ARG)
     For two-dimensional plotting, force the display of a grid on the
     plot.  The argument may be either `"on"' or `"off"'.  If it is
     omitted, `"on"' is assumed.

 - Function File:  title (STRING)
     Specify a title for a plot.  If you already have a plot displayed,
     use the command `replot' to redisplay it with the new title.

 - Function File:  bottom_title (STRING)
     See top_title.

 - Function File:  xlabel (STRING)
 - Function File:  ylabel (STRING)
 - Function File:  zlabel (STRING)
     Specify x, y, and z axis labels for the plot.  If you already have
     a plot displayed, use the command `replot' to redisplay it with
     the new labels.

Multiple Plots on One Page
==========================

   The following functions all require a version of `gnuplot' that
supports the multiplot feature.

 - Function File:  mplot (X, Y)
 - Function File:  mplot (X, Y, FMT)
 - Function File:  mplot (X1, Y1, X2, Y2)
     This is a modified version of the `plot' function that works with
     the multiplot version of `gnuplot' to plot multiple plots per page.
     This plot version automatically advances to the next subplot
     position after each set of arguments are processed.

     See the description of the PLOT function for the various options.

 - Function File:  multiplot (XN, YN)
     Sets and resets multiplot mode.

     If the arguments are non-zero, `multiplot' will set up multiplot
     mode with XN, YN subplots along the X and Y axes.  If both
     arguments are zero, `multiplot' closes multiplot mode.

 - Function File:  oneplot ()
     If in multiplot mode, switches to single plot mode.

 - Function File:  plot_border (...)
     Multiple arguments allowed to specify the sides on which the border
     is shown.  Allowed arguments include:

    `"blank"'
          No borders displayed.

    `"all"'
          All borders displayed

    `"north"'
          North Border

    `"south"'
          South Border

    `"east"'
          East Border

    `"west"'
          West Border

     The arguments may be abbreviated to single characters.  Without any
     arguments, `plot_border' turns borders off.

 - Function File:  subplot (ROWS, COLS, INDEX)
 - Function File:  subplot (RCN)
     Sets `gnuplot' in multiplot mode and plots in location given by
     index (there are COLS by ROWS subwindows).

     Input:

    ROWS
          Number of rows in subplot grid.

    COLUMNS
          Number of columns in subplot grid.

    INDEX
          Index of subplot where to make the next plot.

     If only one argument is supplied, then it must be a three digit
     value specifying the location in digits 1 (rows) and 2 (columns)
     and the plot index in digit 3.

     The plot index runs row-wise.  First all the columns in a row are
     filled and then the next row is filled.

     For example, a plot with 4 by 2 grid will have plot indices
     running as follows:
          +-----+-----+-----+-----+
          |  1  |  2  |  3  |  4  |
          +-----+-----+-----+-----+
          |  5  |  6  |  7  |  8  |
          +-----+-----+-----+-----+

 - Function File:  subwindow (XN, YN)
     Sets the subwindow position in multiplot mode for the next plot.
     The multiplot mode has to be previously initialized using the
     `multiplot' function, otherwise this command just becomes an alias
     to `multiplot'

 - Function File:  top_title (STRING)
 - Function File:  bottom_title (STRING)
     Makes a title with text STRING at the top (bottom) of the plot.

Multiple Plot Windows
=====================

 - Function File:  figure (N)
     Set the current plot window to plot window N.  This function
     currently requires X11 and a version of gnuplot that supports
     multiple frames.

Interaction with `gnuplot'
==========================

 - Built-in Variable: gnuplot_binary
     The name of the program invoked by the plot command.  The default
     value is `"gnuplot"'.  See Installation.

 - Built-in Variable: gnuplot_has_frames
     If the value of this variable is nonzero, Octave assumes that your
     copy of gnuplot has support for multiple frames that is included
     in recent 3.6beta releases.  It's initial value is determined by
     configure, but it can be changed in your startup script or at the
     command line in case configure got it wrong, or if you upgrade
     your gnuplot installation.

 - Built-in Function:  graw (STRING)
     Send STRING directly to gnuplot subprocess.

 - Built-in Variable: gnuplot_command_plot

 - Built-in Variable: gnuplot_command_replot

 - Built-in Variable: gnuplot_command_splot

 - Built-in Variable: gnuplot_command_using

 - Built-in Variable: gnuplot_command_with

 - Built-in Variable: gnuplot_command_axes

 - Built-in Variable: gnuplot_command_title

 - Built-in Variable: gnuplot_command_end

Matrix Manipulation
*******************

   There are a number of functions available for checking to see if the
elements of a matrix meet some condition, and for rearranging the
elements of a matrix.  For example, Octave can easily tell you if all
the elements of a matrix are finite, or are less than some specified
value.  Octave can also rotate the elements, extract the upper- or
lower-triangular parts, or sort the columns of a matrix.

Finding Elements and Checking Conditions
========================================

   The functions `any' and `all' are useful for determining whether any
or all of the elements of a matrix satisfy some condition.  The `find'
function is also useful in determining which elements of a matrix meet
a specified condition.

 - Built-in Function:  any (X, DIM)
     For a vector argument, return 1 if any element of the vector is
     nonzero.

     For a matrix argument, return a row vector of ones and zeros with
     each element indicating whether any of the elements of the
     corresponding column of the matrix are nonzero.  For example,

          any (eye (2, 4))
               => [ 1, 1, 0, 0 ]

     If the optional argument DIM is supplied, work along dimension
     DIM.  For example,

          any (eye (2, 4), 2)
               => [ 1; 1 ]

 - Built-in Function:  all (X, DIM)
     The function `all' behaves like the function `any', except that it
     returns true only if all the elements of a vector, or all the
     elements along dimension DIM of a matrix, are nonzero.

   Since the comparison operators (see Comparison Ops) return
matrices of ones and zeros, it is easy to test a matrix for many
things, not just whether the elements are nonzero.  For example,

     all (all (rand (5) < 0.9))
          => 0

tests a random 5 by 5 matrix to see if all of its elements are less
than 0.9.

   Note that in conditional contexts (like the test clause of `if' and
`while' statements) Octave treats the test as if you had typed `all
(all (condition))'.

 - Mapping Function:  xor (X, Y)
     Return the `exclusive or' of the entries of X and Y.  For boolean
     expressions X and Y, `xor (X, Y)' is true if and only if X or Y is
     true, but not if both X and Y are true.

 - Function File:  is_duplicate_entry (X)
     Return non-zero if any entries in X are duplicates of one another.

 - Function File:  diff (X, K)
     If X is a vector of length N, `diff (X)' is the vector of first
     differences X(2) - X(1), ..., X(n) - X(n-1).

     If X is a matrix, `diff (X)' is the matrix of column differences.

     The second argument is optional.  If supplied, `diff (X, K)',
     where K is a nonnegative integer, returns the K-th differences.

 - Mapping Function:  isinf (X)
     Return 1 for elements of X that are infinite and zero otherwise.
     For example,

          isinf ([13, Inf, NA, NaN])
               => [ 0, 1, 0, 0 ]

 - Mapping Function:  isnan (X)
     Return 1 for elements of X that are NaN values and zero otherwise.
     For example,

          isnan ([13, Inf, NA, NaN])
               => [ 0, 0, 0, 1 ]

 - Mapping Function:  finite (X)
     Return 1 for elements of X that are finite values and zero
     otherwise. For example,

          finite ([13, Inf, NaN])
               => [ 1, 0, 0 ]

 - Loadable Function:  find (X)
     Return a vector of indices of nonzero elements of a matrix.  To
     obtain a single index for each matrix element, Octave pretends
     that the columns of a matrix form one long vector (like Fortran
     arrays are stored).  For example,

          find (eye (2))
               => [ 1; 4 ]

     If two outputs are requested, `find' returns the row and column
     indices of nonzero elements of a matrix.  For example,

          [i, j] = find (2 * eye (2))
               => i = [ 1; 2 ]
               => j = [ 1; 2 ]

     If three outputs are requested, `find' also returns a vector
     containing the nonzero values.  For example,

          [i, j, v] = find (3 * eye (2))
               => i = [ 1; 2 ]
               => j = [ 1; 2 ]
               => v = [ 3; 3 ]

 - Function File: [ERR, Y1, ...] = common_size (X1, ...)
     Determine if all input arguments are either scalar or of common
     size.  If so, ERR is zero, and YI is a matrix of the common size
     with all entries equal to XI if this is a scalar or XI otherwise.
     If the inputs cannot be brought to a common size, errorcode is 1,
     and YI is XI.  For example,

          [errorcode, a, b] = common_size ([1 2; 3 4], 5)
          => errorcode = 0
          => a = [ 1, 2; 3, 4 ]
          => b = [ 5, 5; 5, 5 ]

     This is useful for implementing functions where arguments can
     either be scalars or of common size.

Rearranging Matrices
====================

 - Function File:  fliplr (X)
     Return a copy of X with the order of the columns reversed.  For
     example,

          fliplr ([1, 2; 3, 4])
          =>  2  1
                   4  3

 - Function File:  flipud (X)
     Return a copy of X with the order of the rows reversed.  For
     example,

          flipud ([1, 2; 3, 4])
          =>  3  4
                   1  2

 - Function File:  rot90 (X, N)
     Return a copy of X with the elements rotated counterclockwise in
     90-degree increments.  The second argument is optional, and
     specifies how many 90-degree rotations are to be applied (the
     default value is 1).  Negative values of N rotate the matrix in a
     clockwise direction.  For example,

          rot90 ([1, 2; 3, 4], -1)
          =>  3  1
                   4  2

     rotates the given matrix clockwise by 90 degrees.  The following
     are all equivalent statements:

          rot90 ([1, 2; 3, 4], -1)
          ==
          rot90 ([1, 2; 3, 4], 3)
          ==
          rot90 ([1, 2; 3, 4], 7)

 - Function File:  reshape (A, M, N)
     Return a matrix with M rows and N columns whose elements are taken
     from the matrix A.  To decide how to order the elements, Octave
     pretends that the elements of a matrix are stored in column-major
     order (like Fortran arrays are stored).

     For example,

          reshape ([1, 2, 3, 4], 2, 2)
          =>  1  3
                   2  4

     If the variable `do_fortran_indexing' is nonzero, the `reshape'
     function is equivalent to

          retval = zeros (m, n);
          retval (:) = a;

     but it is somewhat less cryptic to use `reshape' instead of the
     colon operator.  Note that the total number of elements in the
     original matrix must match the total number of elements in the new
     matrix.

 - Function File:  shift (X, B)
     If X is a vector, perform a circular shift of length B of the
     elements of X.

     If X is a matrix, do the same for each column of X.

 - Loadable Function: [S, I] = sort (X)
     Return a copy of X with the elements elements arranged in
     increasing order.  For matrices, `sort' orders the elements in each
     column.

     For example,

          sort ([1, 2; 2, 3; 3, 1])
               =>  1  1
                   2  2
                   3  3

     The `sort' function may also be used to produce a matrix
     containing the original row indices of the elements in the sorted
     matrix.  For example,

          [s, i] = sort ([1, 2; 2, 3; 3, 1])
               => s = 1  1
                      2  2
                      3  3
               => i = 1  3
                      2  1
                      3  2

   Since the `sort' function does not allow sort keys to be specified,
it can't be used to order the rows of a matrix according to the values
of the elements in various columns(1) in a single call.  Using the
second output, however, it is possible to sort all rows based on the
values in a given column.  Here's an example that sorts the rows of a
matrix based on the values in the second column.

     a = [1, 2; 2, 3; 3, 1];
     [s, i] = sort (a (:, 2));
     a (i, :)
          =>  3  1
              1  2
              2  3

 - Function File:  tril (A, K)
 - Function File:  triu (A, K)
     Return a new matrix formed by extracting extract the lower (`tril')
     or upper (`triu') triangular part of the matrix A, and setting all
     other elements to zero.  The second argument is optional, and
     specifies how many diagonals above or below the main diagonal
     should also be set to zero.

     The default value of K is zero, so that `triu' and `tril' normally
     include the main diagonal as part of the result matrix.

     If the value of K is negative, additional elements above (for
     `tril') or below (for `triu') the main diagonal are also selected.

     The absolute value of K must not be greater than the number of
     sub- or super-diagonals.

     For example,

          tril (ones (3), -1)
          =>  0  0  0
                   1  0  0
                   1  1  0

     and

          tril (ones (3), 1)
          =>  1  1  0
                   1  1  1
                   1  1  1

 - Function File:  vec (X)
     Return the vector obtained by stacking the columns of the matrix X
     one above the other.

 - Function File:  vech (X)
     Return the vector obtained by eliminating all supradiagonal
     elements of the square matrix X and stacking the result one column
     above the other.

 - Function File:  prepad (X, L, C)
 - Function File:  postpad (X, L, C)
     Prepends (appends) the scalar value C to the vector X until it is
     of length L.  If the third argument is not supplied, a value of 0
     is used.

     If `length (X) > L', elements from the beginning (end) of X are
     removed until a vector of length L is obtained.

     If X is a matrix, elements are prepended or removed from each row.

   ---------- Footnotes ----------

   (1) For example, to first sort based on the values in column 1, and
then, for any values that are repeated in column 1, sort based on the
values found in column 2, etc.

Special Utility Matrices
========================

 - Built-in Function:  eye (X)
 - Built-in Function:  eye (N, M)
     Return an identity matrix.  If invoked with a single scalar
     argument, `eye' returns a square matrix with the dimension
     specified.  If you supply two scalar arguments, `eye' takes them
     to be the number of rows and columns.  If given a vector with two
     elements, `eye' uses the values of the elements as the number of
     rows and columns, respectively.  For example,

          eye (3)
               =>  1  0  0
                   0  1  0
                   0  0  1

     The following expressions all produce the same result:

          eye (2)
          ==
          eye (2, 2)
          ==
          eye (size ([1, 2; 3, 4])

     For compatibility with MATLAB, calling `eye' with no arguments is
     equivalent to calling it with an argument of 1.

 - Built-in Function:  ones (X)
 - Built-in Function:  ones (N, M)
     Return a matrix whose elements are all 1.  The arguments are
     handled the same as the arguments for `eye'.

     If you need to create a matrix whose values are all the same, you
     should use an expression like

          val_matrix = val * ones (n, m)

 - Built-in Function:  zeros (X)
 - Built-in Function:  zeros (N, M)
     Return a matrix whose elements are all 0.  The arguments are
     handled the same as the arguments for `eye'.

 - Function File:  repmat (A, M, N)
 - Function File:  repmat (A, [M N])
     Form a block matrix of size M by N, with a copy of matrix A as
     each element.  If N is not specified, form an M by M block matrix.

 - Loadable Function:  rand (X)
 - Loadable Function:  rand (N, M)
 - Loadable Function:  rand (`"seed"', X)
     Return a matrix with random elements uniformly distributed on the
     interval (0, 1).  The arguments are handled the same as the
     arguments for `eye'.  In addition, you can set the seed for the
     random number generator using the form

          rand ("seed", X)

     where X is a scalar value.  If called as

          rand ("seed")

     `rand' returns the current value of the seed.

 - Loadable Function:  randn (X)
 - Loadable Function:  randn (N, M)
 - Loadable Function:  randn (`"seed"', X)
     Return a matrix with normally distributed random elements.  The
     arguments are handled the same as the arguments for `eye'.  In
     addition, you can set the seed for the random number generator
     using the form

          randn ("seed", X)

     where X is a scalar value.  If called as

          randn ("seed")

     `randn' returns the current value of the seed.

   The `rand' and `randn' functions use separate generators.  This
ensures that

     rand ("seed", 13);
     randn ("seed", 13);
     u = rand (100, 1);
     n = randn (100, 1);

and

     rand ("seed", 13);
     randn ("seed", 13);
     u = zeros (100, 1);
     n = zeros (100, 1);
     for i = 1:100
       u(i) = rand ();
       n(i) = randn ();
     end

produce equivalent results.

   Normally, `rand' and `randn' obtain their initial seeds from the
system clock, so that the sequence of random numbers is not the same
each time you run Octave.  If you really do need for to reproduce a
sequence of numbers exactly, you can set the seed to a specific value.

   If it is invoked without arguments, `rand' and `randn' return a
single element of a random sequence.

   The `rand' and `randn' functions use Fortran code from RANLIB, a
library of fortran routines for random number generation, compiled by
Barry W. Brown and James Lovato of the Department of Biomathematics at
The University of Texas, M.D. Anderson Cancer Center, Houston, TX 77030.

 - Function File:  randperm (N)
     Return a row vector containing a random permutation of the
     integers from 1 to N.

 - Built-in Function:  diag (V, K)
     Return a diagonal matrix with vector V on diagonal K.  The second
     argument is optional.  If it is positive, the vector is placed on
     the K-th super-diagonal.  If it is negative, it is placed on the
     -K-th sub-diagonal.  The default value of K is 0, and the vector
     is placed on the main diagonal.  For example,

          diag ([1, 2, 3], 1)
               =>  0  1  0  0
                   0  0  2  0
                   0  0  0  3
                   0  0  0  0

   The functions `linspace' and `logspace' make it very easy to create
vectors with evenly or logarithmically spaced elements.  See Ranges.

 - Built-in Function:  linspace (BASE, LIMIT, N)
     Return a row vector with N linearly spaced elements between BASE
     and LIMIT.  The number of elements, N, must be greater than 1.
     The BASE and LIMIT are always included in the range.  If BASE is
     greater than LIMIT, the elements are stored in decreasing order.
     If the number of points is not specified, a value of 100 is used.

     The `linspace' function always returns a row vector, regardless of
     the value of `prefer_column_vectors'.

 - Function File:  logspace (BASE, LIMIT, N)
     Similar to `linspace' except that the values are logarithmically
     spaced from 10^base to 10^limit.

     If LIMIT is equal to pi, the points are between 10^base and pi,
     _not_ 10^base and 10^pi, in order to  be compatible with the
     corresponding MATLAB function.

 - Built-in Variable: treat_neg_dim_as_zero
     If the value of `treat_neg_dim_as_zero' is nonzero, expressions
     like

          eye (-1)

     produce an empty matrix (i.e., row and column dimensions are zero).
     Otherwise, an error message is printed and control is returned to
     the top level.  The default value is 0.

 - Built-in Variable: ok_to_lose_imaginary_part
     If the value of `ok_to_lose_imaginary_part' is nonzero, implicit
     conversions of complex numbers to real numbers are allowed (for
     example, by fsolve).  If the value is `"warn"', the conversion is
     allowed, but a warning is printed.  Otherwise, an error message is
     printed and control is returned to the top level.  The default
     value is `"warn"'.

Famous Matrices
===============

   The following functions return famous matrix forms.

 - Function File:  hankel (C, R)
     Return the Hankel matrix constructed given the first column C, and
     (optionally) the last row R.  If the last element of C is not the
     same as the first element of R, the last element of C is used.  If
     the second argument is omitted, the last row is taken to be the
     same as the first column.

     A Hankel matrix formed from an m-vector C, and an n-vector R, has
     the elements

          H (i, j) = c (i+j-1),  i+j-1 <= m;
          H (i, j) = r (i+j-m),  otherwise

 - Function File:  hilb (N)
     Return the Hilbert matrix of order N.  The i, j element of a
     Hilbert matrix is defined as

          H (i, j) = 1 / (i + j - 1)

 - Function File:  invhilb (N)
     Return the inverse of a Hilbert matrix of order N.  This can be
     computed computed exactly using

                      (i+j)         /n+i-1\  /n+j-1\   /i+j-2\ 2
           A(i,j) = -1      (i+j-1)(       )(       ) (       )
                                    \ n-j /  \ n-i /   \ i-2 /

                  = p(i) p(j) / (i+j-1)
     where
                       k  /k+n-1\   /n\
              p(k) = -1  (       ) (   )
                          \ k-1 /   \k/

     The validity of this formula can easily be checked by expanding
     the binomial coefficients in both formulas as factorials.  It can
     be derived more directly via the theory of Cauchy matrices: see J.
     W. Demmel, Applied Numerical Linear Algebra, page 92.

     Compare this with the numerical calculation of `inverse (hilb
     (n))', which suffers from the ill-conditioning of the Hilbert
     matrix, and the finite precision of your computer's floating point
     arithmetic.

 - Function File:  sylvester_matrix (K)
     Return the Sylvester matrix of order n = 2^k.

 - Function File:  toeplitz (C, R)
     Return the Toeplitz matrix constructed given the first column C,
     and (optionally) the first row R.  If the first element of C is
     not the same as the first element of R, the first element of C is
     used.  If the second argument is omitted, the first row is taken
     to be the same as the first column.

     A square Toeplitz matrix has the form

          c(0)  r(1)   r(2)  ...  r(n)
          c(1)  c(0)   r(1)      r(n-1)
          c(2)  c(1)   c(0)      r(n-2)
           .                       .
           .                       .
           .                       .

          c(n) c(n-1) c(n-2) ...  c(0)

 - Function File:  vander (C)
     Return the Vandermonde matrix whose next to last column is C.

     A Vandermonde matrix has the form

          c(0)^n ... c(0)^2  c(0)  1
          c(1)^n ... c(1)^2  c(1)  1
           .           .      .    .
           .           .      .    .
           .           .      .    .

          c(n)^n ... c(n)^2  c(n)  1

Arithmetic
**********

   Unless otherwise noted, all of the functions described in this
chapter will work for real and complex scalar or matrix arguments.

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

   The following functions are available for working with complex
numbers.  Each expects a single argument.  They are called "mapping
functions" because when given a matrix argument, they apply the given
function to each element of the matrix.

 - Mapping Function:  ceil (X)
     Return the smallest integer not less than X.  If X is complex,
     return `ceil (real (X)) + ceil (imag (X)) * I'.

 - Mapping Function:  exp (X)
     Compute the exponential of X.  To compute the matrix exponential,
     see Linear Algebra.

 - Mapping Function:  fix (X)
     Truncate X toward zero.  If X is complex, return `fix (real (X)) +
     fix (imag (X)) * I'.

 - Mapping Function:  floor (X)
     Return the largest integer not greater than X.  If X is complex,
     return `floor (real (X)) + floor (imag (X)) * I'.

 - Mapping Function:  gcd (X, `...')
     Compute the greatest common divisor of the elements of X, or the
     list of all the arguments.  For example,

          gcd (a1, ..., ak)

     is the same as

          gcd ([a1, ..., ak])

     An optional second return value, V contains an integer vector such
     that

          g = v(1) * a(k) + ... + v(k) * a(k)

 - Mapping Function:  lcm (X, `...')
     Compute the least common multiple of the elements elements of X, or
     the list of all the arguments.  For example,

          lcm (a1, ..., ak)

     is the same as

          lcm ([a1, ..., ak]).

 - Mapping Function:  log (X)
     Compute the natural logarithm for each element of X.  To compute
     the matrix logarithm, see Linear Algebra.

 - Mapping Function:  log10 (X)
     Compute the base-10 logarithm for each element of X.

 - Mapping Function:  log2 (X)
 - Mapping Function: [F, E] log2 (X)
     Compute the base-2 logarithm of X.  With two outputs, returns F
     and E such that  1/2 <= abs(f) < 1 and x = f * 2^e.

   For a vector argument, return the maximum value.  For a matrix
argument, return the maximum value from each column, as a row vector.
Thus,

     max (max (X))

returns the largest element of X.

   For complex arguments, the magnitude of the elements are used for
comparison.

   If called with two output arguments, also returns the index of the
maximum value(s). Thus,
     [x, ix] = max([1, 3, 5, 2, 5])

returns X = 5 and IX = 3.

   For a vector argument, return the minimum value.  For a matrix
argument, return the minimum value from each column, as a row vector.
Thus,

     min (min (X))

returns the smallest element of X.

   For complex arguments, the magnitude of the elements are used for
comparison.

   If called with two output arguments, also returns the index of the
minimum value(s). Thus,
     [x, ix] = min ([1, 3, 0, 2, 5])

returns X = 0 and IX = 3.

 - Mapping Function:  mod (X, Y)
     Compute modulo function, using

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

     Note that this handles negative numbers correctly: `mod (-1, 3)'
     is 2, not -1 as `rem (-1, 3)' returns.  Also, `mod (X, 0)' returns
     X.

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

 - Function File:  nextpow2 (X)
     If X is a scalar, returns the first integer N such that  2^n >=
     abs (x).

     If X is a vector, return `nextpow2 (length (X))'.

 - Mapping Function:  pow2 (X)
 - Mapping Function:  pow2 (F, E)
     With one argument, computes  2 .^ x for each element of X.  With
     two arguments, returns  f .* (2 .^ e).

 - Mapping Function:  rem (X, Y)
     Return the remainder of `X / Y', computed using the expression

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

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

 - Mapping Function:  round (X)
     Return the integer nearest to X.  If X is complex, return `round
     (real (X)) + round (imag (X)) * I'.

 - Mapping Function:  sign (X)
     Compute the "signum" function, which is defined as

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

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

 - Mapping Function:  sqrt (X)
     Compute the square root of X.  If X is negative, a complex result
     is returned.  To compute the matrix square root, see Linear Algebra.

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

   The following functions are available for working with complex
numbers.  Each expects a single argument.  Given a matrix they work on
an element by element basis.  In the descriptions of the following
functions, Z is the complex number X + IY, where I is defined as `sqrt
(-1)'.

 - Mapping Function:  abs (Z)
     Compute the magnitude of Z, defined as |Z| = `sqrt (x^2 + y^2)'.

     For example,

          abs (3 + 4i)
               => 5

 - Mapping Function:  arg (Z)
 - Mapping Function:  angle (Z)
     Compute the argument of Z, defined as THETA = `atan (Y/X)'.

     in radians.

     For example,

          arg (3 + 4i)
               => 0.92730

 - Mapping Function:  conj (Z)
     Return the complex conjugate of Z, defined as `conj (Z)' = X - IY.

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

 - Mapping Function:  real (Z)
     Return the real part of Z.

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

   Octave provides the following trigonometric functions.  Angles are
specified in radians.  To convert from degrees to radians multipy by
`pi/180'  (e.g. `sin (30 * pi/180)' returns the sine of 30 degrees).

 - Mapping Function:  sin (X)
     Compute the sin of each element of X.

 - Mapping Function:  cos (X)
     Compute the cosine of each element of X.

 - Mapping Function:  tan (Z)
     Compute tanget of each element of X.

 - Mapping Function:  sec (X)
     Compute the secant of each element of X.

 - Mapping Function:  csc (X)
     Compute the cosecant of each element of X.

 - Mapping Function:  cot (X)
     Compute the cotangent of each element of X.

 - Mapping Function:  asin (X)
     Compute the inverse sine of each element of X.

 - Mapping Function:  acos (X)
     Compute the inverse cosine of each element of X.

 - Mapping Function:  atan (X)
     Compute the inverse tangent of each element of X.

 - Mapping Function:  asec (X)
     Compute the inverse secant of each element of X.

 - Mapping Function:  acsc (X)
     Compute the inverse cosecant of each element of X.

 - Mapping Function:  acot (X)
     Compute the inverse cotangent of each element of X.

 - Mapping Function:  sinh (X)
     Compute the inverse hyperbolic sin of each element of X.

 - Mapping Function:  cosh (X)
     Compute the hyperbolic cosine of each element of X.

 - Mapping Function:  tanh (X)
     Compute hyperbolic tangent of each element of X.

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

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

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

 - Mapping Function:  asinh (X)
     Ompute the inverse hyperbolic sine of each element of X.

 - Mapping Function:  acosh (X)
     Compute the inverse hyperbolic cosine of each element of X.

 - Mapping Function:  atanh (X)
     Compute the inverse hyperbolic tanget of each element of X.

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

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

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

   Each of these functions expect a single argument.  For matrix
arguments, they work on an element by element basis.  For example,

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

 - Mapping Function:  atan2 (Y, X)
     Compute atan (Y / X) for corresponding elements of Y and X.  The
     result is in range -pi to pi.

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

 - Built-in Function:  sum (X, DIM)
     Sum of elements along dimension DIM.  If DIM is omitted, it
     defaults to 1 (column-wise sum).

 - Built-in Function:  prod (X, DIM)
     Product of elements along dimension DIM.  If DIM is omitted, it
     defaults to 1 (column-wise products).

 - Built-in Function:  cumsum (X, DIM)
     Cumulative sum of elements along dimension DIM.  If DIM is
     omitted, it defaults to 1 (column-wise cumulative sums).

 - Built-in Function:  cumprod (X, DIM)
     Cumulative product of elements along dimension DIM.  If DIM is
     omitted, it defaults to 1 (column-wise cumulative products).

 - Built-in Function:  sumsq (X, DIM)
     Sum of squares of elements along dimension DIM.  If DIM is
     omitted, it defaults to 1 (column-wise sum of squares).

     This function is equivalent to computing
          sum (x .* conj (x), dim)
     but it uses less memory and avoids calling conj if X is real.

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

 - Loadable Function: [J, IERR] = besselj (ALPHA, X, OPT)
 - Loadable Function: [Y, IERR] = bessely (ALPHA, X, OPT)
 - Loadable Function: [I, IERR] = besseli (ALPHA, X, OPT)
 - Loadable Function: [K, IERR] = besselk (ALPHA, X, OPT)
 - Loadable Function: [H, IERR] = besselh (ALPHA, K, X, OPT)
     Compute Bessel or Hankel functions of various kinds:

    `besselj'
          Bessel functions of the first kind.

    `bessely'
          Bessel functions of the second kind.

    `besseli'
          Modified Bessel functions of the first kind.

    `besselk'
          Modified Bessel functions of the second kind.

    `besselh'
          Compute Hankel functions of the first (K = 1) or second (K  =
          2) kind.

     If the argument OPT is supplied, the result is scaled by the `exp
     (-I*X)' for K = 1 or `exp (I*X)' for  K = 2.

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

     The value of ALPHA must be real.  The value of X may be complex.

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

       0. Normal return.

       1. Input error, return `NaN'.

       2. Overflow, return `Inf'.

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

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

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

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

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

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

     The result is the same size as Z.

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

       0. Normal return.

       1. Input error, return `NaN'.

       2. Overflow, return `Inf'.

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

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

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

 - Mapping Function:  beta (A, B)
     Return the Beta function,

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

 - Mapping Function:  betainc (X, A, B)
     Return the incomplete Beta function,

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

     If x has more than one component, both A and B must be scalars.
     If X is a scalar, A and B must be of compatible dimensions.

 - Mapping Function:  bincoeff (N, K)
     Return the binomial coefficient of N and K, defined as

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

     For example,

          bincoeff (5, 2)
          => 10

 - Mapping Function:  erf (Z)
     Computes the error function,

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

 - Mapping Function:  erfc (Z)
     Computes the complementary error function, `1 - erf (Z)'.

 - Mapping Function:  erfinv (Z)
     Computes the inverse of the error function.

 - Mapping Function:  gamma (Z)
     Computes the Gamma function,

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

 - Mapping Function:  gammainc (X, A)
     Computes the incomplete gamma function,

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

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

     If neither X nor A is scalar, the sizes of X and A must agree, and
     GAMMAINC is applied element-by-element.

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

 - Function File:  cross (X, Y)
     Computes the vector cross product of the two 3-dimensional vectors
     X and Y.

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

     If X and Y are two - dimensional matrices the cross product is
     applied along the first dimension with 3 elements.

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

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

     See Magnus and Neudecker (1988), Matrix differential calculus with
     applications in statistics and econometrics.

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

     See Magnus and Neudecker (1988), Matrix differential calculus with
     applications in statistics and econometrics.