~ubuntu-branches/debian/squeeze/maxima/squeeze

This is maxima.info, produced by makeinfo version 4.7 from maxima.texi.

   This is a Texinfo Maxima Manual

   Copyright 1994,2001 William F. Schelter

START-INFO-DIR-ENTRY
* Maxima: (maxima).     A computer algebra system.
END-INFO-DIR-ENTRY


File: maxima.info,  Node: Definitions for Differential Equations,  Prev: Differential Equations,  Up: Differential Equations

22.1 Definitions for Differential Equations
===========================================

 -- Function: bc2 (<solution>, <xval1>, <yval1>, <xval2>, <yval2>)
     Solves boundary value problem for second order differential
     equation.  Here: <solution> is a general solution to the equation,
     as found by `ode2', <xval1> is an equation for the independent
     variable in the form `<x> = <x0>', and <yval1> is an equation for
     the dependent variable in the form `<y> = <y0>'.  The <xval2> and
     <yval2> are equations for these variables at another point.  See
     `ode2' for example of usage.


 -- Function: desolve (<eqn>, <x>)
 -- Function: desolve ([<eqn_1>, ..., <eqn_n>], [<x_1>, ..., <x_n>])
     The function `dsolve' solves systems of linear ordinary
     differential equations using Laplace transform.  Here the <eqn>'s
     are differential equations in the dependent variables <x_1>, ...,
     <x_n>.  The functional relationships must be explicitly indicated
     in both the equations and the variables. For example

          'diff(f,x,2)=sin(x)+'diff(g,x);
          'diff(f,x)+x^2-f=2*'diff(g,x,2);

     is not the proper format.  The correct way is:

          'diff(f(x),x,2)=sin(x)+'diff(g(x),x);
          'diff(f(x),x)+x^2-f=2*'diff(g(x),x,2);

     The call is then `desolve([%o3,%o4],[f(x),g(x)]);' .

     If initial conditions at 0 are known, they should be supplied
     before calling `desolve' by using `atvalue'.

          (%i1) 'diff(f(x),x)='diff(g(x),x)+sin(x);
                           d           d
          (%o1)            -- (f(x)) = -- (g(x)) + sin(x)
                           dx          dx
          (%i2) 'diff(g(x),x,2)='diff(f(x),x)-cos(x);
                            2
                           d            d
          (%o2)            --- (g(x)) = -- (f(x)) - cos(x)
                             2          dx
                           dx
          (%i3) atvalue('diff(g(x),x),x=0,a);
          (%o3)                           a
          (%i4) atvalue(f(x),x=0,1);
          (%o4)                           1
          (%i5) desolve([%o1,%o2],[f(x),g(x)]);
                            x
          (%o5) [f(x) = a %e  - a + 1, g(x) =

                                                          x
                                             cos(x) + a %e  - a + g(0) - 1]
          (%i6) [%o1,%o2],%o5,diff;
                       x       x      x                x
          (%o6)   [a %e  = a %e , a %e  - cos(x) = a %e  - cos(x)]

     If `desolve' cannot obtain a solution, it returns `false'.


 -- Function: ic1 (<solution>, <xval>, <yval>)
     Solves initial value problem for first order differential equation.
     Here: <solution> is a general solution to the equation, as found
     by `ode2', <xval> is an equation for the independent variable in
     the form `<x> = <x0>', and <yval> is an equation for the dependent
     variable in the form `<y> = <y0>'. See `ode2' for example of usage.


 -- Function: ic2 (<solution>, <xval>, <yval>, <dval>)
     Solves initial value problem for second order differential
     equation.  Here: <solution> is a general solution to the equation,
     as found by `ode2', <xval> is an equation for the independent
     variable in the form `<x> = <x0>', <yval> is an equation for the
     dependent variable in the form `<y> = <y0>', and <dval> is an
     equation for the derivative of the dependent variable with respect
     to independent variable evaluated at the point <xval>.  See `ode2'
     for example of usage.


 -- Function: ode2 (<eqn>, <dvar>, <ivar>)
     The function `ode2' solves ordinary differential equations of
     first or second order.  It takes three arguments: an ODE <eqn>,
     the dependent variable <dvar>, and the independent variable <ivar>.
     When successful, it returns either an explicit or implicit
     solution for the dependent variable.  `%c' is used to represent
     the constant in the case of first order equations, and `%k1' and
     `%k2' the constants for second order equations.  If `ode2' cannot
     obtain a solution for whatever reason, it returns `false', after
     perhaps printing out an error message.  The methods implemented
     for first order equations in the order in which they are tested
     are: linear, separable, exact - perhaps requiring an integrating
     factor, homogeneous, Bernoulli's equation, and a generalized
     homogeneous method.  For second order: constant coefficient,
     exact, linear homogeneous with non-constant coefficients which can
     be transformed to constant coefficient, the Euler or
     equidimensional equation, the method of variation of parameters,
     and equations which are free of either the independent or of the
     dependent variable so that they can be reduced to two first order
     linear equations to be solved sequentially.  In the course of
     solving ODEs, several variables are set purely for informational
     purposes: `method' denotes the method of solution used e.g.
     `linear', `intfactor' denotes any integrating factor used,
     `odeindex' denotes the index for Bernoulli's method or for the
     generalized homogeneous method, and `yp' denotes the particular
     solution for the variation of parameters technique.

     In order to solve initial value problems (IVPs) and boundary value
     problems (BVPs), the routine `ic1' is available for first order
     equations, and `ic2' and `bc2' for second order IVPs and BVPs,
     respectively.

     Example:

          (%i1) x^2*'diff(y,x) + 3*y*x = sin(x)/x;
                                2 dy           sin(x)
          (%o1)                x  -- + 3 x y = ------
                                  dx             x
          (%i2) ode2(%,y,x);
                                       %c - cos(x)
          (%o2)                    y = -----------
                                            3
                                           x
          (%i3) ic1(%o2,x=%pi,y=0);
                                        cos(x) + 1
          (%o3)                   y = - ----------
                                             3
                                            x
          (%i4) 'diff(y,x,2) + y*'diff(y,x)^3 = 0;
                                   2
                                  d y      dy 3
          (%o4)                   --- + y (--)  = 0
                                    2      dx
                                  dx
          (%i5) ode2(%,y,x);
                                3
                               y  + 6 %k1 y
          (%o5)                ------------ = x + %k2
                                    6
          (%i6) ratsimp(ic2(%o5,x=0,y=0,'diff(y,x)=2));
                                       3
                                    2 y  - 3 y
          (%o6)                   - ---------- = x
                                        6
          (%i7) bc2(%o5,x=0,y=1,x=1,y=3);
                                   3
                                  y  - 10 y       3
          (%o7)                   --------- = x - -
                                      6           2



File: maxima.info,  Node: Numerical,  Next: Statistics,  Prev: Differential Equations,  Up: Top

23 Numerical
************

* Menu:

* Introduction to Numerical::
* Fourier packages::
* Definitions for Numerical::
* Definitions for Fourier Series::


File: maxima.info,  Node: Introduction to Numerical,  Next: Fourier packages,  Prev: Numerical,  Up: Numerical

23.1 Introduction to Numerical
==============================


File: maxima.info,  Node: Fourier packages,  Next: Definitions for Numerical,  Prev: Introduction to Numerical,  Up: Numerical

23.2 Fourier packages
=====================

The `fft' package comprises functions for the numerical (not symbolic)
computation of the fast Fourier transform.  `load ("fft")' loads this
package.  See `fft'.

   The `fourie' package comprises functions for the symbolic computation
of Fourier series.  `load ("fourie")' loads this package.  There are
functions in the `fourie' package to calculate Fourier integral
coefficients and some functions for manipulation of expressions.  See
`Definitions for Fourier Series'.


File: maxima.info,  Node: Definitions for Numerical,  Next: Definitions for Fourier Series,  Prev: Fourier packages,  Up: Numerical

23.3 Definitions for Numerical
==============================

 -- Function: polartorect (<magnitude_array>, <phase_array>)
     Translates complex values of the form `r %e^(%i t)' to the form `a
     + b %i'.  `load ("fft")' loads this function into Maxima. See also
     `fft'.

     The magnitude and phase, `r' and `t', are taken from
     <magnitude_array> and <phase_array>, respectively. The original
     values of the input arrays are replaced by the real and imaginary
     parts, `a' and `b', on return. The outputs are calculated as

          a: r cos (t)
          b: r sin (t)

     The input arrays must be the same size and 1-dimensional.  The
     array size need not be a power of 2.

     `polartorect' is the inverse function of `recttopolar'.


 -- Function: recttopolar (<real_array>, <imaginary_array>)
     Translates complex values of the form `a + b %i' to the form `r
     %e^(%i t)'.  `load ("fft")' loads this function into Maxima. See
     also `fft'.

     The real and imaginary parts, `a' and `b', are taken from
     <real_array> and <imaginary_array>, respectively. The original
     values of the input arrays are replaced by the magnitude and
     angle, `r' and `t', on return. The outputs are calculated as

          r: sqrt (a^2 + b^2)
          t: atan2 (b, a)

     The computed angle is in the range `-%pi' to `%pi'.

     The input arrays must be the same size and 1-dimensional.  The
     array size need not be a power of 2.

     `recttopolar' is the inverse function of `polartorect'.


 -- Function: ift (<real_array>, <imaginary_array>)
     Fast inverse discrete Fourier transform. `load ("fft")' loads this
     function into Maxima.

     `ift' carries out the inverse complex fast Fourier transform on
     1-dimensional floating point arrays. The inverse transform is
     defined as

          x[j]: sum (y[j] exp (+2 %i %pi j k / n), k, 0, n-1)

     See `fft' for more details.


 -- Function: fft (<real_array>, <imaginary_array>)
 -- Function: ift (<real_array>, <imaginary_array>)
 -- Function: recttopolar (<real_array>, <imaginary_array>)
 -- Function: polartorect (<magnitude_array>, <phase_array>)
     Fast Fourier transform and related functions. `load ("fft")' loads
     these functions into Maxima.

     `fft' and `ift' carry out the complex fast Fourier transform and
     inverse transform, respectively, on 1-dimensional floating point
     arrays. The size of <imaginary_array> must equal the size of
     <real_array>.

     `fft' and `ift' operate in-place. That is, on return from `fft' or
     `ift', the original content of the input arrays is replaced by the
     output.  The `fillarray' function can make a copy of an array,
     should it be necessary.

     The discrete Fourier transform and inverse transform are defined
     as follows. Let `x' be the original data, with

          x[i]: real_array[i] + %i imaginary_array[i]

     Let `y' be the transformed data. The forward and inverse
     transforms are

          y[k]: (1/n) sum (x[j] exp (-2 %i %pi j k / n), j, 0, n-1)

          x[j]:       sum (y[j] exp (+2 %i %pi j k / n), k, 0, n-1)

     Suitable arrays can be allocated by the `array' function. For
     example:

          array (my_array, float, n-1)$

     declares a 1-dimensional array with n elements, indexed from 0
     through n-1 inclusive. The number of elements n must be equal to
     2^m for some m.

     `fft' can be applied to real data (imaginary array all zeros) to
     obtain sine and cosine coefficients. After calling `fft', the sine
     and cosine coefficients, say `a' and `b', can be calculated as

          a[0]: real_array[0]
          b[0]: 0

     and

          a[j]: real_array[j] + real_array[n-j]
          b[j]: imaginary_array[j] - imaginary_array[n-j]

     for j equal to 1 through n/2-1, and

          a[n/2]: real_array[n/2]
          b[n/2]: 0

     `recttopolar' translates complex values of the form `a + b %i' to
     the form `r %e^(%i t)'. See `recttopolar'.

     `polartorect' translates complex values of the form `r %e^(%i t)'
     to the form `a + b %i'. See `polartorect'.

     `demo ("fft")' displays a demonstration of the `fft' package.


 -- Option variable: fortindent
     Default value: 0

     `fortindent' controls the left margin indentation of expressions
     printed out by the `fortran' command.  0 gives normal printout
     (i.e., 6 spaces), and positive values will causes the expressions
     to be printed farther to the right.


 -- Function: fortran (<expr>)
     Prints <expr> as a Fortran statement.  The output line is indented
     with spaces.  If the line is too long, `fortran' prints
     continuation lines.  `fortran' prints the exponentiation operator
     `^' as `**', and prints a complex number `a + b %i' in the form
     `(a,b)'.

     <expr> may be an equation. If so, `fortran' prints an assignment
     statement, assigning the right-hand side of the equation to the
     left-hand side.  In particular, if the right-hand side of <expr>
     is the name of a matrix, then `fortran' prints an assignment
     statement for each element of the matrix.

     If <expr> is not something recognized by `fortran', the expression
     is printed in `grind' format without complaint.  `fortran' does
     not know about lists, arrays, or functions.

     `fortindent' controls the left margin of the printed lines.  0 is
     the normal margin (i.e., indented 6 spaces). Increasing
     `fortindent' causes expressions to be printed further to the right.

     When `fortspaces' is `true', `fortran' fills out each printed line
     with spaces to 80 columns.

     `fortran' evaluates its arguments; quoting an argument defeats
     evaluation.  `fortran' always returns `done'.

     Examples:

          (%i1) expr: (a + b)^12$
          (%i2) fortran (expr);
                (b+a)**12
          (%o2)                         done
          (%i3) fortran ('x=expr);
                x = (b+a)**12
          (%o3)                         done
          (%i4) fortran ('x=expand (expr));
                x = b**12+12*a*b**11+66*a**2*b**10+220*a**3*b**9+495*a**4*b**8+792
               1   *a**5*b**7+924*a**6*b**6+792*a**7*b**5+495*a**8*b**4+220*a**9*b
               2   **3+66*a**10*b**2+12*a**11*b+a**12
          (%o4)                         done
          (%i5) fortran ('x=7+5*%i);
                x = (7,5)
          (%o5)                         done
          (%i6) fortran ('x=[1,2,3,4]);
                x = [1,2,3,4]
          (%o6)                         done
          (%i7) f(x) := x^2$
          (%i8) fortran (f);
                f
          (%o8)                         done


 -- Option variable: fortspaces
     Default value: `false'

     When `fortspaces' is `true', `fortran' fills out each printed line
     with spaces to 80 columns.


 -- Function: horner (<expr>, <x>)
 -- Function: horner (<expr>)
     Returns a rearranged representation of <expr> as in Horner's rule,
     using <x> as the main variable if it is specified.  `x' may be
     omitted in which case the main variable of the canonical rational
     expression form of <expr> is used.

     `horner' sometimes improves stability if `expr' is to be
     numerically evaluated.  It is also useful if Maxima is used to
     generate programs to be run in Fortran. See also `stringout'.

          (%i1) expr: 1e-155*x^2 - 5.5*x + 5.2e155;
                                     2
          (%o1)            1.0E-155 x  - 5.5 x + 5.2E+155
          (%i2) expr2: horner (%, x), keepfloat: true;
          (%o2)            (1.0E-155 x - 5.5) x + 5.2E+155
          (%i3) ev (expr, x=1e155);
          Maxima encountered a Lisp error:

           floating point overflow

          Automatically continuing.
          To reenable the Lisp debugger set *debugger-hook* to nil.
          (%i4) ev (expr2, x=1e155);
          (%o4)                       7.0E+154


 -- Function: find_root (<f>(<x>), <x>, <a>, <b>)
 -- Function: find_root (<f>, <a>, <b>)
     Finds the zero of function <f> as variable <x> varies over the
     range `[<a>, <b>]'.  The function must have a different sign at
     each endpoint.  If this condition is not met, the action of the
     function is governed by `find_root_error'.  If `find_root_error'
     is `true' then an error occurs, otherwise the value of
     `find_root_error' is returned (thus for plotting `find_root_error'
     might be set to 0.0).  Otherwise (given that Maxima can evaluate
     the first argument in the specified range, and that it is
     continuous) `find_root' is guaranteed to come up with the zero (or
     one of them if there is more than one zero).  The accuracy of
     `find_root' is governed by `find_root_abs' and `find_root_rel'
     which must be non-negative floating point numbers.  `find_root'
     will stop when the first arg evaluates to something less than or
     equal to `find_root_abs' or if successive approximants to the root
     differ by no more than `find_root_rel * <one of the approximants>'.
     The default values of `find_root_abs' and `find_root_rel' are 0.0
     so `find_root' gets as good an answer as is possible with the
     single precision arithmetic we have.  The first arg may be an
     equation.  The order of the last two args is irrelevant.  Thus

          find_root (sin(x) = x/2, x, %pi, 0.1);

     is equivalent to

          find_root (sin(x) = x/2, x, 0.1, %pi);

     The method used is a binary search in the range specified by the
     last two args.  When it thinks the function is close enough to
     being linear, it starts using linear interpolation.

     Examples:

          (%i1) f(x) := sin(x) - x/2;
                                                  x
          (%o1)                  f(x) := sin(x) - -
                                                  2
          (%i2) find_root (sin(x) - x/2, x, 0.1, %pi);
          (%o2)                   1.895494267033981
          (%i3) find_root (sin(x) = x/2, x, 0.1, %pi);
          (%o3)                   1.895494267033981
          (%i4) find_root (f(x), x, 0.1, %pi);
          (%o4)                   1.895494267033981
          (%i5) find_root (f, 0.1, %pi);
          (%o5)                   1.895494267033981


 -- Option variable: find_root_abs
     Default value: 0.0

     `find_root_abs' is the accuracy of the `find_root' command is
     governed by `find_root_abs' and `find_root_rel' which must be
     non-negative floating point numbers.  `find_root' will stop when
     the first arg evaluates to something less than or equal to
     `find_root_abs' or if successive approximants to the root differ
     by no more than `find_root_rel * <one of the approximants>'.  The
     default values of `find_root_abs' and `find_root_rel' are 0.0 so
     `find_root' gets as good an answer as is possible with the single
     precision arithmetic we have.


 -- Option variable: find_root_error
     Default value: `true'

     `find_root_error' governs the behavior of `find_root'.  When
     `find_root' is called, it determines whether or not the function
     to be solved satisfies the condition that the values of the
     function at the endpoints of the interpolation interval are
     opposite in sign.  If they are of opposite sign, the interpolation
     proceeds.  If they are of like sign, and `find_root_error' is
     `true', then an error is signaled.  If they are of like sign and
     `find_root_error' is not `true', the value of `find_root_error' is
     returned.  Thus for plotting, `find_root_error' might be set to
     0.0.


 -- Option variable: find_root_rel
     Default value: 0.0

     `find_root_rel' is the accuracy of the `find_root' command is
     governed by `find_root_abs' and `find_root_rel' which must be
     non-negative floating point numbers.  `find_root' will stop when
     the first arg evaluates to something less than or equal to
     `find_root_abs' or if successive approximants to the root differ
     by no more than `find_root_rel * <one of the approximants>'.  The
     default values of `find_root_abs' and `find_root_rel' are 0.0 so
     `find_root' gets as good an answer as is possible with the single
     precision arithmetic we have.



File: maxima.info,  Node: Definitions for Fourier Series,  Prev: Definitions for Numerical,  Up: Numerical

23.4 Definitions for Fourier Series
===================================

 -- Function: equalp (<x>, <y>)
     Returns `true' if `equal (<x>, <y>)' otherwise `false' (doesn't
     give an error message like `equal (x, y)' would do in this case).


 -- Function: remfun (<f>, <expr>)
 -- Function: remfun (<f>, <expr>, <x>)
     `remfun (<f>, <expr>)' replaces all occurrences of `<f> (<arg>)'
     by <arg> in <expr>.

     `remfun (<f>, <expr>, <x>)' replaces all occurrences of `<f>
     (<arg>)' by <arg> in <expr> only if <arg> contains the variable
     <x>.


 -- Function: funp (<f>, <expr>)
 -- Function: funp (<f>, <expr>, <x>)
     `funp (<f>, <expr>)' returns `true' if <expr> contains the
     function <f>.

     `funp (<f>, <expr>, <x>)' returns `true' if <expr> contains the
     function <f> and the variable <x> is somewhere in the argument of
     one of the instances of <f>.


 -- Function: absint (<f>, <x>, <halfplane>)
 -- Function: absint (<f>, <x>)
 -- Function: absint (<f>, <x>, <a>, <b>)
     `absint (<f>, <x>, <halfplane>)' returns the indefinite integral
     of <f> with respect to <x> in the given halfplane (`pos', `neg',
     or `both').  <f> may contain expressions of the form `abs (x)',
     `abs (sin (x))', `abs (a) * exp (-abs (b) * abs (x))'.

     `absint (<f>, <x>)' is equivalent to `absint (<f>, <x>, pos)'.

     `absint (<f>, <x>, <a>, <b>)' returns the definite integral of <f>
     with respect to <x> from <a> to <b>.  <f> may include absolute
     values.


 -- Function: fourier (<f>, <x>, <p>)
     Returns a list of the Fourier coefficients of `<f>(<x>)' defined
     on the interval `[-%pi, %pi]'.


 -- Function: foursimp (<l>)
     Simplifies `sin (n %pi)' to 0 if `sinnpiflag' is `true' and `cos
     (n %pi)' to `(-1)^n' if `cosnpiflag' is `true'.


 -- Option variable: sinnpiflag
     Default value: `true'

     See `foursimp'.


 -- Option variable: cosnpiflag
     Default value: `true'

     See `foursimp'.


 -- Function: fourexpand (<l>, <x>, <p>, <limit>)
     Constructs and returns the Fourier series from the list of Fourier
     coefficients <l> up through <limit> terms (<limit> may be `inf').
     <x> and <p> have same meaning as in `fourier'.


 -- Function: fourcos (<f>, <x>, <p>)
     Returns the Fourier cosine coefficients for `<f>(<x>)' defined on
     `[0, %pi]'.


 -- Function: foursin (<f>, <x>, <p>)
     Returns the Fourier sine coefficients for `<f>(<x>)' defined on
     `[0, %pi]'.


 -- Function: totalfourier (<f>, <x>, <p>)
     Returns `fourexpand (foursimp (fourier (<f>, <x>, <p>)), <x>, <p>,
     'inf)'.


 -- Function: fourint (<f>, <x>)
     Constructs and returns a list of the Fourier integral coefficients
     of `<f>(<x>)' defined on `[minf, inf]'.


 -- Function: fourintcos (<f>, <x>)
     Returns the Fourier cosine integral coefficients for `<f>(<x>)' on
     `[0, inf]'.


 -- Function: fourintsin (<f>, <x>)
     Returns the Fourier sine integral coefficients for `<f>(<x>)' on
     `[0, inf]'.



File: maxima.info,  Node: Statistics,  Next: Arrays,  Prev: Numerical,  Up: Top

24 Statistics
*************

* Menu:

* Definitions for Statistics::


File: maxima.info,  Node: Definitions for Statistics,  Prev: Statistics,  Up: Statistics

24.1 Definitions for Statistics
===============================

 -- Function: gauss (<mean>, <sd>)
     Returns a random floating point number from a normal distribution
     with mean <mean> and standard deviation <sd>.



File: maxima.info,  Node: Arrays,  Next: Matrices and Linear Algebra,  Prev: Statistics,  Up: Top

25 Arrays
*********

* Menu:

* Definitions for Arrays::


File: maxima.info,  Node: Definitions for Arrays,  Prev: Arrays,  Up: Arrays

25.1 Definitions for Arrays
===========================

 -- Function: array (<name>, <dim_1>, ..., <dim_n>)
 -- Function: array (<name>, <type>, <dim_1>, ..., <dim_n>)
 -- Function: array ([<name_1>, ..., <name_m>], <dim_1>, ..., <dim_n>)
     Creates an n-dimensional array.  n may be less than or equal to 5.
     The subscripts for the i'th dimension are the integers running
     from 0 to <dim_i>.

     `array (<name>, <dim_1>, ..., <dim_n>)' creates a general array.

     `array (<name>, <type>, <dim_1>, ..., <dim_n>)' creates an array,
     with elements of a specified type.  <type> can be `fixnum' for
     integers of limited size or `flonum' for floating-point numbers.

     `array ([<name_1>, ..., <name_m>], <dim_1>, ..., <dim_n>)' creates
     m arrays, all of the same dimensions.

     If the user assigns to a subscripted variable before declaring the
     corresponding array, an undeclared array is created.  Undeclared
     arrays, otherwise known as hashed arrays (because hash coding is
     done on the subscripts), are more general than declared arrays.
     The user does not declare their maximum size, and they grow
     dynamically by hashing as more elements are assigned values.  The
     subscripts of undeclared arrays need not even be numbers.  However,
     unless an array is rather sparse, it is probably more efficient to
     declare it when possible than to leave it undeclared.  The `array'
     function can be used to transform an undeclared array into a
     declared array.


 -- Function: arrayapply (<A>, [<i_1>, ..., <i_n>])
     Evaluates `<A> [<i_1>, ..., <i_n>]', where <A> is an array and
     <i_1>, ..., <i_n> are integers.

     This is reminiscent of `apply', except the first argument is an
     array instead of a function.


 -- Function: arrayinfo (<A>)
     Returns information about the array <A>.  The argument <A> may be
     a declared array, an undeclared (hashed) array, an array function,
     or a subscripted function.

     For declared arrays, `arrayinfo' returns a list comprising the
     atom `declared', the number of dimensions, and the size of each
     dimension.  The elements of the array, both bound and unbound, are
     returned by `listarray'.

     For undeclared arrays (hashed arrays), `arrayinfo' returns a list
     comprising the atom `hashed', the number of subscripts, and the
     subscripts of every element which has a value.  The values are
     returned by `listarray'.

     For array functions, `arrayinfo' returns a list comprising the
     atom `hashed', the number of subscripts, and any subscript values
     for which there are stored function values.  The stored function
     values are returned by `listarray'.

     For subscripted functions, `arrayinfo' returns a list comprising
     the atom `hashed', the number of subscripts, and any subscript
     values for which there are lambda expressions.  The lambda
     expressions are returned by `listarray'.

     Examples:

     `arrayinfo' and `listarray' applied to a declared array.

          (%i1) array (aa, 2, 3);
          (%o1)                          aa
          (%i2) aa [2, 3] : %pi;
          (%o2)                          %pi
          (%i3) aa [1, 2] : %e;
          (%o3)                          %e
          (%i4) arrayinfo (aa);
          (%o4)                 [declared, 2, [2, 3]]
          (%i5) listarray (aa);
          (%o5) [#####, #####, #####, #####, #####, #####, %e, #####,
                                                  #####, #####, #####, %pi]

     `arrayinfo' and `listarray' applied to an undeclared (hashed)
     array.

          (%i1) bb [FOO] : (a + b)^2;
                                             2
          (%o1)                       (b + a)
          (%i2) bb [BAR] : (c - d)^3;
                                             3
          (%o2)                       (c - d)
          (%i3) arrayinfo (bb);
          (%o3)               [hashed, 1, [BAR], [FOO]]
          (%i4) listarray (bb);
                                        3         2
          (%o4)                 [(c - d) , (b + a) ]

     `arrayinfo' and `listarray' applied to an array function.

          (%i1) cc [x, y] := y / x;
                                               y
          (%o1)                      cc     := -
                                       x, y    x
          (%i2) cc [u, v];
                                          v
          (%o2)                           -
                                          u
          (%i3) cc [4, z];
                                          z
          (%o3)                           -
                                          4
          (%i4) arrayinfo (cc);
          (%o4)              [hashed, 2, [4, z], [u, v]]
          (%i5) listarray (cc);
                                        z  v
          (%o5)                        [-, -]
                                        4  u

     `arrayinfo' and `listarray' applied to a subscripted function.

          (%i1) dd [x] (y) := y ^ x;
                                               x
          (%o1)                     dd (y) := y
                                      x
          (%i2) dd [a + b];
                                              b + a
          (%o2)                  lambda([y], y     )
          (%i3) dd [v - u];
                                              v - u
          (%o3)                  lambda([y], y     )
          (%i4) arrayinfo (dd);
          (%o4)             [hashed, 1, [b + a], [v - u]]
          (%i5) listarray (dd);
                                   b + a                v - u
          (%o5)      [lambda([y], y     ), lambda([y], y     )]


 -- Function: arraymake (<name>, [<i_1>, ..., <i_n>])
     Returns the expression `<name> [<i_1>, ..., <i_n>]'.

     This is reminiscent of `funmake', except the return value is an
     unevaluated array reference instead of an unevaluated function
     call.


 -- System variable: arrays
     Default value: `[]'

     `arrays' is a list of arrays that have been allocated.  These
     comprise arrays declared by `array', hashed arrays constructed by
     implicit definition (assigning something to an array element), and
     array functions defined by `:=' and `define'.  Arrays defined by
     `make_array' are not included.

     See also `array', `arrayapply', `arrayinfo', `arraymake',
     `fillarray', `listarray', and `rearray'.

     Examples:

          (%i1) array (aa, 5, 7);
          (%o1)                          aa
          (%i2) bb [FOO] : (a + b)^2;
                                             2
          (%o2)                       (b + a)
          (%i3) cc [x] := x/100;
                                             x
          (%o3)                      cc  := ---
                                       x    100
          (%i4) dd : make_array ('any, 7);
          (%o4)       {Array:  #(NIL NIL NIL NIL NIL NIL NIL)}
          (%i5) arrays;
          (%o5)                     [aa, bb, cc]


 -- Function: bashindices (<expr>)
     Transforms the expression <expr> by giving each summation and
     product a unique index. This gives `changevar' greater precision
     when it is working with summations or products.  The form of the
     unique index is `j<number>'. The quantity <number> is determined by
     referring to `gensumnum', which can be changed by the user.  For
     example, `gensumnum:0$' resets it.


 -- Function: fillarray (<A>, <B>)
     Fills array <A> from <B>, which is a list or an array.

     If <A> is a floating-point (integer) array then <B> should be
     either a list of floating-point (integer) numbers or another
     floating-point (integer) array.

     If the dimensions of the arrays are different <A> is filled in
     row-major order.  If there are not enough elements in <B> the last
     element is used to fill out the rest of <A>.  If there are too
     many the remaining ones are thrown away.

     `fillarray' returns its first argument.


 -- Function: listarray (<A>)
     Returns a list of the elements of the array <A>.  The argument <A>
     may be a declared array, an undeclared (hashed) array, an array
     function, or a subscripted function.

     Elements are listed in row-major order.  That is, elements are
     sorted according to the first index, then according to the second
     index, and so on.  The sorting order of index values is the same
     as the order established by `orderless'.

     For undeclared arrays, array functions, and subscripted functions,
     the elements correspond to the index values returned by
     `arrayinfo'.

     Unbound elements of declared general arrays (that is, not `fixnum'
     and not `flonum') are returned as `#####'.  Unbound elements of
     declared `fixnum' or `flonum' arrays are returned as 0 or 0.0,
     respectively.  Unbound elements of undeclared arrays, array
     functions, and subscripted functions are not returned.

     Examples:

     `listarray' and `arrayinfo' applied to a declared array.

          (%i1) array (aa, 2, 3);
          (%o1)                          aa
          (%i2) aa [2, 3] : %pi;
          (%o2)                          %pi
          (%i3) aa [1, 2] : %e;
          (%o3)                          %e
          (%i4) listarray (aa);
          (%o4) [#####, #####, #####, #####, #####, #####, %e, #####,
                                                  #####, #####, #####, %pi]
          (%i5) arrayinfo (aa);
          (%o5)                 [declared, 2, [2, 3]]

     `listarray' and `arrayinfo' applied to an undeclared (hashed)
     array.

          (%i1) bb [FOO] : (a + b)^2;
                                             2
          (%o1)                       (b + a)
          (%i2) bb [BAR] : (c - d)^3;
                                             3
          (%o2)                       (c - d)
          (%i3) listarray (bb);
                                        3         2
          (%o3)                 [(c - d) , (b + a) ]
          (%i4) arrayinfo (bb);
          (%o4)               [hashed, 1, [BAR], [FOO]]

     `listarray' and `arrayinfo' applied to an array function.

          (%i1) cc [x, y] := y / x;
                                               y
          (%o1)                      cc     := -
                                       x, y    x
          (%i2) cc [u, v];
                                          v
          (%o2)                           -
                                          u
          (%i3) cc [4, z];
                                          z
          (%o3)                           -
                                          4
          (%i4) listarray (cc);
                                        z  v
          (%o4)                        [-, -]
                                        4  u
          (%i5) arrayinfo (cc);
          (%o5)              [hashed, 2, [4, z], [u, v]]

     `listarray' and `arrayinfo' applied to a subscripted function.

          (%i1) dd [x] (y) := y ^ x;
                                               x
          (%o1)                     dd (y) := y
                                      x
          (%i2) dd [a + b];
                                              b + a
          (%o2)                  lambda([y], y     )
          (%i3) dd [v - u];
                                              v - u
          (%o3)                  lambda([y], y     )
          (%i4) listarray (dd);
                                   b + a                v - u
          (%o4)      [lambda([y], y     ), lambda([y], y     )]
          (%i5) arrayinfo (dd);
          (%o5)             [hashed, 1, [b + a], [v - u]]


 -- Function: make_array (<type>, <dim_1>, ..., <dim_n>)
     Creates and returns a Lisp array.  <type> may be `any', `flonum',
     `fixnum', `hashed' or `functional'.  There are n indices, and the
     i'th index runs from 0 to <dim_i> - 1.

     The advantage of `make_array' over `array' is that the return
     value doesn't have a name, and once a pointer to it goes away, it
     will also go away.  For example, if `y: make_array (...)' then `y'
     points to an object which takes up space, but after `y: false',
     `y' no longer points to that object, so the object can be garbage
     collected.


 -- Function: rearray (<A>, <dim_1>, ..., <dim_n>)
     Changes the dimensions of an array.  The new array will be filled
     with the elements of the old one in row-major order.  If the old
     array was too small, the remaining elements are filled with
     `false', `0.0' or `0', depending on the type of the array.  The
     type of the array cannot be changed.


 -- Function: remarray (<A_1>, ..., <A_n>)
 -- Function: remarray (all)
     Removes arrays and array associated functions and frees the
     storage occupied.  The arguments may be declared arrays,
     undeclared (hashed) arrays, array functions, and subscripted
     functions.

     `remarray (all)' removes all items in the global list `arrays'.

     It may be necessary to use this function if it is desired to
     redefine the values in a hashed array.

     `remarray' returns the list of arrays removed.


 -- Function: subvar (<x>, <i>)
     Evaluates the subscripted expression `<x>[<i>]'.

     `subvar' evaluates its arguments.

     `arraymake (<x>, [<i>]' constructs the expression `<x>[<i>]', but
     does not evaluate it.

     Examples:

          (%i1) x : foo $

          (%i2) i : 3 $

          (%i3) subvar (x, i);
          (%o3)                         foo
                                           3
          (%i4) foo : [aa, bb, cc, dd, ee]$

          (%i5) subvar (x, i);
          (%o5)                          cc
          (%i6) arraymake (x, [i]);
          (%o6)                         foo
                                           3
          (%i7) ''%;
          (%o7)                          cc


 -- Option variable: use_fast_arrays
     - if `true' then only two types of arrays are recognized.

     1) The art-q array (t in Common Lisp) which may have several
     dimensions indexed by integers, and may hold any Lisp or Maxima
     object as an entry.  To construct such an array, enter
     `a:make_array(any,3,4);' then `a' will have as value, an array
     with twelve slots, and the indexing is zero based.

     2) The Hash_table array which is the default type of array created
     if one does `b[x+1]:y^2' (and `b' is not already an array, a list,
     or a matrix - if it were one of these an error would be caused
     since `x+1' would not be a valid subscript for an art-q array, a
     list or a matrix). Its indices (also known as keys) may be any
     object.  It only takes one key at a time (`b[x+1,u]:y' would
     ignore the `u').  Referencing is done by `b[x+1] ==> y^2'.  Of
     course the key may be a list, e.g. `b[[x+1,u]]:y' would be valid.
     This is incompatible with the old Maxima hash arrays, but saves
     consing.

     An advantage of storing the arrays as values of the symbol is that
     the usual conventions about local variables of a function apply to
     arrays as well.  The Hash_table type also uses less consing and is
     more efficient than the old type of Maxima hashar.  To obtain
     consistent behaviour in translated and compiled code set
     `translate_fast_arrays' to be `true'.



File: maxima.info,  Node: Matrices and Linear Algebra,  Next: Affine,  Prev: Arrays,  Up: Top

26 Matrices and Linear Algebra
******************************

* Menu:

* Introduction to Matrices and Linear Algebra::
* Definitions for Matrices and Linear Algebra::


File: maxima.info,  Node: Introduction to Matrices and Linear Algebra,  Next: Definitions for Matrices and Linear Algebra,  Prev: Matrices and Linear Algebra,  Up: Matrices and Linear Algebra

26.1 Introduction to Matrices and Linear Algebra
================================================

* Menu:

* Dot::
* Vectors::
* eigen::


File: maxima.info,  Node: Dot,  Next: Vectors,  Prev: Introduction to Matrices and Linear Algebra,  Up: Introduction to Matrices and Linear Algebra

26.1.1 Dot
----------

The operator `.' represents noncommutative multiplication and scalar
product.  When the operands are 1-column or 1-row matrices `a' and `b',
the expression `a.b' is equivalent to `sum (a[i]*b[i], i, 1,
length(a))'.  If `a' and `b' are not complex, this is the scalar
product, also called the inner product or dot product, of `a' and `b'.
The scalar product is defined as `conjugate(a).b' when `a' and `b' are
complex; `innerproduct' in the `eigen' package provides the complex
scalar product.

   When the operands are more general matrices, the product is the
matrix product `a' and `b'.  The number of rows of `b' must equal the
number of columns of `a', and the result has number of rows equal to
the number of rows of `a' and number of columns equal to the number of
columns of `b'.

   To distinguish `.' as an arithmetic operator from the decimal point
in a floating point number, it may be necessary to leave spaces on
either side.  For example, `5.e3' is `5000.0' but `5 . e3' is `5' times
`e3'.

   There are several flags which govern the simplification of
expressions involving `.', namely `dot', `dot0nscsimp', `dot0simp',
`dot1simp', `dotassoc', `dotconstrules', `dotdistrib', `dotexptsimp',
`dotident', and `dotscrules'.


File: maxima.info,  Node: Vectors,  Next: eigen,  Prev: Dot,  Up: Introduction to Matrices and Linear Algebra

26.1.2 Vectors
--------------

`vect' is a package of functions for vector analysis.  `load ("vect")'
loads this package, and `demo ("vect")' displays a demonstration.

   The vector analysis package can combine and simplify symbolic
expressions including dot products and cross products, together with
the gradient, divergence, curl, and Laplacian operators.  The
distribution of these operators over sums or products is governed by
several flags, as are various other expansions, including expansion
into components in any specific orthogonal coordinate systems.  There
are also functions for deriving the scalar or vector potential of a
field.

   The `vect' package contains these functions: `vectorsimp',
`scalefactors', `express', `potential', and `vectorpotential'.

   Warning: the `vect' package declares the dot operator `.' to be a
commutative operator.


File: maxima.info,  Node: eigen,  Prev: Vectors,  Up: Introduction to Matrices and Linear Algebra

26.1.3 eigen
------------

The package `eigen' contains several functions devoted to the symbolic
computation of eigenvalues and eigenvectors.  Maxima loads the package
automatically if one of the functions `eigenvalues' or `eigenvectors'
is invoked.  The package may be loaded explicitly as `load ("eigen")'.

   `demo ("eigen")' displays a demonstration of the capabilities of
this package.  `batch ("eigen")' executes the same demonstration, but
without the user prompt between successive computations.

   The functions in the `eigen' package are `innerproduct',
`unitvector', `columnvector', `gramschmidt', `eigenvalues',
`eigenvectors', `uniteigenvectors', and `similaritytransform'.


File: maxima.info,  Node: Definitions for Matrices and Linear Algebra,  Prev: Introduction to Matrices and Linear Algebra,  Up: Matrices and Linear Algebra

26.2 Definitions for Matrices and Linear Algebra
================================================

 -- Function: addcol (<M>, <list_1>, ..., <list_n>)
     Appends the column(s) given by the one or more lists (or matrices)
     onto the matrix <M>.


 -- Function: addrow (<M>, <list_1>, ..., <list_n>)
     Appends the row(s) given by the one or more lists (or matrices)
     onto the matrix <M>.


 -- Function: adjoint (<M>)
     Returns the adjoint of the matrix <M>.  The adjoint matrix is the
     transpose of the matrix of cofactors of <M>.


 -- Function: augcoefmatrix ([<eqn_1>, ..., <eqn_m>], [<x_1>, ...,
          <x_n>])
     Returns the augmented coefficient matrix for the variables <x_1>,
     ..., <x_n> of the system of linear equations <eqn_1>, ...,
     <eqn_m>.  This is the coefficient matrix with a column adjoined for
     the constant terms in each equation (i.e., those terms not
     dependent upon <x_1>, ..., <x_n>).

          (%i1) m: [2*x - (a - 1)*y = 5*b, c + b*y + a*x = 0]$
          (%i2) augcoefmatrix (m, [x, y]);
                                 [ 2  1 - a  - 5 b ]
          (%o2)                  [                 ]
                                 [ a    b      c   ]


 -- Function: charpoly (<M>, <x>)
     Returns the characteristic polynomial for the matrix <M> with
     respect to variable <x>.  That is, `determinant (<M> - diagmatrix
     (length (<M>), <x>))'.

          (%i1) a: matrix ([3, 1], [2, 4]);
                                      [ 3  1 ]
          (%o1)                       [      ]
                                      [ 2  4 ]
          (%i2) expand (charpoly (a, lambda));
                                     2
          (%o2)                lambda  - 7 lambda + 10
          (%i3) (programmode: true, solve (%));
          (%o3)               [lambda = 5, lambda = 2]
          (%i4) matrix ([x1], [x2]);
                                       [ x1 ]
          (%o4)                        [    ]
                                       [ x2 ]
          (%i5) ev (a . % - lambda*%, %th(2)[1]);
                                    [ x2 - 2 x1 ]
          (%o5)                     [           ]
                                    [ 2 x1 - x2 ]
          (%i6) %[1, 1] = 0;
          (%o6)                     x2 - 2 x1 = 0
          (%i7) x2^2 + x1^2 = 1;
                                      2     2
          (%o7)                     x2  + x1  = 1
          (%i8) solve ([%th(2), %], [x1, x2]);
                            1               2
          (%o8) [[x1 = - -------, x2 = - -------],
                         sqrt(5)         sqrt(5)

                                                       1             2
                                              [x1 = -------, x2 = -------]]
                                                    sqrt(5)       sqrt(5)


 -- Function: coefmatrix ([<eqn_1>, ..., <eqn_m>], [<x_1>, ..., <x_n>])
     Returns the coefficient matrix for the variables <x_1>, ..., <x_n>
     of the system of linear equations <eqn_1>, ..., <eqn_m>.

          (%i1) coefmatrix([2*x-(a-1)*y+5*b = 0, b*y+a*x = 3], [x,y]);
                                           [ 2  1 - a ]
          (%o1)                            [          ]
                                           [ a    b   ]


 -- Function: col (<M>, <i>)
     Returns the <i>'th column of the matrix <M>.  The return value is
     a matrix.


 -- Function: columnvector (<L>)
 -- Function: covect (<L>)
     Returns a matrix of one column and `length (<L>)' rows, containing
     the elements of the list <L>.

     `covect' is a synonym for `columnvector'.

     `load ("eigen")' loads this function.

     This is useful if you want to use parts of the outputs of the
     functions in this package in matrix calculations.

     Example:

          (%i1) load ("eigen")$
          Warning - you are redefining the Macsyma function eigenvalues
          Warning - you are redefining the Macsyma function eigenvectors
          (%i2) columnvector ([aa, bb, cc, dd]);
                                       [ aa ]
                                       [    ]
                                       [ bb ]
          (%o2)                        [    ]
                                       [ cc ]
                                       [    ]
                                       [ dd ]


 -- Function: conjugate (<x>)
     Returns the complex conjugate of <x>.

          (%i1) declare ([aa, bb], real, cc, complex, ii, imaginary);

          (%o1)                         done
          (%i2) conjugate (aa + bb*%i);

          (%o2)                      aa - %i bb
          (%i3) conjugate (cc);

          (%o3)                     conjugate(cc)
          (%i4) conjugate (ii);

          (%o4)                         - ii
          (%i5) conjugate (xx + yy);

          (%o5)             conjugate(yy) + conjugate(xx)


 -- Function: copymatrix (<M>)
     Returns a copy of the matrix <M>.  This is the only way to make a
     copy aside from copying <M> element by element.

     Note that an assignment of one matrix to another, as in `m2: m1',
     does not copy `m1'.  An assignment `m2 [i,j]: x' or `setelmx (x,
     i, j, m2' also modifies `m1 [i,j]'.  Creating a copy with
     `copymatrix' and then using assignment creates a separate,
     modified copy.


 -- Function: determinant (<M>)
     Computes the determinant of <M> by a method similar to Gaussian
     elimination.

     The form of the result depends upon the setting of the switch
     `ratmx'.

     There is a special routine for computing sparse determinants which
     is called when the switches `ratmx' and `sparse' are both `true'.


 -- Option variable: detout
     Default value: `false'

     When `detout' is `true', the determinant of a matrix whose inverse
     is computed is factored out of the inverse.

     For this switch to have an effect `doallmxops' and `doscmxops'
     should be `false' (see their descriptions).  Alternatively this
     switch can be given to `ev' which causes the other two to be set
     correctly.

     Example:

          (%i1) m: matrix ([a, b], [c, d]);
                                      [ a  b ]
          (%o1)                       [      ]
                                      [ c  d ]
          (%i2) detout: true$
          (%i3) doallmxops: false$
          (%i4) doscmxops: false$
          (%i5) invert (m);
                                    [  d   - b ]
                                    [          ]
                                    [ - c   a  ]
          (%o5)                     ------------
                                     a d - b c


 -- Function: diagmatrix (<n>, <x>)
     Returns a diagonal matrix of size <n> by <n> with the diagonal
     elements all equal to <x>.  `diagmatrix (<n>, 1)' returns an
     identity matrix (same as `ident (<n>)').

     <n> must evaluate to an integer, otherwise `diagmatrix' complains
     with an error message.

     <x> can be any kind of expression, including another matrix.  If
     <x> is a matrix, it is not copied; all diagonal elements refer to
     the same instance, <x>.


 -- Option variable: doallmxops
     Default value: `true'

     When `doallmxops' is `true', all operations relating to matrices
     are carried out.  When it is `false' then the setting of the
     individual `dot' switches govern which operations are performed.


 -- Option variable: domxexpt
     Default value: `true'

     When `domxexpt' is `true', a matrix exponential, `exp (<M>)' where
     <M> is a matrix, is interpreted as a matrix with element `[i,j'
     equal to `exp (m[i,j])'.  Otherwise `exp (<M>)' evaluates to `exp
     (<ev(M)>'.

     `domxexpt' affects all expressions of the form `<base>^<power>'
     where <base> is an expression assumed scalar or constant, and
     <power> is a list or matrix.

     Example:

          (%i1) m: matrix ([1, %i], [a+b, %pi]);
                                   [   1    %i  ]
          (%o1)                    [            ]
                                   [ b + a  %pi ]
          (%i2) domxexpt: false$
          (%i3) (1 - c)^m;
                                       [   1    %i  ]
                                       [            ]
                                       [ b + a  %pi ]
          (%o3)                 (1 - c)
          (%i4) domxexpt: true$
          (%i5) (1 - c)^m;
                            [                      %i  ]
                            [    1 - c      (1 - c)    ]
          (%o5)             [                          ]
                            [        b + a         %pi ]
                            [ (1 - c)       (1 - c)    ]


 -- Option variable: domxmxops
     Default value: `true'

     When `domxmxops' is `true', all matrix-matrix or matrix-list
     operations are carried out (but not scalar-matrix operations); if
     this switch is `false' such operations are not carried out.


 -- Option variable: domxnctimes
     Default value: `false'

     When `domxnctimes' is `true', non-commutative products of matrices
     are carried out.


 -- Option variable: dontfactor
     Default value: `[]'

     `dontfactor' may be set to a list of variables with respect to
     which factoring is not to occur.  (The list is initially empty.)
     Factoring also will not take place with respect to any variables
     which are less important, according the variable ordering assumed
     for canonical rational expression (CRE) form, than those on the
     `dontfactor' list.


 -- Option variable: doscmxops
     Default value: `false'

     When `doscmxops' is `true', scalar-matrix operations are carried
     out.


 -- Option variable: doscmxplus
     Default value: `false'

     When `doscmxplus' is `true', scalar-matrix operations yield a
     matrix result.  This switch is not subsumed under `doallmxops'.


 -- Option variable: dot0nscsimp
     Default value: `true'

     When `dot0nscsimp' is `true', a non-commutative product of zero
     and a nonscalar term is simplified to a commutative product.


 -- Option variable: dot0simp
     Default value: `true'

     When `dot0simp' is `true', a non-commutative product of zero and a
     scalar term is simplified to a commutative product.


 -- Option variable: dot1simp
     Default value: `true'

     When `dot1simp' is `true', a non-commutative product of one and
     another term is simplified to a commutative product.


 -- Option variable: dotassoc
     Default value: `true'

     When `dotassoc' is `true', an expression `(A.B).C' simplifies to
     `A.(B.C)'.


 -- Option variable: dotconstrules
     Default value: `true'

     When `dotconstrules' is `true', a non-commutative product of a
     constant and another term is simplified to a commutative product.
     Turning on this flag effectively turns on `dot0simp',
     `dot0nscsimp', and `dot1simp' as well.


 -- Option variable: dotdistrib
     Default value: `false'

     When `dotdistrib' is `true', an expression `A.(B + C)' simplifies
     to `A.B + A.C'.


 -- Option variable: dotexptsimp
     Default value: `true'

     When `dotexptsimp' is `true', an expression `A.A' simplifies to
     `A^^2'.


 -- Option variable: dotident
     Default value: 1

     `dotident' is the value returned by `X^^0'.


 -- Option variable: dotscrules
     Default value: `false'

     When `dotscrules' is `true', an expression `A.SC' or `SC.A'
     simplifies to `SC*A' and `A.(SC*B)' simplifies to `SC*(A.B)'.


 -- Function: echelon (<M>)
     Returns the echelon form of the matrix <M>, as produced by
     Gaussian elimination.  The echelon form is computed from <M> by
     elementary row operations such that the first non-zero element in
     each row in the resulting matrix is one and the column elements
     under the first one in each row are all zero.

     `triangularize' also carries out Gaussian elimination, but it does
     not normalize the leading non-zero element in each row.

     `lu_factor' and `cholesky' are other functions which yield
     triangularized matrices.

          (%i1) M: matrix ([3, 7, aa, bb], [-1, 8, 5, 2], [9, 2, 11, 4]);
                                 [  3   7  aa  bb ]
                                 [                ]
          (%o1)                  [ - 1  8  5   2  ]
                                 [                ]
                                 [  9   2  11  4  ]
          (%i2) echelon (M);
                            [ 1  - 8  - 5      - 2     ]
                            [                          ]
                            [         28       11      ]
                            [ 0   1   --       --      ]
          (%o2)             [         37       37      ]
                            [                          ]
                            [              37 bb - 119 ]
                            [ 0   0    1   ----------- ]
                            [              37 aa - 313 ]


 -- Function: eigenvalues (<M>)
 -- Function: eivals (<M>)
     Returns a list of two lists containing the eigenvalues of the
     matrix <M>.  The first sublist of the return value is the list of
     eigenvalues of the matrix, and the second sublist is the list of
     the multiplicities of the eigenvalues in the corresponding order.

     `eivals' is a synonym for `eigenvalues'.

     `eigenvalues' calls the function `solve' to find the roots of the
     characteristic polynomial of the matrix.  Sometimes `solve' may
     not be able to find the roots of the polynomial; in that case some
     other functions in this package (except `innerproduct',
     `unitvector', `columnvector' and `gramschmidt') will not work.

     In some cases the eigenvalues found by `solve' may be complicated
     expressions.  (This may happen when `solve' returns a
     not-so-obviously real expression for an eigenvalue which is known
     to be real.)  It may be possible to simplify the eigenvalues using
     some other functions.

     The package `eigen.mac' is loaded automatically when `eigenvalues'
     or `eigenvectors' is referenced.  If `eigen.mac' is not already
     loaded, `load ("eigen")' loads it.  After loading, all functions
     and variables in the package are available.


 -- Function: eigenvectors (<M>)
 -- Function: eivects (<M>)
     takes a matrix <M> as its argument and returns a list of lists the
     first sublist of which is the output of `eigenvalues' and the
     other sublists of which are the eigenvectors of the matrix
     corresponding to those eigenvalues respectively.  The calculated
     eigenvectors and the unit eigenvectors of the matrix are the right
     eigenvectors and the right unit eigenvectors respectively.

     `eivects' is a synonym for `eigenvectors'.

     The package `eigen.mac' is loaded automatically when `eigenvalues'
     or `eigenvectors' is referenced.  If `eigen.mac' is not already
     loaded, `load ("eigen")' loads it.  After loading, all functions
     and variables in the package are available.

     The flags that affect this function are:

     `nondiagonalizable' is set to `true' or `false' depending on
     whether the matrix is nondiagonalizable or diagonalizable after
     `eigenvectors' returns.

     `hermitianmatrix' when `true', causes the degenerate eigenvectors
     of the Hermitian matrix to be orthogonalized using the
     Gram-Schmidt algorithm.

     `knowneigvals' when `true' causes the `eigen' package to assume the
     eigenvalues of the matrix are known to the user and stored under
     the global name `listeigvals'.  `listeigvals' should be set to a
     list similar to the output `eigenvalues'.

     The function `algsys' is used here to solve for the eigenvectors.
     Sometimes if the eigenvalues are messy, `algsys' may not be able
     to find a solution.  In some cases, it may be possible to simplify
     the eigenvalues by first finding them using `eigenvalues' command
     and then using other functions to reduce them to something simpler.
     Following simplification, `eigenvectors' can be called again with
     the `knowneigvals' flag set to `true'.


 -- Function: ematrix (<m>, <n>, <x>, <i>, <j>)
     Returns an <m> by <n> matrix, all elements of which are zero
     except for the `[<i>, <j>]' element which is <x>.


 -- Function: entermatrix (<m>, <n>)
     Returns an <m> by <n> matrix, reading the elements interactively.

     If <n> is equal to <m>, Maxima prompts for the type of the matrix
     (diagonal, symmetric, antisymmetric, or general) and for each
     element.  Each response is terminated by a semicolon `;' or dollar
     sign `$'.

     If <n> is not equal to <m>, Maxima prompts for each element.

     The elements may be any expressions, which are evaluated.
     `entermatrix' evaluates its arguments.

          (%i1) n: 3$
          (%i2) m: entermatrix (n, n)$

          Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
          Answer 1, 2, 3 or 4 :
          1$
          Row 1 Column 1:
          (a+b)^n$
          Row 2 Column 2:
          (a+b)^(n+1)$
          Row 3 Column 3:
          (a+b)^(n+2)$

          Matrix entered.
          (%i3) m;
                          [        3                     ]
                          [ (b + a)      0         0     ]
                          [                              ]
          (%o3)           [                  4           ]
                          [    0      (b + a)      0     ]
                          [                              ]
                          [                            5 ]
                          [    0         0      (b + a)  ]


 -- Function: genmatrix (<a>, <i_2>, <j_2>, <i_1>, <j_1>)
 -- Function: genmatrix (<a>, <i_2>, <j_2>, <i_1>)
 -- Function: genmatrix (<a>, <i_2>, <j_2>)
     Returns a matrix generated from <a>, taking element
     `<a>[<i_1>,<j_1>]' as the upper-left element and `<a>[<i_2>,<j_2>]'
     as the lower-right element of the matrix.  Here <a> is a declared
     array (created by `array' but not by `make_array') or an
     undeclared array, or an array function, or a lambda expression of
     two arguments.  (An array function is created like other functions
     with `:=' or `define', but arguments are enclosed in square
     brackets instead of parentheses.)

     If <j_1> is omitted, it is assumed equal to <i_1>.  If both <j_1>
     and <i_1> are omitted, both are assumed equal to 1.

     If a selected element `i,j' of the array is undefined, the matrix
     will contain a symbolic element `<a>[i,j]'.

     Examples:

          (%i1) h [i, j] := 1 / (i + j - 1);
                                              1
          (%o1)                  h     := ---------
                                  i, j    i + j - 1
          (%i2) genmatrix (h, 3, 3);
                                     [    1  1 ]
                                     [ 1  -  - ]
                                     [    2  3 ]
                                     [         ]
                                     [ 1  1  1 ]
          (%o2)                      [ -  -  - ]
                                     [ 2  3  4 ]
                                     [         ]
                                     [ 1  1  1 ]
                                     [ -  -  - ]
                                     [ 3  4  5 ]
          (%i3) array (a, fixnum, 2, 2);
          (%o3)                           a
          (%i4) a [1, 1] : %e;
          (%o4)                          %e
          (%i5) a [2, 2] : %pi;
          (%o5)                          %pi
          (%i6) genmatrix (a, 2, 2);
                                     [ %e   0  ]
          (%o6)                      [         ]
                                     [ 0   %pi ]
          (%i7) genmatrix (lambda ([i, j], j - i), 3, 3);
                                   [  0    1   2 ]
                                   [             ]
          (%o7)                    [ - 1   0   1 ]
                                   [             ]
                                   [ - 2  - 1  0 ]
          (%i8) genmatrix (B, 2, 2);
                                  [ B      B     ]
                                  [  1, 1   1, 2 ]
          (%o8)                   [              ]
                                  [ B      B     ]
                                  [  2, 1   2, 2 ]


 -- Function: gramschmidt (<x>)
 -- Function: gschmit (<x>)
     Carries out the Gram-Schmidt orthogonalization algorithm on <x>,
     which is either a matrix or a list of lists.  <x> is not modified
     by `gramschmidt'.

     If <x> is a matrix, the algorithm is applied to the rows of <x>.
     If <x> is a list of lists, the algorithm is applied to the
     sublists, which must have equal numbers of elements.  In either
     case, the return value is a list of lists, the sublists of which
     are orthogonal and span the same space as <x>.  If the dimension
     of the span of <x> is less than the number of rows or sublists,
     some sublists of the return value are zero.

     `factor' is called at each stage of the algorithm to simplify
     intermediate results.  As a consequence, the return value may
     contain factored integers.

     `gschmit' (note spelling) is a synonym for `gramschmidt'.

     `load ("eigen")' loads this function.

     Example:

          (%i1) load ("eigen")$
          Warning - you are redefining the Macsyma function eigenvalues
          Warning - you are redefining the Macsyma function eigenvectors
          (%i2) x: matrix ([1, 2, 3], [9, 18, 30], [12, 48, 60]);
                                   [ 1   2   3  ]
                                   [            ]
          (%o2)                    [ 9   18  30 ]
                                   [            ]
                                   [ 12  48  60 ]
          (%i3) y: gramschmidt (x);
                                 2      2            4     3
                                3      3   3 5      2  3  2  3
          (%o3)  [[1, 2, 3], [- ---, - --, ---], [- ----, ----, 0]]
                                2 7    7   2 7       5     5
          (%i4) i: innerproduct$
          (%i5) [i (y[1], y[2]), i (y[2], y[3]), i (y[3], y[1])];
          (%o5)                       [0, 0, 0]


 -- Function: hach (<a>, <b>, <m>, <n>, <l>)
     `hach' is an implementation of Hacijan's linear programming
     algorithm.

     `load ("kach")' loads this function.  `demo ("kach")' executes a
     demonstration of this function.


 -- Function: ident (<n>)
     Returns an <n> by <n> identity matrix.


 -- Function: innerproduct (<x>, <y>)
 -- Function: inprod (<x>, <y>)
     Returns the inner product (also called the scalar product or dot
     product) of <x> and <y>, which are lists of equal length, or both
     1-column or 1-row matrices of equal length.  The return value is
     `conjugate (x) . y', where `.' is the noncommutative
     multiplication operator.

     `load ("eigen")' loads this function.

     `inprod' is a synonym for `innerproduct'.


 -- Function: invert (<M>)
     Returns the inverse of the matrix <M>.  The inverse is computed by
     the adjoint method.

     This allows a user to compute the inverse of a matrix with bfloat
     entries or polynomials with floating pt. coefficients without
     converting to cre-form.

     Cofactors are computed by the `determinant' function, so if
     `ratmx' is `false' the inverse is computed without changing the
     representation of the elements.

     The current implementation is inefficient for matrices of high
     order.

     When `detout' is `true', the determinant is factored out of the
     inverse.

     The elements of the inverse are not automatically expanded.  If
     <M> has polynomial elements, better appearing output can be
     generated by `expand (invert (m)), detout'.  If it is desirable to
     then divide through by the determinant this can be accomplished by
     `xthru (%)' or alternatively from scratch by

          expand (adjoint (m)) / expand (determinant (m))
          invert (m) := adjoint (m) / determinant (m)

     See `^^' (noncommutative exponent) for another method of inverting
     a matrix.


 -- Option variable: lmxchar
     Default value: `['

     `lmxchar' is the character displayed as the left delimiter of a
     matrix.  See also `rmxchar'.

     Example:

          (%i1) lmxchar: "|"$
          (%i2) matrix ([a, b, c], [d, e, f], [g, h, i]);
                                     | a  b  c ]
                                     |         ]
          (%o2)                      | d  e  f ]
                                     |         ]
                                     | g  h  i ]


 -- Function: matrix (<row_1>, ..., <row_n>)
     Returns a rectangular matrix which has the rows <row_1>, ...,
     <row_n>.  Each row is a list of expressions.  All rows must be the
     same length.

     The operations `+' (addition), `-' (subtraction), `*'
     (multiplication), and `/' (division), are carried out element by
     element when the operands are two matrices, a scalar and a matrix,
     or a matrix and a scalar.  The operation `^' (exponentiation,
     equivalently `**') is carried out element by element if the
     operands are a scalar and a matrix or a matrix and a scalar, but
     not if the operands are two matrices.  All operations are normally
     carried out in full, including `.' (noncommutative multiplication).

     Matrix multiplication is represented by the noncommutative
     multiplication operator `.'.  The corresponding noncommutative
     exponentiation operator is `^^'.  For a matrix `<A>', `<A>.<A> =
     <A>^^2' and `<A>^^-1' is the inverse of <A>, if it exists.

     There are switches for controlling simplification of expressions
     involving dot and matrix-list operations.  These are `doallmxops',
     `domxexpt' `domxmxops', `doscmxops', and `doscmxplus'.

     There are additional options which are related to matrices. These
     are: `lmxchar', `rmxchar', `ratmx', `listarith', `detout',
     `scalarmatrix', and `sparse'.

     There are a number of functions which take matrices as arguments
     or yield matrices as return values.  See `eigenvalues',
     `eigenvectors', `determinant', `charpoly', `genmatrix', `addcol',
     `addrow', `copymatrix', `transpose', `echelon', and `rank'.

     Examples:

        * Construction of matrices from lists.

          (%i1) x: matrix ([17, 3], [-8, 11]);
                                     [ 17   3  ]
          (%o1)                      [         ]
                                     [ - 8  11 ]
          (%i2) y: matrix ([%pi, %e], [a, b]);
                                     [ %pi  %e ]
          (%o2)                      [         ]
                                     [  a   b  ]

        * Addition, element by element.

          (%i3) x + y;
                                [ %pi + 17  %e + 3 ]
          (%o3)                 [                  ]
                                [  a - 8    b + 11 ]

        * Subtraction, element by element.

          (%i4) x - y;
                                [ 17 - %pi  3 - %e ]
          (%o4)                 [                  ]
                                [ - a - 8   11 - b ]

        * Multiplication, element by element.

          (%i5) x * y;
                                  [ 17 %pi  3 %e ]
          (%o5)                   [              ]
                                  [ - 8 a   11 b ]

        * Division, element by element.

          (%i6) x / y;
                                  [ 17       - 1 ]
                                  [ ---  3 %e    ]
                                  [ %pi          ]
          (%o6)                   [              ]
                                  [   8    11    ]
                                  [ - -    --    ]
                                  [   a    b     ]

        * Matrix to a scalar exponent, element by element.

          (%i7) x ^ 3;
                                   [ 4913    27  ]
          (%o7)                    [             ]
                                   [ - 512  1331 ]

        * Scalar base to a matrix exponent, element by element.

          (%i8) exp(y);
                                   [   %pi    %e ]
                                   [ %e     %e   ]
          (%o8)                    [             ]
                                   [    a     b  ]
                                   [  %e    %e   ]

        * Matrix base to a matrix exponent. This is not carried out
          element by element.

          (%i9) x ^ y;
                                          [ %pi  %e ]
                                          [         ]
                                          [  a   b  ]
                               [ 17   3  ]
          (%o9)                [         ]
                               [ - 8  11 ]

        * Noncommutative matrix multiplication.

          (%i10) x . y;
                            [ 3 a + 17 %pi  3 b + 17 %e ]
          (%o10)            [                           ]
                            [ 11 a - 8 %pi  11 b - 8 %e ]
          (%i11) y . x;
                          [ 17 %pi - 8 %e  3 %pi + 11 %e ]
          (%o11)          [                              ]
                          [  17 a - 8 b     11 b + 3 a   ]

        * Noncommutative matrix exponentiation.  A scalar base <b> to a
          matrix power <M> is carried out element by element and so
          `b^^m' is the same as `b^m'.

          (%i12) x ^^ 3;
                                  [  3833   1719 ]
          (%o12)                  [              ]
                                  [ - 4584  395  ]
          (%i13) %e ^^ y;
                                   [   %pi    %e ]
                                   [ %e     %e   ]
          (%o13)                   [             ]
                                   [    a     b  ]
                                   [  %e    %e   ]

        * A matrix raised to a -1 exponent with noncommutative
          exponentiation is the matrix inverse, if it exists.

          (%i14) x ^^ -1;
                                   [ 11      3  ]
                                   [ ---  - --- ]
                                   [ 211    211 ]
          (%o14)                   [            ]
                                   [  8    17   ]
                                   [ ---   ---  ]
                                   [ 211   211  ]
          (%i15) x . (x ^^ -1);
                                      [ 1  0 ]
          (%o15)                      [      ]
                                      [ 0  1 ]


 -- Function: matrixmap (<f>, <M>)
     Returns a matrix with element `i,j' equal to `<f>(<M>[i,j])'.

     See also `map', `fullmap', `fullmapl', and `apply'.


 -- Function: matrixp (<expr>)
     Returns `true' if <expr> is a matrix, otherwise `false'.


 -- Option variable: matrix_element_add
     Default value: `+'

     `matrix_element_add' is the operation invoked in place of addition
     in a matrix multiplication.  `matrix_element_add' can be assigned
     any n-ary operator (that is, a function which handles any number
     of arguments).  The assigned value may be the name of an operator
     enclosed in quote marks, the name of a function, or a lambda
     expression.

     See also `matrix_element_mult' and `matrix_element_transpose'.

     Example:

          (%i1) matrix_element_add: "*"$
          (%i2) matrix_element_mult: "^"$
          (%i3) aa: matrix ([a, b, c], [d, e, f]);
                                     [ a  b  c ]
          (%o3)                      [         ]
                                     [ d  e  f ]
          (%i4) bb: matrix ([u, v, w], [x, y, z]);
                                     [ u  v  w ]
          (%o4)                      [         ]
                                     [ x  y  z ]
          (%i5) aa . transpose (bb);
                               [  u  v  w   x  y  z ]
                               [ a  b  c   a  b  c  ]
          (%o5)                [                    ]
                               [  u  v  w   x  y  z ]
                               [ d  e  f   d  e  f  ]


 -- Option variable: matrix_element_mult
     Default value: `*'

     `matrix_element_mult' is the operation invoked in place of
     multiplication in a matrix multiplication.  `matrix_element_mult'
     can be assigned any binary operator.  The assigned value may be
     the name of an operator enclosed in quote marks, the name of a
     function, or a lambda expression.

     The dot operator `.' is a useful choice in some contexts.

     See also `matrix_element_add' and `matrix_element_transpose'.

     Example:

          (%i1) matrix_element_add: lambda ([[x]], sqrt (apply ("+", x)))$
          (%i2) matrix_element_mult: lambda ([x, y], (x - y)^2)$
          (%i3) [a, b, c] . [x, y, z];
                                    2          2          2
          (%o3)         sqrt((c - z)  + (b - y)  + (a - x) )
          (%i4) aa: matrix ([a, b, c], [d, e, f]);
                                     [ a  b  c ]
          (%o4)                      [         ]
                                     [ d  e  f ]
          (%i5) bb: matrix ([u, v, w], [x, y, z]);
                                     [ u  v  w ]
          (%o5)                      [         ]
                                     [ x  y  z ]
          (%i6) aa . transpose (bb);
                         [             2          2          2  ]
                         [ sqrt((c - w)  + (b - v)  + (a - u) ) ]
          (%o6)  Col 1 = [                                      ]
                         [             2          2          2  ]
                         [ sqrt((f - w)  + (e - v)  + (d - u) ) ]

                                   [             2          2          2  ]
                                   [ sqrt((c - z)  + (b - y)  + (a - x) ) ]
                           Col 2 = [                                      ]
                                   [             2          2          2  ]
                                   [ sqrt((f - z)  + (e - y)  + (d - x) ) ]


 -- Option variable: matrix_element_transpose
     Default value: `false'

     `matrix_element_transpose' is the operation applied to each
     element of a matrix when it is transposed.  `matrix_element_mult'
     can be assigned any unary operator.  The assigned value may be the
     name of an operator enclosed in quote marks, the name of a
     function, or a lambda expression.

     When `matrix_element_transpose' equals `transpose', the
     `transpose' function is applied to every element.  When
     `matrix_element_transpose' equals `nonscalars', the `transpose'
     function is applied to every nonscalar element.  If some element
     is an atom, the `nonscalars' option applies `transpose' only if
     the atom is declared nonscalar, while the `transpose' option
     always applies `transpose'.

     The default value, `false', means no operation is applied.

     See also `matrix_element_add' and `matrix_element_mult'.

     Examples:

          (%i1) declare (a, nonscalar)$
          (%i2) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o2)                   [              ]
                                  [      b       ]
          (%i3) matrix_element_transpose: nonscalars$
          (%i4) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o4)                   [              ]
                                  [      b       ]
          (%i5) matrix_element_transpose: transpose$
          (%i6) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o6)                   [              ]
                                  [ transpose(b) ]
          (%i7) matrix_element_transpose: lambda ([x], realpart(x) - %i*imagpart(x))$
          (%i8) m: matrix ([1 + 5*%i, 3 - 2*%i], [7*%i, 11]);
                               [ 5 %i + 1  3 - 2 %i ]
          (%o8)                [                    ]
                               [   7 %i       11    ]
          (%i9) transpose (m);
                                [ 1 - 5 %i  - 7 %i ]
          (%o9)                 [                  ]
                                [ 2 %i + 3    11   ]


 -- Function: mattrace (<M>)
     Returns the trace (that is, the sum of the elements on the main
     diagonal) of the square matrix <M>.

     `mattrace' is called by `ncharpoly', an alternative to Maxima's
     `charpoly'.

     `load ("nchrpl")' loads this function.


 -- Function: minor (<M>, <i>, <j>)
     Returns the <i>, <j> minor of the matrix <M>.  That is, <M> with
     row <i> and column <j> removed.


 -- Function: ncexpt (<a>, <b>)
     If a non-commutative exponential expression is too wide to be
     displayed as `<a>^^<b>' it appears as `ncexpt (<a>,<b>)'.

     `ncexpt' is not the name of a function or operator; the name only
     appears in output, and is not recognized in input.


 -- Function: ncharpoly (<M>, <x>)
     Returns the characteristic polynomial of the matrix <M> with
     respect to <x>.  This is an alternative to Maxima's `charpoly'.

     `ncharpoly' works by computing traces of powers of the given
     matrix, which are known to be equal to sums of powers of the roots
     of the characteristic polynomial.  From these quantities the
     symmetric functions of the roots can be calculated, which are
     nothing more than the coefficients of the characteristic
     polynomial.  `charpoly' works by forming the determinant of `<x> *
     ident [n] - a'.  Thus `ncharpoly' wins, for example, in the case
     of large dense matrices filled with integers, since it avoids
     polynomial arithmetic altogether.

     `load ("nchrpl")' loads this file.


 -- Function: newdet (<M>, <n>)
     Computes the determinant of the matrix or array <M> by the
     Johnson-Gentleman tree minor algorithm.  The argument <n> is the
     order; it is optional if <M> is a matrix.


 -- Declaration: nonscalar
     Makes atoms behave as does a list or matrix with respect to the
     dot operator.


 -- Function: nonscalarp (<expr>)
     Returns `true' if <expr> is a non-scalar, i.e., it contains atoms
     declared as non-scalars, lists, or matrices.


 -- Function: permanent (<M>, <n>)
     Computes the permanent of the matrix <M>.  A permanent is like a
     determinant but with no sign changes.


 -- Function: rank (<M>)
     Computes the rank of the matrix <M>.  That is, the order of the
     largest non-singular subdeterminant of <M>.

     <rank> may return the wrong answer if it cannot determine that a
     matrix element that is equivalent to zero is indeed so.


 -- Option variable: ratmx
     Default value: `false'

     When `ratmx' is `false', determinant and matrix addition,
     subtraction, and multiplication are performed in the
     representation of the matrix elements and cause the result of
     matrix inversion to be left in general representation.

     When `ratmx' is `true', the 4 operations mentioned above are
     performed in CRE form and the result of matrix inverse is in CRE
     form.  Note that this may cause the elements to be expanded
     (depending on the setting of `ratfac') which might not always be
     desired.


 -- Function: row (<M>, <i>)
     Returns the <i>'th row of the matrix <M>.  The return value is a
     matrix.


 -- Option variable: scalarmatrixp
     Default value: `true'

     When `scalarmatrixp' is `true', then whenever a 1 x 1 matrix is
     produced as a result of computing the dot product of matrices it
     is simplified to a scalar, namely the sole element of the matrix.

     When `scalarmatrixp' is `all', then all 1 x 1 matrices are
     simplified to scalars.

     When `scalarmatrixp' is `false', 1 x 1 matrices are not simplified
     to scalars.


 -- Function: scalefactors (<coordinatetransform>)
     Here coordinatetransform evaluates to the form [[expression1,
     expression2, ...], indeterminate1, indeterminat2, ...], where
     indeterminate1, indeterminate2, etc. are the curvilinear
     coordinate variables and where a set of rectangular Cartesian
     components is given in terms of the curvilinear coordinates by
     [expression1, expression2, ...].  `coordinates' is set to the
     vector [indeterminate1, indeterminate2,...], and `dimension' is
     set to the length of this vector.  SF[1], SF[2], ...,
     SF[DIMENSION] are set to the coordinate scale factors, and `sfprod'
     is set to the product of these scale factors.  Initially,
     `coordinates' is [X, Y, Z], `dimension' is 3, and
     SF[1]=SF[2]=SF[3]=SFPROD=1, corresponding to 3-dimensional
     rectangular Cartesian coordinates.  To expand an expression into
     physical components in the current coordinate system, there is a
     function with usage of the form


 -- Function: setelmx (<x>, <i>, <j>, <M>)
     Assigns <x> to the (<i>, <j>)'th element of the matrix <M>, and
     returns the altered matrix.

     `<M> [<i>, <j>]: <x>' has the same effect, but returns <x> instead
     of <M>.


 -- Function: similaritytransform (<M>)
 -- Function: simtran (<M>)
     `similaritytransform' computes a similarity transform of the
     matrix `M'.  It returns a list which is the output of the
     `uniteigenvectors' command.  In addition if the flag
     `nondiagonalizable' is `false' two global matrices `leftmatrix'
     and `rightmatrix' are computed.  These matrices have the property
     that `leftmatrix . <M> . rightmatrix' is a diagonal matrix with
     the eigenvalues of <M> on the diagonal.  If `nondiagonalizable' is
     `true' the left and right matrices are not computed.

     If the flag `hermitianmatrix' is `true' then `leftmatrix' is the
     complex conjugate of the transpose of `rightmatrix'.  Otherwise
     `leftmatrix' is the inverse of `rightmatrix'.

     `rightmatrix' is the matrix the columns of which are the unit
     eigenvectors of <M>.  The other flags (see `eigenvalues' and
     `eigenvectors') have the same effects since `similaritytransform'
     calls the other functions in the package in order to be able to
     form `rightmatrix'.

     `load ("eigen")' loads this function.

     `simtran' is a synonym for `similaritytransform'.


 -- Option variable: sparse
     Default value: `false'

     When `sparse' is `true', and if `ratmx' is `true', then
     `determinant' will use special routines for computing sparse
     determinants.


 -- Function: submatrix (<i_1>, ..., <i_m>, <M>, <j_1>, ..., <j_n>)
 -- Function: submatrix (<i_1>, ..., <i_m>, <M>)
 -- Function: submatrix (<M>, <j_1>, ..., <j_n>)
     Returns a new matrix composed of the matrix <M> with rows <i_1>,
     ..., <i_m> deleted, and columns <j_1>, ..., <j_n> deleted.


 -- Function: transpose (<M>)
     Returns the transpose of <M>.

     If <M> is a matrix, the return value is another matrix <N> such
     that `N[i,j] = M[j,i]'.

     Otherwise <M> is a list, and the return value is a matrix <N> of
     `length (m)' rows and 1 column, such that `N[i,1] = M[i]'.


 -- Function: triangularize (<M>)
     Returns the upper triangular form of the matrix `M', as produced
     by Gaussian elimination.  The return value is the same as
     `echelon', except that the leading nonzero coefficient in each row
     is not normalized to 1.

     `lu_factor' and `cholesky' are other functions which yield
     triangularized matrices.

          (%i1) M: matrix ([3, 7, aa, bb], [-1, 8, 5, 2], [9, 2, 11, 4]);
                                 [  3   7  aa  bb ]
                                 [                ]
          (%o1)                  [ - 1  8  5   2  ]
                                 [                ]
                                 [  9   2  11  4  ]
          (%i2) triangularize (M);
                       [ - 1   8         5            2      ]
                       [                                     ]
          (%o2)        [  0   - 74     - 56         - 22     ]
                       [                                     ]
                       [  0    0    626 - 74 aa  238 - 74 bb ]


 -- Function: uniteigenvectors (<M>)
 -- Function: ueivects (<M>)
     Computes unit eigenvectors of the matrix <M>.  The return value is
     a list of lists, the first sublist of which is the output of the
     `eigenvalues' command, and the other sublists of which are the
     unit eigenvectors of the matrix corresponding to those eigenvalues
     respectively.

     The flags mentioned in the description of the `eigenvectors'
     command have the same effects in this one as well.

     When `knowneigvects' is `true', the `eigen' package assumes that
     the eigenvectors of the matrix are known to the user and are
     stored under the global name `listeigvects'.  `listeigvects'
     should be set to a list similar to the output of the
     `eigenvectors' command.

     If `knowneigvects' is set to `true' and the list of eigenvectors
     is given the setting of the flag `nondiagonalizable' may not be
     correct.  If that is the case please set it to the correct value.
     The author assumes that the user knows what he is doing and will
     not try to diagonalize a matrix the eigenvectors of which do not
     span the vector space of the appropriate dimension.

     `load ("eigen")' loads this function.

     `ueivects' is a synonym for `uniteigenvectors'.


 -- Function: unitvector (<x>)
 -- Function: uvect (<x>)
     Returns <x>/norm(<x>); this is a unit vector in the same direction
     as <x>.

     `load ("eigen")' loads this function.

     `uvect' is a synonym for `unitvector'.


 -- Function: vectorsimp (<expr>)
     Applies simplifications and expansions according to the following
     global flags:

     `expandall', `expanddot', `expanddotplus', `expandcross',
     `expandcrossplus', `expandcrosscross', `expandgrad',
     `expandgradplus', `expandgradprod', `expanddiv', `expanddivplus',
     `expanddivprod', `expandcurl', `expandcurlplus', `expandcurlcurl',
     `expandlaplacian', `expandlaplacianplus', and
     `expandlaplacianprod'.

     All these flags have default value `false'. The `plus' suffix
     refers to employing additivity or distributivity.  The `prod'
     suffix refers to the expansion for an operand that is any kind of
     product.

    `expandcrosscross'
          Simplifies p ~ (q ~ r) to (p . r)*q - (p . q)*r.

    `expandcurlcurl'
          Simplifies curl curl p to grad div p + div grad p.

    `expandlaplaciantodivgrad'
          Simplifies laplacian p to div grad p.

    `expandcross'
          Enables `expandcrossplus' and `expandcrosscross'.

    `expandplus'
          Enables `expanddotplus', `expandcrossplus', `expandgradplus',
          `expanddivplus', `expandcurlplus', and `expandlaplacianplus'.

    `expandprod'
          Enables `expandgradprod', `expanddivprod', and
          `expandlaplacianprod'.

     These flags have all been declared `evflag'.


 -- Option variable: vect_cross
     Default value: `false'

     When `vect_cross' is `true', it allows DIFF(X~Y,T) to work where ~
     is defined in SHARE;VECT (where VECT_CROSS is set to `true',
     anyway.)


 -- Function: zeromatrix (<m>, <n>)
     Returns an <m> by <n> matrix, all elements of which are zero.


 -- Special symbol: [
 -- Special symbol: ]
     `[' and `]' mark the beginning and end, respectively, of a list.

     `[' and `]' also enclose the subscripts of a list, array, hash
     array, or array function.

     Examples:

          (%i1) x: [a, b, c];
          (%o1)                       [a, b, c]
          (%i2) x[3];
          (%o2)                           c
          (%i3) array (y, fixnum, 3);
          (%o3)                           y
          (%i4) y[2]: %pi;
          (%o4)                          %pi
          (%i5) y[2];
          (%o5)                          %pi
          (%i6) z['foo]: 'bar;
          (%o6)                          bar
          (%i7) z['foo];
          (%o7)                          bar
          (%i8) g[k] := 1/(k^2+1);
                                            1
          (%o8)                     g  := ------
                                     k     2
                                          k  + 1
          (%i9) g[10];
                                          1
          (%o9)                          ---
                                         101



File: maxima.info,  Node: Affine,  Next: itensor,  Prev: Matrices and Linear Algebra,  Up: Top

27 Affine
*********

* Menu:

* Definitions for Affine::


File: maxima.info,  Node: Definitions for Affine,  Prev: Affine,  Up: Affine

27.1 Definitions for Affine
===========================

 -- Function: fast_linsolve ([<expr_1>, ..., <expr_m>], [<x_1>, ...,
          <x_n>])
     Solves the simultaneous linear equations <expr_1>, ..., <expr_m>
     for the variables <x_1>, ..., <x_n>.  Each <expr_i> may be an
     equation or a general expression; if given as a general
     expression, it is treated as an equation of the form `<expr_i> =
     0'.

     The return value is a list of equations of the form `[<x_1> =
     <a_1>, ..., <x_n> = <a_n>]' where <a_1>, ..., <a_n> are all free
     of <x_1>, ..., <x_n>.

     `fast_linsolve' is faster than `linsolve' for system of equations
     which are sparse.


 -- Function: grobner_basis ([<expr_1>, ..., <expr_m>])
     Returns a Groebner basis for the equations <expr_1>, ..., <expr_m>.
     The function `polysimp' can then be used to simplify other
     functions relative to the equations.

          grobner_basis ([3*x^2+1, y*x])$

          polysimp (y^2*x + x^3*9 + 2) ==> -3*x + 2

     `polysimp(f)' yields 0 if and only if <f> is in the ideal
     generated by <expr_1>, ..., <expr_m>, that is, if and only if <f>
     is a polynomial combination of the elements of <expr_1>, ...,
     <expr_m>.


 -- Function: set_up_dot_simplifications (<eqns>,
          <check_through_degree>)
 -- Function: set_up_dot_simplifications (<eqns>)
     The <eqns> are polynomial equations in non commutative variables.
     The value of `current_variables' is the list of variables used for
     computing degrees.  The equations must be homogeneous, in order
     for the procedure to terminate.

     If you have checked overlapping simplifications in
     `dot_simplifications' above the degree of <f>, then the following
     is true: `dotsimp (<f>)' yields 0 if and only if <f> is in the
     ideal generated by the equations, i.e., if and only if <f> is a
     polynomial combination of the elements of the equations.

     The degree is that returned by `nc_degree'.   This in turn is
     influenced by the weights of individual variables.


 -- Function: declare_weight (<x_1>, <w_1>, ..., <x_n>, <w_n>)
     Assigns weights <w_1>, ..., <w_n> to <x_1>, ..., <x_n>,
     respectively.  These are the weights used in computing `nc_degree'.


 -- Function: nc_degree (<p>)
     Returns the degree of a noncommutative polynomial <p>.  See
     `declare_weights'.


 -- Function: dotsimp (<f>)
     Returns 0 if and only if <f> is in the ideal generated by the
     equations, i.e., if and only if <f> is a polynomial combination of
     the elements of the equations.


 -- Function: fast_central_elements ([<x_1>, ..., <x_n>], <n>)
     If `set_up_dot_simplifications' has been previously done, finds
     the central polynomials in the variables <x_1>, ..., <x_n> in the
     given degree, <n>.

     For example:
          set_up_dot_simplifications ([y.x + x.y], 3);
          fast_central_elements ([x, y], 2);
          [y.y, x.x];


 -- Function: check_overlaps (<n>, <add_to_simps>)
     Checks the overlaps thru degree <n>, making sure that you have
     sufficient simplification rules in each degree, for `dotsimp' to
     work correctly.  This process can be speeded up if you know before
     hand what the dimension of the space of monomials is.  If it is of
     finite global dimension, then `hilbert' should be used.  If you
     don't know the monomial dimensions, do not specify a
     `rank_function'.  An optional third argument `reset', `false' says
     don't bother to query about resetting things.


 -- Function: mono ([<x_1>, ..., <x_n>], <n>)
     Returns the list of independent monomials relative to the current
     dot simplifications of degree <n> in the variables <x_1>, ...,
     <x_n>.


 -- Function: monomial_dimensions (<n>)
     Compute the Hilbert series through degree <n> for the current
     algebra.


 -- Function: extract_linear_equations ([<p_1>, ..., <p_n>], [<m_1>,
          ..., <m_n>])
     Makes a list of the coefficients of the noncommutative polynomials
     <p_1>, ..., <p_n> of the noncommutative monomials <m_1>, ...,
     <m_n>.  The coefficients should be scalars.   Use
     `list_nc_monomials' to build the list of monomials.


 -- Function: list_nc_monomials ([<p_1>, ..., <p_n>])
 -- Function: list_nc_monomials (<p>)
     Returns a list of the non commutative monomials occurring in a
     polynomial <p> or a list of polynomials <p_1>, ..., <p_n>.


 -- Option variable: all_dotsimp_denoms
     Default value: `false'

     When `all_dotsimp_denoms' is a list, the denominators encountered
     by `dotsimp' are appended to the list.  `all_dotsimp_denoms' may
     be initialized to an empty list `[]' before calling `dotsimp'.

     By default, denominators are not collected by `dotsimp'.



File: maxima.info,  Node: itensor,  Next: ctensor,  Prev: Affine,  Up: Top

28 itensor
**********

* Menu:

* Introduction to itensor::
* Definitions for itensor::


File: maxima.info,  Node: Introduction to itensor,  Next: Definitions for itensor,  Prev: itensor,  Up: itensor

28.1 Introduction to itensor
============================

Maxima implements symbolic tensor manipulation of two distinct types:
component tensor manipulation (`ctensor' package) and indicial tensor
manipulation (`itensor' package).

   Nota bene: Please see the note on 'new tensor notation' below.

   Component tensor manipulation means that geometrical tensor objects
are represented as arrays or matrices. Tensor operations such as
contraction or covariant differentiation are carried out by actually
summing over repeated (dummy) indices with `do' statements.  That is,
one explicitly performs operations on the appropriate tensor components
stored in an array or matrix.

   Indicial tensor manipulation is implemented by representing tensors
as functions of their covariant, contravariant and derivative indices.
Tensor operations such as contraction or covariant differentiation are
performed by manipulating the indices themselves rather than the
components to which they correspond.

   These two approaches to the treatment of differential, algebraic and
analytic processes in the context of Riemannian geometry have various
advantages and disadvantages which reveal themselves only through the
particular nature and difficulty of the user's problem.  However, one
should keep in mind the following characteristics of the two
implementations:

   The representation of tensors and tensor operations explicitly in
terms of their components makes `ctensor' easy to use. Specification of
the metric and the computation of the induced tensors and invariants is
straightforward. Although all of Maxima's powerful simplification
capacity is at hand, a complex metric with intricate functional and
coordinate dependencies can easily lead to expressions whose size is
excessive and whose structure is hidden. In addition, many calculations
involve intermediate expressions which swell causing programs to
terminate before completion. Through experience, a user can avoid many
of these difficulties.

   Because of the special way in which tensors and tensor operations
are represented in terms of symbolic operations on their indices,
expressions which in the component representation would be unmanageable
can sometimes be greatly simplified by using the special routines for
symmetrical objects in `itensor'. In this way the structure of a large
expression may be more transparent. On the other hand, because of the
the special indicial representation in `itensor', in some cases the
user may find difficulty with the specification of the metric, function
definition, and the evaluation of differentiated "indexed" objects.

28.1.1 New tensor notation
--------------------------

Until now, the `itensor' package in Maxima has used a notation that
sometimes led to incorrect index ordering. Consider the following, for
instance:

     (%i2) imetric(g);
     (%o2)                                done
     (%i3) ishow(g([],[j,k])*g([],[i,l])*a([i,j],[]))$
                                      i l  j k
     (%t3)                           g    g    a
                                                i j
     (%i4) ishow(contract(%))$
                                           k l
     (%t4)                                a

   This result is incorrect unless `a' happens to be a symmetric tensor.
The reason why this happens is that although `itensor' correctly
maintains the order within the set of covariant and contravariant
indices, once an index is raised or lowered, its position relative to
the other set of indices is lost.

   To avoid this problem, a new notation has been developed that
remains fully compatible with the existing notation and can be used
interchangeably. In this notation, contravariant indices are inserted
in the appropriate positions in the covariant index list, but with a
minus sign prepended.  Functions like `contract' and `ishow' are now
aware of this new index notation and can process tensors appropriately.

   In this new notation, the previous example yields a correct result:

     (%i5) ishow(g([-j,-k],[])*g([-i,-l],[])*a([i,j],[]))$
                                      i l       j k
     (%t5)                           g    a    g
                                           i j
     (%i6) ishow(contract(%))$
                                           l k
     (%t6)                                a

   Presently, the only code that makes use of this notation is the
`lc2kdt' function. Through this notation, it achieves consistent
results as it applies the metric tensor to resolve Levi-Civita symbols
without resorting to numeric indices.

   Since this code is brand new, it probably contains bugs. While it
has been tested to make sure that it doesn't break anything using the
"old" tensor notation, there is a considerable chance that "new"
tensors will fail to interoperate with certain functions or features.
These bugs will be fixed as they are encountered... until then, caveat
emptor!

28.1.2 Indicial tensor manipulation
-----------------------------------

The indicial tensor manipulation package may be loaded by
`load(itensor)'. Demos are also available: try `demo(tensor)'.

   In `itensor' a tensor is represented as an "indexed object" .  This
is a function of 3 groups of indices which represent the covariant,
contravariant and derivative indices.  The covariant indices are
specified by a list as the first argument to the indexed object, and
the contravariant indices by a list as the second argument. If the
indexed object lacks either of these groups of indices then the empty
list `[]' is given as the corresponding argument.  Thus, `g([a,b],[c])'
represents an indexed object called `g' which has two covariant indices
`(a,b)', one contravariant index (`c') and no derivative indices.

   The derivative indices, if they are present, are appended as
additional arguments to the symbolic function representing the tensor.
They can be explicitly specified by the user or be created in the
process of differentiation with respect to some coordinate variable.
Since ordinary differentiation is commutative, the derivative indices
are sorted alphanumerically, unless `iframe_flag' is set to `true',
indicating that a frame metric is being used. This canonical ordering
makes it possible for Maxima to recognize that, for example,
`t([a],[b],i,j)' is the same as `t([a],[b],j,i)'.  Differentiation of
an indexed object with respect to some coordinate whose index does not
appear as an argument to the indexed object would normally yield zero.
This is because Maxima would not know that the tensor represented by
the indexed object might depend implicitly on the corresponding
coordinate.  By modifying the existing Maxima function `diff' in
`itensor', Maxima now assumes that all indexed objects depend on any
variable of differentiation unless otherwise stated.  This makes it
possible for the summation convention to be extended to derivative
indices. It should be noted that `itensor' does not possess the
capabilities of raising derivative indices, and so they are always
treated as covariant.

   The following functions are available in the tensor package for
manipulating indexed objects.  At present, with respect to the
simplification routines, it is assumed that indexed objects do not by
default possess symmetry properties. This can be overridden by setting
the variable `allsym[false]' to `true', which will result in treating
all indexed objects completely symmetric in their lists of covariant
indices and symmetric in their lists of contravariant indices.

   The `itensor' package generally treats tensors as opaque objects.
Tensorial equations are manipulated based on algebraic rules,
specifically symmetry and contraction rules. In addition, the `itensor'
package understands covariant differentiation, curvature, and torsion.
Calculations can be performed relative to a metric of moving frame,
depending on the setting of the `iframe_flag' variable.

   A sample session below demonstrates how to load the `itensor'
package, specify the name of the metric, and perform some simple
calculations.

     (%i1) load(itensor);
     (%o1)      /share/tensor/itensor.lisp
     (%i2) imetric(g);
     (%o2)                                done
     (%i3) components(g([i,j],[]),p([i,j],[])*e([],[]))$
     (%i4) ishow(g([k,l],[]))$
     (%t4)                               e p
                                            k l
     (%i5) ishow(diff(v([i],[]),t))$
     (%t5)                                  0
     (%i6) depends(v,t);
     (%o6)                               [v(t)]
     (%i7) ishow(diff(v([i],[]),t))$
                                         d
     (%t7)                               -- (v )
                                         dt   i
     (%i8) ishow(idiff(v([i],[]),j))$
     (%t8)                                v
                                           i,j
     (%i9) ishow(extdiff(v([i],[]),j))$
     (%t9)                             v    - v
                                        j,i    i,j
                                       -----------
                                            2
     (%i10) ishow(liediff(v,w([i],[])))$
                                    %3          %3
     (%t10)                        v   w     + v   w
                                        i,%3    ,i  %3
     (%i11) ishow(covdiff(v([i],[]),j))$
                                                   %4
     (%t11)                        v    - v   ichr2
                                    i,j    %4      i j
     (%i12) ishow(ev(%,ichr2))$
                    %4 %5
     (%t12) v    - g      v   (e p       + e   p     - e p       - e    p
             i,j           %4     j %5,i    ,i  j %5      i j,%5    ,%5  i j

                                                     + e p       + e   p    )/2
                                                          i %5,j    ,j  i %5
     (%i13) iframe_flag:true;
     (%o13)                               true
     (%i14) ishow(covdiff(v([i],[]),j))$
                                                  %6
     (%t14)                        v    - v   icc2
                                    i,j    %6     i j
     (%i15) ishow(ev(%,icc2))$
                                                  %6
     (%t15)                        v    - v   ifc2
                                    i,j    %6     i j
     (%i16) ishow(radcan(ev(%,ifc2,ifc1)))$
                  %6 %8                    %6 %8
     (%t16) - (ifg      v   ifb       + ifg      v   ifb       - 2 v
                         %6    j %8 i             %6    i j %8      i,j

                                                         %6 %8
                                                    - ifg      v   ifb      )/2
                                                                %6    %8 i j
     (%i17) ishow(canform(s([i,j],[])-s([j,i])))$
     (%t17)                            s    - s
                                        i j    j i
     (%i18) decsym(s,2,0,[sym(all)],[]);
     (%o18)                               done
     (%i19) ishow(canform(s([i,j],[])-s([j,i])))$
     (%t19)                                 0
     (%i20) ishow(canform(a([i,j],[])+a([j,i])))$
     (%t20)                            a    + a
                                        j i    i j
     (%i21) decsym(a,2,0,[anti(all)],[]);
     (%o21)                               done
     (%i22) ishow(canform(a([i,j],[])+a([j,i])))$
     (%t22)                                 0


File: maxima.info,  Node: Definitions for itensor,  Prev: Introduction to itensor,  Up: itensor

28.2 Definitions for itensor
============================

28.2.1 Managing indexed objects
-------------------------------

 -- Function: entertensor (<name>)
     is a function which, by prompting, allows one to create an indexed
     object called <name> with any number of tensorial and derivative
     indices. Either a single index or a list of indices (which may be
     null) is acceptable input (see the example under `covdiff').


 -- Function: changename (<old>, <new>, <expr>)
     will change the name of all indexed objects called <old> to <new>
     in <expr>. <old> may be either a symbol or a list of the form
     `[<name>, <m>, <n>]' in which case only those indexed objects
     called <name> with <m> covariant and <n> contravariant indices
     will be renamed to <new>.


 -- Function: listoftens
     Lists all tensors in a tensorial expression, complete with their
     indices. E.g.,


          (%i6) ishow(a([i,j],[k])*b([u],[],v)+c([x,y],[])*d([],[])*e)$
                                                   k
          (%t6)                        d e c    + a    b
                                            x y    i j  u,v
          (%i7) ishow(listoftens(%))$
                                         k
          (%t7)                        [a   , b   , c   , d]
                                         i j   u,v   x y


 -- Function: ishow (<expr>)
     displays <expr> with the indexed objects in it shown having their
     covariant indices as subscripts and contravariant indices as
     superscripts. The derivative indices are displayed as subscripts,
     separated from the covariant indices by a comma (see the examples
     throughout this document).


 -- Function: indices (<expr>)
     Returns a list of two elements.  The first is a list of the free
     indices in <expr> (those that occur only once). The second is the
     list of the dummy indices in <expr> (those that occur exactly
     twice) as the following example demonstrates.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(a([i,j],[k,l],m,n)*b([k,o],[j,m,p],q,r))$
                                          k l      j m p
          (%t2)                          a        b
                                          i j,m n  k o,q r
          (%i3) indices(%);
          (%o3)                 [[l, p, i, n, o, q, r], [k, j, m]]

     A tensor product containing the same index more than twice is
     syntactically illegal. `indices' attempts to deal with these
     expressions in a reasonable manner; however, when it is called to
     operate upon such an illegal expression, its behavior should be
     considered undefined.


 -- Function: rename (<expr>)
 -- Function: rename (<expr>, <count>)
     Returns an expression equivalent to <expr> but with the dummy
     indices in each term chosen from the set `[%1, %2,...]', if the
     optional second argument is omitted. Otherwise, the dummy indices
     are indexed beginning at the value of <count>.  Each dummy index
     in a product will be different. For a sum, `rename' will operate
     upon each term in the sum resetting the counter with each term. In
     this way `rename' can serve as a tensorial simplifier. In
     addition, the indices will be sorted alphanumerically (if `allsym'
     is `true') with respect to covariant or contravariant indices
     depending upon the value of `flipflag'.  If `flipflag' is `false'
     then the indices will be renamed according to the order of the
     contravariant indices. If `flipflag' is `true' the renaming will
     occur according to the order of the covariant indices. It often
     happens that the combined effect of the two renamings will reduce
     an expression more than either one by itself.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) allsym:true;
          (%o2)                                true
          (%i3) g([],[%4,%5])*g([],[%6,%7])*ichr2([%1,%4],[%3])*
          ichr2([%2,%3],[u])*ichr2([%5,%6],[%1])*ichr2([%7,r],[%2])-
          g([],[%4,%5])*g([],[%6,%7])*ichr2([%1,%2],[u])*
          ichr2([%3,%5],[%1])*ichr2([%4,%6],[%3])*ichr2([%7,r],[%2]),noeval$
          (%i4) expr:ishow(%)$

                 %4 %5  %6 %7      %3         u          %1         %2
          (%t4) g      g      ichr2      ichr2      ichr2      ichr2
                                   %1 %4      %2 %3      %5 %6      %7 r

                        %4 %5  %6 %7      u          %1         %3         %2
                     - g      g      ichr2      ichr2      ichr2      ichr2
                                          %1 %2      %3 %5      %4 %6      %7 r
          (%i5) flipflag:true;
          (%o5)                                true
          (%i6) ishow(rename(expr))$
                 %2 %5  %6 %7      %4         u          %1         %3
          (%t6) g      g      ichr2      ichr2      ichr2      ichr2
                                   %1 %2      %3 %4      %5 %6      %7 r

                        %4 %5  %6 %7      u          %1         %3         %2
                     - g      g      ichr2      ichr2      ichr2      ichr2
                                          %1 %2      %3 %4      %5 %6      %7 r
          (%i7) flipflag:false;
          (%o7)                                false
          (%i8) rename(%th(2));
          (%o8)                                  0
          (%i9) ishow(rename(expr))$
                 %1 %2  %3 %4      %5         %6         %7        u
          (%t9) g      g      ichr2      ichr2      ichr2     ichr2
                                   %1 %6      %2 %3      %4 r      %5 %7

                        %1 %2  %3 %4      %6         %5         %7        u
                     - g      g      ichr2      ichr2      ichr2     ichr2
                                          %1 %3      %2 %6      %4 r      %5 %7


 -- Option variable: flipflag
     Default: `false'. If `false' then the indices will be renamed
     according to the order of the contravariant indices, otherwise
     according to the order of the covariant indices.

     If `flipflag' is `false' then `rename' forms a list of the
     contravariant indices as they are encountered from left to right
     (if `true' then of the covariant indices). The first dummy index
     in the list is renamed to `%1', the next to `%2', etc.  Then
     sorting occurs after the `rename'-ing (see the example under
     `rename').


 -- Function: defcon (<tensor_1>)
 -- Function: defcon (<tensor_1>, <tensor_2>, <tensor_3>)
     gives <tensor_1> the property that the contraction of a product of
     <tensor_1> and <tensor_2> results in <tensor_3> with the
     appropriate indices.  If only one argument, <tensor_1>, is given,
     then the contraction of the product of <tensor_1> with any indexed
     object having the appropriate indices (say `my_tensor') will yield
     an indexed object with that name, i.e. `my_tensor', and with a new
     set of indices reflecting the contractions performed.      For
     example, if `imetric:g', then `defcon(g)' will implement the
     raising and lowering of indices through contraction with the metric
     tensor.      More than one `defcon' can be given for the same
     indexed object; the latest one given which applies in a particular
     contraction will be used.  `contractions' is a list of those
     indexed objects which have been given contraction properties with
     `defcon'.


 -- Function: remcon (<tensor_1>, ..., <tensor_n>)
 -- Function: remcon (all)
     removes all the contraction properties from the <tensor_1>, ...,
     <tensor_n>). `remcon(all)' removes all contraction properties from
     all indexed objects.


 -- Function: contract (<expr>)
     Carries out the tensorial contractions in <expr> which may be any
     combination of sums and products. This function uses the
     information given to the `defcon' function. For best results,
     `expr' should be fully expanded. `ratexpand' is the fastest way to
     expand products and powers of sums if there are no variables in
     the denominators of the terms. The `gcd' switch should be `false'
     if GCD cancellations are unnecessary.


 -- Function: indexed_tensor (<tensor>)
     Must be executed before assigning components to a <tensor> for
     which a built in value already exists as with `ichr1', `ichr2',
     `icurvature'. See the example under `icurvature'.


 -- Function: components (<tensor>, <expr>)
     permits one to assign an indicial value to an expression <expr>
     giving the values of the components of <tensor>. These are
     automatically substituted for the tensor whenever it occurs with
     all of its indices. The tensor must be of the form `t([...],[...])'
     where either list may be empty. <expr> can be any indexed
     expression involving other objects with the same free indices as
     <tensor>. When used to assign values to the metric tensor wherein
     the components contain dummy indices one must be careful to define
     these indices to avoid the generation of multiple dummy indices.
     Removal of this assignment is given to the function `remcomps'.

     It is important to keep in mind that `components' cares only about
     the valence of a tensor, not about any particular index ordering.
     Thus assigning components to, say, `x([i,-j],[])', `x([-j,i],[])',
     or `x([i],[j])' all produce the same result, namely components
     being assigned to a tensor named `x' with valence `(1,1)'.

     Components can be assigned to an indexed expression in four ways,
     two of which involve the use of the `components' command:

     1) As an indexed expression. For instance:


          (%i2) components(g([],[i,j]),e([],[i])*p([],[j]))$
          (%i3) ishow(g([],[i,j]))$
                                                i  j
          (%t3)                                e  p

     2) As a matrix:


          (%i6) components(g([i,j],[]),lg);
          (%o6)                                done
          (%i7) ishow(g([i,j],[]))$
          (%t7)                                g
                                                i j
          (%i8) g([3,3],[]);
          (%o8)                                  1
          (%i9) g([4,4],[]);
          (%o9)                                 - 1

     3) As a function. You can use a Maxima function to specify the
     components of a tensor based on its indices. For instance, the
     following code assigns `kdelta' to `h' if `h' has the same number
     of covariant and contravariant indices and no derivative indices,
     and `g' otherwise:


          (%i4) h(l1,l2,[l3]):=if length(l1)=length(l2) and length(l3)=0
            then kdelta(l1,l2) else apply(g,append([l1,l2], l3))$
          (%i5) ishow(h([i],[j]))$
                                                    j
          (%t5)                               kdelta
                                                    i
          (%i6) ishow(h([i,j],[k],l))$
                                               k
          (%t6)                               g
                                               i j,l

     4) Using Maxima's pattern matching capabilities, specifically the
     `defrule' and `applyb1' commands:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) matchdeclare(l1,listp);
          (%o2)                                done
          (%i3) defrule(r1,m(l1,[]),(i1:idummy(),
                g([l1[1],l1[2]],[])*q([i1],[])*e([],[i1])))$

          (%i4) defrule(r2,m([],l1),(i1:idummy(),
                w([],[l1[1],l1[2]])*e([i1],[])*q([],[i1])))$

          (%i5) ishow(m([i,n],[])*m([],[i,m]))$
                                              i m
          (%t5)                              m    m
                                                   i n
          (%i6) ishow(rename(applyb1(%,r1,r2)))$
                                     %1  %2  %3 m
          (%t6)                     e   q   w     q   e   g
                                                   %1  %2  %3 n


 -- Function: remcomps (<tensor>)
     Unbinds all values from <tensor> which were assigned with the
     `components' function.


 -- Function: showcomps (<tensor>)
     Shows component assignments of a tensor, as made using the
     `components' command. This function can be particularly useful
     when a matrix is assigned to an indicial tensor using
     `components', as demonstrated by the following example:


          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) load(itensor);
          (%o2)      /share/tensor/itensor.lisp
          (%i3) lg:matrix([sqrt(r/(r-2*m)),0,0,0],[0,r,0,0],
                          [0,0,sin(theta)*r,0],[0,0,0,sqrt((r-2*m)/r)]);
                         [         r                                     ]
                         [ sqrt(-------)  0       0              0       ]
                         [      r - 2 m                                  ]
                         [                                               ]
                         [       0        r       0              0       ]
          (%o3)          [                                               ]
                         [       0        0  r sin(theta)        0       ]
                         [                                               ]
                         [                                      r - 2 m  ]
                         [       0        0       0        sqrt(-------) ]
                         [                                         r     ]
          (%i4) components(g([i,j],[]),lg);
          (%o4)                                done
          (%i5) showcomps(g([i,j],[]));
                            [         r                                     ]
                            [ sqrt(-------)  0       0              0       ]
                            [      r - 2 m                                  ]
                            [                                               ]
                            [       0        r       0              0       ]
          (%t5)      g    = [                                               ]
                      i j   [       0        0  r sin(theta)        0       ]
                            [                                               ]
                            [                                      r - 2 m  ]
                            [       0        0       0        sqrt(-------) ]
                            [                                         r     ]
          (%o5)                                false

     The `showcomps' command can also display components of a tensor of
     rank higher than 2.


 -- Function: idummy ()
     Increments `icounter' and returns as its value an index of the form
     `%n' where n is a positive integer.  This guarantees that dummy
     indices which are needed in forming expressions will not conflict
     with indices already in use (see the example under `indices').


 -- Option variable: idummyx
     Default value: `%'

     Is the prefix for dummy indices (see the example under `indices').


 -- Option variable: icounter
     Default value: `1'

     Determines the numerical suffix to be used in generating the next
     dummy index in the tensor package.  The prefix is determined by
     the option `idummy' (default: `%').

 -- Function: kdelta (<L1>, <L2>)
     is the generalized Kronecker delta function defined in the
     `itensor' package with <L1> the list of covariant indices and <L2>
     the list of contravariant indices.  `kdelta([i],[j])' returns the
     ordinary Kronecker delta.  The command `ev(<expr>,kdelta)' causes
     the evaluation of an expression containing `kdelta([],[])' to the
     dimension of the manifold.

     In what amounts to an abuse of this notation, `itensor' also allows
     `kdelta' to have 2 covariant and no contravariant, or 2
     contravariant and no covariant indices, in effect providing a
     co(ntra)variant "unit matrix" capability. This is strictly
     considered a programming aid and not meant to imply that
     `kdelta([i,j],[])' is a valid tensorial object.


 -- Function: kdels (<L1>, <L2>)
     Symmetricized Kronecker delta, used in some calculations. For
     instance:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) kdelta([1,2],[2,1]);
          (%o2)                                 - 1
          (%i3) kdels([1,2],[2,1]);
          (%o3)                                  1
          (%i4) ishow(kdelta([a,b],[c,d]))$
                                       c       d         d       c
          (%t4)                  kdelta  kdelta  - kdelta  kdelta
                                       a       b         a       b
          (%i4) ishow(kdels([a,b],[c,d]))$
                                       c       d         d       c
          (%t4)                  kdelta  kdelta  + kdelta  kdelta
                                       a       b         a       b


 -- Function: levi_civita (<L>)
     is the permutation (or Levi-Civita) tensor which yields 1 if the
     list <L> consists of an even permutation of integers, -1 if it
     consists of an odd permutation, and 0 if some indices in <L> are
     repeated.


 -- Function: lc2kdt (<expr>)
     Simplifies expressions containing the Levi-Civita symbol,
     converting these to Kronecker-delta expressions when possible. The
     main difference between this function and simply evaluating the
     Levi-Civita symbol is that direct evaluation often results in
     Kronecker expressions containing numerical indices. This is often
     undesirable as it prevents further simplification.  The `lc2kdt'
     function avoids this problem, yielding expressions that are more
     easily simplified with `rename' or `contract'.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) expr:ishow('levi_civita([],[i,j])*'levi_civita([k,l],[])*a([j],[k]))$
                                            i j  k
          (%t2)                  levi_civita    a  levi_civita
                                                 j            k l
          (%i3) ishow(ev(expr,levi_civita))$
                                            i j  k       1 2
          (%t3)                       kdelta    a  kdelta
                                            1 2  j       k l
          (%i4) ishow(ev(%,kdelta))$
                       i       j         j       i   k
          (%t4) (kdelta  kdelta  - kdelta  kdelta ) a
                       1       2         1       2   j

                                         1       2         2       1
                                  (kdelta  kdelta  - kdelta  kdelta )
                                         k       l         k       l
          (%i5) ishow(lc2kdt(expr))$
                               k       i       j    k       j       i
          (%t5)               a  kdelta  kdelta  - a  kdelta  kdelta
                               j       k       l    j       k       l
          (%i6) ishow(contract(expand(%)))$
                                           i           i
          (%t6)                           a  - a kdelta
                                           l           l

     The `lc2kdt' function sometimes makes use of the metric tensor.
     If the metric tensor was not defined previously with `imetric',
     this results in an error.


          (%i7) expr:ishow('levi_civita([],[i,j])*'levi_civita([],[k,l])*a([j,k],[]))$
                                           i j            k l
          (%t7)                 levi_civita    levi_civita    a
                                                               j k
          (%i8) ishow(lc2kdt(expr))$
          Maxima encountered a Lisp error:

           Error in $IMETRIC [or a callee]:
           $IMETRIC [or a callee] requires less than two arguments.

          Automatically continuing.
          To reenable the Lisp debugger set *debugger-hook* to nil.
          (%i9) imetric(g);
          (%o9)                                done
          (%i10) ishow(lc2kdt(expr))$
                   %3 i       k   %4 j       l     %3 i       l   %4 j       k
          (%t10) (g     kdelta   g     kdelta   - g     kdelta   g     kdelta  ) a
                              %3             %4               %3             %4   j k
          (%i11) ishow(contract(expand(%)))$
                                            l i      l i
          (%t11)                           a    - a g


 -- Function: lc_l
     Simplification rule used for expressions containing the
     unevaluated Levi-Civita symbol (`levi_civita'). Along with `lc_u',
     it can be used to simplify many expressions more efficiently than
     the evaluation of `levi_civita'.  For example:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2)  el1:ishow('levi_civita([i,j,k],[])*a([],[i])*a([],[j]))$
                                       i  j
          (%t2)                       a  a  levi_civita
                                                       i j k
          (%i3) el2:ishow('levi_civita([],[i,j,k])*a([i])*a([j]))$
                                                 i j k
          (%t3)                       levi_civita      a  a
                                                        i  j
          (%i4) ishow(canform(contract(expand(applyb1(el1,lc_l,lc_u)))))$
          (%t4)                                  0
          (%i5) ishow(canform(contract(expand(applyb1(el2,lc_l,lc_u)))))$
          (%t5)                                  0


 -- Function: lc_u
     Simplification rule used for expressions containing the
     unevaluated Levi-Civita symbol (`levi_civita'). Along with `lc_u',
     it can be used to simplify many expressions more efficiently than
     the evaluation of `levi_civita'.  For details, see `lc_l'.


 -- Function: canten (<expr>)
     Simplifies <expr> by renaming (see `rename') and permuting dummy
     indices. `rename' is restricted to sums of tensor products in
     which no derivatives are present. As such it is limited and should
     only be used if `canform' is not capable of carrying out the
     required simplification.

     The `canten' function returns a mathematically correct result only
     if its argument is an expression that is fully symmetric in its
     indices.  For this reason, `canten' returns an error if `allsym'
     is not set to `true'.


 -- Function: concan (<expr>)
     Similar to `canten' but also performs index contraction.


28.2.2 Tensor symmetries
------------------------

 -- Option variable: allsym
     Default: `false'. if `true' then all indexed objects are assumed
     symmetric in all of their covariant and contravariant indices. If
     `false' then no symmetries of any kind are assumed in these
     indices. Derivative indices are always taken to be symmetric
     unless `iframe_flag' is set to `true'.


 -- Function: decsym (<tensor>, <m>, <n>, [<cov_1>, <cov_2>, ...],
          [<contr_1>, <contr_2>, ...])
     Declares symmetry properties for <tensor> of <m> covariant and <n>
     contravariant indices. The <cov_i> and <contr_i> are
     pseudofunctions expressing symmetry relations among the covariant
     and contravariant indices respectively.  These are of the form
     `symoper(<index_1>, <index_2>,...)' where `symoper' is one of
     `sym', `anti' or `cyc' and the <index_i> are integers indicating
     the position of the index in the <tensor>.  This will declare
     <tensor> to be symmetric, antisymmetric or cyclic respectively in
     the <index_i>. `symoper(all)' is also an allowable form which
     indicates all indices obey the symmetry condition. For example,
     given an object `b' with 5 covariant indices,
     `decsym(b,5,3,[sym(1,2),anti(3,4)],[cyc(all)])' declares `b'
     symmetric in its first and second and antisymmetric in its third
     and fourth covariant indices, and cyclic in all of its
     contravariant indices.  Either list of symmetry declarations may
     be null.  The function which performs the simplifications is
     `canform' as the example below illustrates.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) expr:contract(expand(a([i1,j1,k1],[])*kdels([i,j,k],[i1,j1,k1])))$
          (%i3) ishow(expr)$
          (%t3)         a      + a      + a      + a      + a      + a
                         k j i    k i j    j k i    j i k    i k j    i j k
          (%i4) decsym(a,3,0,[sym(all)],[]);
          (%o4)                                done
          (%i5) ishow(canform(expr))$
          (%t5)                              6 a
                                                i j k
          (%i6) remsym(a,3,0);
          (%o6)                                done
          (%i7) decsym(a,3,0,[anti(all)],[]);
          (%o7)                                done
          (%i8) ishow(canform(expr))$
          (%t8)                                  0
          (%i9) remsym(a,3,0);
          (%o9)                                done
          (%i10) decsym(a,3,0,[cyc(all)],[]);
          (%o10)                               done
          (%i11) ishow(canform(expr))$
          (%t11)                        3 a      + 3 a
                                           i k j      i j k
          (%i12) dispsym(a,3,0);
          (%o12)                     [[cyc, [[1, 2, 3]], []]]


 -- Function: remsym (<tensor>, <m>, <n>)
     Removes all symmetry properties from <tensor> which has <m>
     covariant indices and <n> contravariant indices.

 -- Function: canform (<expr>)
     Simplifies <expr> by renaming dummy indices and reordering all
     indices as dictated by symmetry conditions imposed on them. If
     `allsym' is `true' then all indices are assumed symmetric,
     otherwise symmetry information provided by `decsym' declarations
     will be used. The dummy indices are renamed in the same manner as
     in the `rename' function. When `canform' is applied to a large
     expression the calculation may take a considerable amount of time.
     This time can be shortened by calling `rename' on the expression
     first.  Also see the example under `decsym'. Note: `canform' may
     not be able to reduce an expression completely to its simplest
     form although it will always return a mathematically correct
     result.

28.2.3 Indicial tensor calculus
-------------------------------

 -- Function: diff (<expr>, <v_1>, [<n_1>, [<v_2>, <n_2>] ...])
     is the usual Maxima differentiation function which has been
     expanded in its abilities for `itensor'. It takes the derivative
     of <expr> with respect to <v_1> <n_1> times, with respect to <v_2>
     <n_2> times, etc. For the tensor package, the function has been
     modified so that the <v_i> may be integers from 1 up to the value
     of the variable `dim'.  This will cause the differentiation to be
     carried out with respect to the <v_i>th member of the list
     `vect_coords'.  If `vect_coords' is bound to an atomic variable,
     then that variable subscripted by <v_i> will be used for the
     variable of differentiation.  This permits an array of coordinate
     names or subscripted names like `x[1]', `x[2]', ...  to be used.

 -- Function: idiff (<expr>, <v_1>, [<n_1>, [<v_2>, <n_2>] ...])
     Indicial differentiation. Unlike `diff', which differentiates with
     respect to an independent variable, `idiff)' can be used to
     differentiate with respect to a coordinate. For an indexed object,
     this amounts to appending the <v_i> as derivative indices.
     Subsequently, derivative indices will be sorted, unless
     `iframe_flag' is set to `true'.

     `idiff' can also differentiate the determinant of the metric
     tensor. Thus, if `imetric' has been bound to `G' then
     `idiff(determinant(g),k)' will return
     `2*determinant(g)*ichr2([%i,k],[%i])' where the dummy index `%i'
     is chosen appropriately.


 -- Function: liediff (<v>, <ten>)
     Computes the Lie-derivative of the tensorial expression <ten> with
     respect to the vector field <v>. <ten> should be any indexed
     tensor expression; <v> should be the name (without indices) of a
     vector field. For example:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(liediff(v,a([i,j],[])*b([],[k],l)))$
                 k    %2            %2          %2
          (%t2) b   (v   a       + v   a     + v   a    )
                 ,l       i j,%2    ,j  i %2    ,i  %2 j

                                          %1  k        %1  k      %1  k
                                      + (v   b      - b   v    + v   b   ) a
                                              ,%1 l    ,l  ,%1    ,l  ,%1   i j


 -- Function: rediff (<ten>)
     Evaluates all occurrences of the `idiff' command in the tensorial
     expression <ten>.


 -- Function: undiff (<expr>)
     Returns an expression equivalent to <expr> but with all derivatives
     of indexed objects replaced by the noun form of the `idiff'
     function. Its arguments would yield that indexed object if the
     differentiation were carried out.  This is useful when it is
     desired to replace a differentiated indexed object with some
     function definition resulting in <expr> and then carry out the
     differentiation by saying `ev(<expr>, idiff)'.


 -- Function: evundiff (<expr>)
     Equivalent to the execution of `undiff', followed by `ev' and
     `rediff'.

     The point of this operation is to easily evalute expressions that
     cannot be directly evaluated in derivative form. For instance, the
     following causes an error:

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) icurvature([i,j,k],[l],m);
          Maxima encountered a Lisp error:

           Error in $ICURVATURE [or a callee]:
           $ICURVATURE [or a callee] requires less than three arguments.

          Automatically continuing.
          To reenable the Lisp debugger set *debugger-hook* to nil.

     However, if `icurvature' is entered in noun form, it can be
     evaluated using `evundiff':

          (%i3) ishow('icurvature([i,j,k],[l],m))$
                                                   l
          (%t3)                          icurvature
                                                   i j k,m
          (%i4) ishow(evundiff(%))$
                       l              l         %1           l           %1
          (%t4) - ichr2        - ichr2     ichr2      - ichr2       ichr2
                       i k,j m        %1 j      i k,m        %1 j,m      i k

                          l              l         %1           l           %1
                   + ichr2        + ichr2     ichr2      + ichr2       ichr2
                          i j,k m        %1 k      i j,m        %1 k,m      i j

     Note: In earlier versions of Maxima, derivative forms of the
     Christoffel-symbols also could not be evaluated. This has been
     fixed now, so `evundiff' is no longer necessary for expressions
     like this:

          (%i5) imetric(g);
          (%o5)                                done
          (%i6) ishow(ichr2([i,j],[k],l))$
                 k %3
                g     (g         - g         + g        )
                        j %3,i l    i j,%3 l    i %3,j l
          (%t6) -----------------------------------------
                                    2

                                   k %3
                                  g     (g       - g       + g      )
                                   ,l     j %3,i    i j,%3    i %3,j
                                + -----------------------------------
                                                   2


 -- Function: flush (<expr>, <tensor_1>, <tensor_2>, ...)
     Set to zero, in <expr>, all occurrences of the <tensor_i> that
     have no derivative indices.


 -- Function: flushd (<expr>, <tensor_1>, <tensor_2>, ...)
     Set to zero, in <expr>, all occurrences of the <tensor_i> that
     have derivative indices.


 -- Function: flushnd (<expr>, <tensor>, <n>)
     Set to zero, in <expr>, all occurrences of the differentiated
     object <tensor> that have <n> or more derivative indices as the
     following example demonstrates.

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(a([i],[J,r],k,r)+a([i],[j,r,s],k,r,s))$
                                          J r      j r s
          (%t2)                          a      + a
                                          i,k r    i,k r s
          (%i3) ishow(flushnd(%,a,3))$
                                               J r
          (%t3)                               a
                                               i,k r

 -- Function: coord (<tensor_1>, <tensor_2>, ...)
     Gives <tensor_i> the coordinate differentiation property that the
     derivative of contravariant vector whose name is one of the
     <tensor_i> yields a Kronecker delta. For example, if `coord(x)' has
     been done then `idiff(x([],[i]),j)' gives `kdelta([i],[j])'.
     `coord' is a list of all indexed objects having this property.


 -- Function: remcoord (<tensor_1>, <tensor_2>, ...)
 -- Function: remcoord (all)
     Removes the coordinate differentiation property from the `tensor_i'
     that was established by the function `coord'.  `remcoord(all)'
     removes this property from all indexed objects.


 -- Function: makebox (<expr>)
     Display <expr> in the same manner as `show'; however, any tensor
     d'Alembertian occurring in <expr> will be indicated using the
     symbol `[]'.  For example, `[]p([m],[n])' represents
     `g([],[i,j])*p([m],[n],i,j)'.


 -- Function: conmetderiv (<expr>, <tensor>)
     Simplifies expressions containing ordinary derivatives of both
     covariant and contravariant forms of the metric tensor (the
     current restriction).  For example, `conmetderiv' can relate the
     derivative of the contravariant metric tensor with the Christoffel
     symbols as seen from the following:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(g([],[a,b],c))$
                                                a b
          (%t2)                                g
                                                ,c
          (%i3) ishow(conmetderiv(%,g))$
                                   %1 b      a       %1 a      b
          (%t3)                 - g     ichr2     - g     ichr2
                                             %1 c              %1 c

 -- Function: simpmetderiv (<expr>)
 -- Function: simpmetderiv (<expr>[, <stop>])
     Simplifies expressions containing products of the derivatives of
     the metric tensor. Specifically, `simpmetderiv' recognizes two
     identities:


             ab        ab           ab                 a
            g   g   + g   g     = (g   g  )   = (kdelta )   = 0
             ,d  bc        bc,d         bc ,d          c ,d

     hence


             ab          ab
            g   g   = - g   g
             ,d  bc          bc,d

     and


            ab          ab
           g   g     = g   g
            ,j  ab,i    ,i  ab,j

     which follows from the symmetries of the Christoffel symbols.

     The `simpmetderiv' function takes one optional parameter which,
     when present, causes the function to stop after the first
     successful substitution in a product expression. The
     `simpmetderiv' function also makes use of the global variable
     <flipflag> which determines how to apply a "canonical" ordering to
     the product indices.

     Put together, these capabilities can be used to achieve powerful
     simplifications that are difficult or impossible to accomplish
     otherwise.  This is demonstrated through the following example
     that explicitly uses the partial simplification features of
     `simpmetderiv' to obtain a contractible expression:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) imetric(g);
          (%o2)                                done
          (%i3) ishow(g([],[a,b])*g([],[b,c])*g([a,b],[],d)*g([b,c],[],e))$
                                       a b  b c
          (%t3)                       g    g    g      g
                                                 a b,d  b c,e
          (%i4) ishow(canform(%))$

          errexp1 has improper indices
           -- an error.  Quitting.  To debug this try debugmode(true);
          (%i5) ishow(simpmetderiv(%))$
                                       a b  b c
          (%t5)                       g    g    g      g
                                                 a b,d  b c,e
          (%i6) flipflag:not flipflag;
          (%o6)                                true
          (%i7) ishow(simpmetderiv(%th(2)))$
                                         a b  b c
          (%t7)                         g    g    g    g
                                         ,d   ,e   a b  b c
          (%i8) flipflag:not flipflag;
          (%o8)                                false
          (%i9) ishow(simpmetderiv(%th(2),stop))$
                                         a b  b c
          (%t9)                       - g    g    g      g
                                              ,e   a b,d  b c
          (%i10) ishow(contract(%))$
                                              b c
          (%t10)                           - g    g
                                              ,e   c b,d

     See also `weyl.dem' for an example that uses `simpmetderiv' and
     `conmetderiv' together to simplify contractions of the Weyl tensor.


 -- Function: flush1deriv (<expr>, <tensor>)
     Set to zero, in `expr', all occurrences of `tensor' that have
     exactly one derivative index.


28.2.4 Tensors in curved spaces
-------------------------------

 -- Function: imetric (<g>)
 -- System variable: imetric
     Specifies the metric by assigning the variable `imetric:<g>' in
     addition, the contraction properties of the metric <g> are set up
     by executing the commands `defcon(<g>),defcon(<g>,<g>,kdelta)'.
     The variable `imetric' (unbound by default), is bound to the
     metric, assigned by the `imetric(<g>)' command.


 -- Function: idim (<n>)
     Sets the dimensions of the metric. Also initializes the
     antisymmetry properties of the Levi-Civita symbols for the given
     dimension.


 -- Function: ichr1 ([<i>, <j>, <k>])
     Yields the Christoffel symbol of the first kind via the definition
                 (g      + g      - g     )/2 .
                   ik,j     jk,i     ij,k
     To evaluate the Christoffel symbols for a particular metric, the
     variable `imetric' must be assigned a name as in the example under
     `chr2'.


 -- Function: ichr2 ([<i>, <j>], [<k>])
     Yields the Christoffel symbol of the second kind defined by the
     relation
                                 ks
             ichr2([i,j],[k]) = g    (g      + g      - g     )/2
                                       is,j     js,i     ij,s

 -- Function: icurvature ([<i>, <j>, <k>], [<h>])
     Yields the Riemann curvature tensor in terms of the Christoffel
     symbols of the second kind (`ichr2').  The following notation is
     used:
                         h             h            h         %1         h
               icurvature     = - ichr2      - ichr2     ichr2    + ichr2
                         i j k         i k,j        %1 j      i k        i j,k
                                         h          %1
                                  + ichr2      ichr2
                                         %1 k       i j

 -- Function: covdiff (<expr>, <v_1>, <v_2>, ...)
     Yields the covariant derivative of <expr> with respect to the
     variables <v_i> in terms of the Christoffel symbols of the second
     kind (`ichr2').  In order to evaluate these, one should use
     `ev(<expr>,ichr2)'.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) entertensor()$
          Enter tensor name: a;
          Enter a list of the covariant indices: [i,j];
          Enter a list of the contravariant indices: [k];
          Enter a list of the derivative indices: [];
                                                k
          (%t2)                                a
                                                i j
          (%i3) ishow(covdiff(%,s))$
                       k         %1     k         %1     k            k     %1
          (%t3)     - a     ichr2    - a     ichr2    + a      + ichr2     a
                       i %1      j s    %1 j      i s    i j,s        %1 s  i j
          (%i4) imetric:g;
          (%o4)                                  g
          (%i5) ishow(ev(%th(2),ichr2))$
                   %1 %4  k
                  g      a     (g       - g       + g      )
                          i %1   s %4,j    j s,%4    j %4,s
          (%t5) - ------------------------------------------
                                      2
              %1 %3  k
             g      a     (g       - g       + g      )
                     %1 j   s %3,i    i s,%3    i %3,s
           - ------------------------------------------
                                 2
              k %2  %1
             g     a    (g        - g        + g       )
                    i j   s %2,%1    %1 s,%2    %1 %2,s     k
           + ------------------------------------------- + a
                                  2                         i j,s
          (%i6)


 -- Function: lorentz_gauge (<expr>)
     Imposes the Lorentz condition by substituting 0 for all indexed
     objects in <expr> that have a derivative index identical to a
     contravariant index.


 -- Function: igeodesic_coords (<expr>, <name>)
     Causes undifferentiated Christoffel symbols and first derivatives
     of the metric tensor vanish in <expr>. The <name> in the
     `igeodesic_coords' function refers to the metric <name> (if it
     appears in <expr>) while the connection coefficients must be
     called with the names `ichr1' and/or `ichr2'. The following example
     demonstrates the verification of the cyclic identity satisfied by
     the Riemann curvature tensor using the `igeodesic_coords' function.


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(icurvature([r,s,t],[u]))$
                       u            u         %1         u            u         %1
          (%t2) - ichr2      - ichr2     ichr2    + ichr2      + ichr2     ichr2
                       r t,s        %1 s      r t        r s,t        %1 t      r s
          (%i3) ishow(igeodesic_coords(%,ichr2))$
                                           u            u
          (%t3)                       ichr2      - ichr2
                                           r s,t        r t,s
          (%i4) ishow(igeodesic_coords(icurvature([r,s,t],[u]),ichr2)+
                      igeodesic_coords(icurvature([s,t,r],[u]),ichr2)+
                      igeodesic_coords(icurvature([t,r,s],[u]),ichr2))$
                       u            u            u            u            u
          (%t4) - ichr2      + ichr2      + ichr2      - ichr2      - ichr2
                       t s,r        t r,s        s t,r        s r,t        r t,s

                                                                            u
                                                                     + ichr2
                                                                            r s,t
          (%i5) canform(%);
          (%o5)                                  0


28.2.5 Moving frames
--------------------

Maxima now has the ability to perform calculations using moving frames.
These can be orthonormal frames (tetrads, vielbeins) or an arbitrary
frame.

   To use frames, you must first set `iframe_flag' to `true'. This
causes the Christoffel-symbols, `ichr1' and `ichr2', to be replaced by
the more general frame connection coefficients `icc1' and `icc2' in
calculations. Speficially, the behavior of `covdiff' and `icurvature'
is changed.

   The frame is defined by two tensors: the inverse frame field (`ifri',
the dual basis tetrad), and the frame metric `ifg'. The frame metric is
the identity matrix for orthonormal frames, or the Lorentz metric for
orthonormal frames in Minkowski spacetime. The inverse frame field
defines the frame base (unit vectors).  Contraction properties are
defined for the frame field and the frame metric.

   When `iframe_flag' is true, many `itensor' expressions use the frame
metric `ifg' instead of the metric defined by `imetric' for raising and
lowerind indices.

   IMPORTANT: Setting the variable `iframe_flag' to `true' does NOT
undefine the contraction properties of a metric defined by a call to
`defcon' or `imetric'. If a frame field is used, it is best to define
the metric by assigning its name to the variable `imetric' and NOT
invoke the `imetric' function.

   Maxima uses these two tensors to define the frame coefficients
(`ifc1' and `ifc2') which form part of the connection coefficients
(`icc1' and `icc2'), as the following example demonstrates:


     (%i1) load(itensor);
     (%o1)      /share/tensor/itensor.lisp
     (%i2) iframe_flag:true;
     (%o2)                                true
     (%i3) ishow(covdiff(v([],[i]),j))$
                                    i        i     %1
     (%t3)                         v   + icc2     v
                                    ,j       %1 j
     (%i4) ishow(ev(%,icc2))$
                             %1      i           i        i
     (%t4)                  v   (ifc2     + ichr2    ) + v
                                     %1 j        %1 j     ,j
     (%i5) ishow(ev(%,ifc2))$
                 %1    i %2
                v   ifg     (ifb        - ifb        + ifb       )
                                j %2 %1      %2 %1 j      %1 j %2     i
     (%t5)      -------------------------------------------------- + v
                                        2                             ,j
     (%i6) ishow(ifb([a,b,c]))$
                            %5    %4
     (%t6)               ifr   ifr   (ifri        - ifri       )
                            a     b       c %4,%5       c %5,%4

   An alternate method is used to compute the frame bracket (`ifb') if
the `iframe_bracket_form' flag is set to `false':


     (%i8) block([iframe_bracket_form:false],ishow(ifb([a,b,c])))$
                            %7    %6        %6      %7
     (%t8)              (ifr   ifr     - ifr     ifr  ) ifri
                            a     b,%7      a,%7    b       c %6

 -- Function: iframes ()
     Since in this version of Maxima, contraction identities for `ifr'
     and `ifri' are always defined, as is the frame bracket (`ifb'),
     this function does nothing.


 -- Variable: ifb
     The frame bracket. The contribution of the frame metric to the
     connection coefficients is expressed using the frame bracket:


                    - ifb      + ifb      + ifb
                         c a b      b c a      a b c
          ifc1    = --------------------------------
              abc                  2

     The frame bracket itself is defined in terms of the frame field
     and frame metric. Two alternate methods of computation are used
     depending on the value of `frame_bracket_form'. If true (the
     default) or if the `itorsion_flag' is `true':


                    d      e                                      f
          ifb =  ifr    ifr   (ifri      - ifri      - ifri    itr   )
             abc    b      c       a d,e       a e,d       a f    d e

     Otherwise:


                       e      d        d      e
          ifb    = (ifr    ifr    - ifr    ifr   ) ifri
             abc       b      c,e      b,e    c        a d


 -- Variable: icc1
     Connection coefficients of the first kind. In `itensor', defined as


          icc1    = ichr1    - ikt1    - inmc1
              abc        abc       abc        abc

     In this expression, if `iframe_flag' is true, the
     Christoffel-symbol `ichr1' is replaced with the frame connection
     coefficient `ifc1'.  If `itorsion_flag' is `false', `ikt1' will be
     omitted. It is also omitted if a frame base is used, as the
     torsion is already calculated as part of the frame bracket.
     Lastly, of `inonmet_flag' is `false', `inmc1' will not be present.


 -- Variable: icc2
     Connection coefficients of the second kind. In `itensor', defined
     as


              c         c        c         c
          icc2   = ichr2   - ikt2   - inmc2
              ab        ab       ab        ab

     In this expression, if `iframe_flag' is true, the
     Christoffel-symbol `ichr2' is replaced with the frame connection
     coefficient `ifc2'.  If `itorsion_flag' is `false', `ikt2' will be
     omitted. It is also omitted if a frame base is used, as the
     torsion is already calculated as part of the frame bracket.
     Lastly, of `inonmet_flag' is `false', `inmc2' will not be present.


 -- Variable: ifc1
     Frame coefficient of the first kind (also known as Ricci-rotation
     coefficients.) This tensor represents the contribution of the
     frame metric to the connection coefficient of the first kind.
     Defined as:


                    - ifb      + ifb      + ifb
                         c a b      b c a      a b c
          ifc1    = --------------------------------
              abc                   2


 -- Variable: ifc2
     Frame coefficient of the first kind. This tensor represents the
     contribution of the frame metric to the connection coefficient of
     the first kind. Defined as a permutation of the frame bracket
     (`ifb') with the appropriate indices raised and lowered as
     necessary:


              c       cd
          ifc2   = ifg   ifc1
              ab             abd


 -- Variable: ifr
     The frame field. Contracts with the inverse frame field (`ifri') to
     form the frame metric (`ifg').


 -- Variable: ifri
     The inverse frame field. Specifies the frame base (dual basis
     vectors). Along with the frame metric, it forms the basis of all
     calculations based on frames.


 -- Variable: ifg
     The frame metric. Defaults to `kdelta', but can be changed using
     `components'.


 -- Variable: ifgi
     The inverse frame metric. Contracts with the frame metric (`ifg')
     to `kdelta'.


 -- Option variable: iframe_bracket_form
     Default value: `true'

     Specifies how the frame bracket (`ifb') is computed.


28.2.6 Torsion and nonmetricity
-------------------------------

Maxima can now take into account torsion and nonmetricity. When the flag
`itorsion_flag' is set to `true', the contribution of torsion is added
to the connection coefficients. Similarly, when the flag `inonmet_flag'
is true, nonmetricity components are included.

 -- Variable: inm
     The nonmetricity vector. Conformal nonmetricity is defined through
     the covariant derivative of the metric tensor. Normally zero, the
     metric tensor's covariant derivative will evaluate to the
     following when `inonmet_flag' is set to `true':


          g     =- g  inm
           ij;k     ij   k


 -- Variable: inmc1
     Covariant permutation of the nonmetricity vector components.
     Defined as


                     g   inm  - inm  g   - g   inm
                      ab    c      a  bc    ac    b
          inmc1    = ------------------------------
               abc                 2

     (Substitute `ifg' in place of `g' if a frame metric is used.)


 -- Variable: inmc2
     Contravariant permutation of the nonmetricity vector components.
     Used in the connection coefficients if `inonmet_flag' is `true'.
     Defined as:


                                c         c         cd
                    -inm  kdelta  - kdelta  inm  + g   inm  g
               c        a       b         a    b          d  ab
          inmc2   = -------------------------------------------
               ab                        2

     (Substitute `ifg' in place of `g' if a frame metric is used.)


 -- Variable: ikt1
     Covariant permutation of the torsion tensor (also known as
     contorsion).  Defined as:


                            d           d       d
                    -g   itr  - g    itr   - itr   g
                      ad    cb    bd    ca      ab  cd
          ikt1    = ----------------------------------
              abc                   2

     (Substitute `ifg' in place of `g' if a frame metric is used.)


 -- Variable: ikt2
     Contravariant permutation of the torsion tensor (also known as
     contorsion).  Defined as:


              c     cd
          ikt2   = g   ikt1
              ab           abd

     (Substitute `ifg' in place of `g' if a frame metric is used.)


 -- Variable: itr
     The torsion tensor. For a metric with torsion, repeated covariant
     differentiation on a scalar function will not commute, as
     demonstrated by the following example:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) imetric:g;
          (%o2)                                  g
          (%i3) covdiff(covdiff(f([],[]),i),j)-covdiff(covdiff(f([],[]),j),i)$
          (%i4) ishow(%)$
                                             %4              %2
          (%t4)                    f    ichr2    - f    ichr2
                                    ,%4      j i    ,%2      i j
          (%i5) canform(%);
          (%o5)                                  0
          (%i6) itorsion_flag:true;
          (%o6)                                true
          (%i7) covdiff(covdiff(f([],[]),i),j)-covdiff(covdiff(f([],[]),j),i)$
          (%i8) ishow(%)$
                                     %8             %6
          (%t8)             f    icc2    - f    icc2    - f     + f
                             ,%8     j i    ,%6     i j    ,j i    ,i j
          (%i9) ishow(canform(%))$
                                             %1             %1
          (%t9)                     f    icc2    - f    icc2
                                     ,%1     j i    ,%1     i j
          (%i10) ishow(canform(ev(%,icc2)))$
                                             %1             %1
          (%t10)                    f    ikt2    - f    ikt2
                                     ,%1     i j    ,%1     j i
          (%i11) ishow(canform(ev(%,ikt2)))$
                                %2 %1                    %2 %1
          (%t11)          f    g      ikt1       - f    g      ikt1
                           ,%2            i j %1    ,%2            j i %1
          (%i12) ishow(factor(canform(rename(expand(ev(%,ikt1))))))$
                                     %3 %2            %1       %1
                               f    g      g      (itr    - itr   )
                                ,%3         %2 %1     j i      i j
          (%t12)               ------------------------------------
                                                2
          (%i13) decsym(itr,2,1,[anti(all)],[]);
          (%o13)                               done
          (%i14) defcon(g,g,kdelta);
          (%o14)                               done
          (%i15) subst(g,nounify(g),%th(3))$
          (%i16) ishow(canform(contract(%)))$
                                                     %1
          (%t16)                           - f    itr
                                              ,%1    i j


28.2.7 Exterior algebra
-----------------------

The `itensor' package can perform operations on totally antisymmetric
covariant tensor fields. A totally antisymmetric tensor field of rank
(0,L) corresponds with a differential L-form. On these objects, a
multiplication operation known as the exterior product, or wedge
product, is defined.

   Unfortunately, not all authors agree on the definition of the wedge
product. Some authors prefer a definition that corresponds with the
notion of antisymmetrization: in these works, the wedge product of two
vector fields, for instance, would be defined as

                 a a  - a a
                  i j    j i
      a  /\ a  = -----------
       i     j        2

   More generally, the product of a p-form and a q-form would be
defined as

                            1     k1..kp l1..lq
     A       /\ B       = ------ D              A       B
      i1..ip     j1..jq   (p+q)!  i1..ip j1..jq  k1..kp  l1..lq

   where `D' stands for the Kronecker-delta.

   Other authors, however, prefer a "geometric" definition that
corresponds with the notion of the volume element:

     a  /\ a  = a a  - a a
      i     j    i j    j i

   and, in the general case

                            1    k1..kp l1..lq
     A       /\ B       = ----- D              A       B
      i1..ip     j1..jq   p! q!  i1..ip j1..jq  k1..kp  l1..lq

   Since `itensor' is a tensor algebra package, the first of these two
definitions appears to be the more natural one. Many applications,
however, utilize the second definition. To resolve this dilemma, a flag
has been implemented that controls the behavior of the wedge product: if
`igeowedge_flag' is `false' (the default), the first, "tensorial"
definition is used, otherwise the second, "geometric" definition will
be applied.

 -- Operator: ~
     The wedge product operator is denoted by the tilde `~'. This is a
     binary operator. Its arguments should be expressions involving
     scalars, covariant tensors of rank one, or covariant tensors of
     rank `l' that have been declared antisymmetric in all covariant
     indices.

     The behavior of the wedge product operator is controlled by the
     `igeowedge_flag' flag, as in the following example:

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(a([i])~b([j]))$
                                           a  b  - b  a
                                            i  j    i  j
          (%t2)                            -------------
                                                 2
          (%i3) decsym(a,2,0,[anti(all)],[]);
          (%o3)                                done
          (%i4) ishow(a([i,j])~b([k]))$
                                    a    b  + b  a    - a    b
                                     i j  k    i  j k    i k  j
          (%t4)                     ---------------------------
                                                 3
          (%i5) igeowedge_flag:true;
          (%o5)                                true
          (%i6) ishow(a([i])~b([j]))$
          (%t6)                            a  b  - b  a
                                            i  j    i  j
          (%i7) ishow(a([i,j])~b([k]))$
          (%t7)                     a    b  + b  a    - a    b
                                     i j  k    i  j k    i k  j


 -- Operator: |
     The vertical bar `|' denotes the "contraction with a vector" binary
     operation. When a totally antisymmetric covariant tensor is
     contracted with a contravariant vector, the result is the same
     regardless which index was used for the contraction. Thus, it is
     possible to define the contraction operation in an index-free
     manner.

     In the `itensor' package, contraction with a vector is always
     carried out with respect to the first index in the literal sorting
     order. This ensures better simplification of expressions involving
     the `|' operator. For instance:

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) decsym(a,2,0,[anti(all)],[]);
          (%o2)                                done
          (%i3) ishow(a([i,j],[])|v)$
                                              %1
          (%t3)                              v   a
                                                  %1 j
          (%i4) ishow(a([j,i],[])|v)$
                                               %1
          (%t4)                             - v   a
                                                   %1 j

     Note that it is essential that the tensors used with the `|'
     operator be declared totally antisymmetric in their covariant
     indices. Otherwise, the results will be incorrect.


 -- Function: extdiff (<expr>, <i>)
     Computes the exterior derivative of <expr> with respect to the
     index <i>. The exterior derivative is formally defined as the wedge
     product of the partial derivative operator and a differential
     form. As such, this operation is also controlled by the setting of
     `igeowedge_flag'.  For instance:

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) ishow(extdiff(v([i]),j))$
                                            v    - v
                                             j,i    i,j
          (%t2)                             -----------
                                                 2
          (%i3) decsym(a,2,0,[anti(all)],[]);
          (%o3)                                done
          (%i4) ishow(extdiff(a([i,j]),k))$
                                     a      - a      + a
                                      j k,i    i k,j    i j,k
          (%t4)                      ------------------------
                                                3
          (%i5) igeowedge_flag:true;
          (%o5)                                true
          (%i6) ishow(extdiff(v([i]),j))$
          (%t6)                             v    - v
                                             j,i    i,j
          (%i7) ishow(extdiff(a([i,j]),k))$
          (%t7)                      a      - a      + a
                                      j k,i    i k,j    i j,k


 -- Function: hodge (<expr>)
     Compute the Hodge-dual of <expr>. For instance:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) imetric(g);
          (%o2)                            done
          (%i3) idim(4);
          (%o3)                            done
          (%i4) icounter:100;
          (%o4)                             100
          (%i5) decsym(A,3,0,[anti(all)],[])$

          (%i6) ishow(A([i,j,k],[]))$
          (%t6)                           A
                                           i j k
          (%i7) ishow(canform(hodge(%)))$
                                    %1 %2 %3 %4
                         levi_civita            g        A
                                                 %1 %102  %2 %3 %4
          (%t7)          -----------------------------------------
                                             6
          (%i8) ishow(canform(hodge(%)))$
                           %1 %2 %3 %8            %4 %5 %6 %7
          (%t8) levi_civita            levi_civita            g        g
                                                               %1 %106  %2 %107
                                                      g        g      A        /6
                                                       %3 %108  %4 %8  %5 %6 %7
          (%i9) lc2kdt(%)$

          (%i10) %,kdelta$

          (%i11) ishow(canform(contract(expand(%))))$
          (%t11)                     - A
                                        %106 %107 %108


 -- Option variable: igeowedge_flag
     Default value: `false'

     Controls the behavior of the wedge product and exterior
     derivative. When set to `false' (the default), the notion of
     differential forms will correspond with that of a totally
     antisymmetric covariant tensor field.  When set to `true',
     differential forms will agree with the notion of the volume
     element.


28.2.8 Exporting TeX expressions
--------------------------------

The `itensor' package provides limited support for exporting tensor
expressions to TeX. Since `itensor' expressions appear as function
calls, the regular Maxima `tex' command will not produce the expected
output. You can try instead the `tentex' command, which attempts to
translate tensor expressions into appropriately indexed TeX objects.

 -- Function: tentex (<expr>)
     To use the `tentex' function, you must first load `tentex', as in
     the following example:


          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) load(tentex);
          (%o2)       /share/tensor/tentex.lisp
          (%i3) idummyx:m;
          (%o3)                                  m
          (%i4) ishow(icurvature([j,k,l],[i]))$
                      m1       i           m1       i           i            i
          (%t4)  ichr2    ichr2     - ichr2    ichr2     - ichr2      + ichr2
                      j k      m1 l        j l      m1 k        j l,k        j k,l
          (%i5) tentex(%)$
          $$\Gamma_{j\,k}^{m_1}\,\Gamma_{l\,m_1}^{i}-\Gamma_{j\,l}^{m_1}\,
           \Gamma_{k\,m_1}^{i}-\Gamma_{j\,l,k}^{i}+\Gamma_{j\,k,l}^{i}$$

     Note the use of the `idummyx' assignment, to avoid the appearance
     of the percent sign in the TeX expression, which may lead to
     compile errors.

     NB: This version of the `tentex' function is somewhat experimental.


28.2.9 Interfacing with ctensor
-------------------------------

The `itensor' package has the ability to generate Maxima code that can
then be executed in the context of the `ctensor' package. The function
that performs this task is `ic_convert'.

 -- Function: ic_convert (<eqn>)
     Converts the `itensor' equation <eqn> to a `ctensor' assignment
     statement.  Implied sums over dummy indices are made explicit
     while indexed objects are transformed into arrays (the array
     subscripts are in the order of covariant followed by contravariant
     indices of the indexed objects). The derivative of an indexed
     object will be replaced by the noun form of `diff' taken with
     respect to `ct_coords' subscripted by the derivative index. The
     Christoffel symbols `ichr1' and `ichr2' will be translated to
     `lcs' and `mcs', respectively and if `metricconvert' is `true'
     then all occurrences of the metric with two covariant
     (contravariant) indices will be renamed to `lg' (`ug'). In
     addition, `do' loops will be introduced summing over all free
     indices so that the transformed assignment statement can be
     evaluated by just doing `ev'. The following examples demonstrate
     the features of this function.

          (%i1) load(itensor);
          (%o1)      /share/tensor/itensor.lisp
          (%i2) eqn:ishow(t([i,j],[k])=f([],[])*g([l,m],[])*a([],[m],j)*b([i],[l,k]))$
                                       k        m   l k
          (%t2)                       t    = f a   b    g
                                       i j      ,j  i    l m
          (%i3) ic_convert(eqn);
          (%o3) for i thru dim do (for j thru dim

          do (for k thru dim do t        : f sum(sum(diff(a , ct_coords ) b
                                 i, j, k                   m           j   i, l, k

           g    , l, 1, dim), m, 1, dim)))
            l, m
          (%i4) imetric(g);
          (%o4)                                done
          (%i5) metricconvert:true;
          (%o5)                                true
          (%i6) ic_convert(eqn);
          (%o6) for i thru dim do (for j thru dim

          do (for k thru dim do t        : f sum(sum(diff(a , ct_coords ) b
                                 i, j, k                   m           j   i, l, k

           lg    , l, 1, dim), m, 1, dim)))
             l, m


28.2.10 Reserved words
----------------------

The following Maxima words are used by the `itensor' package internally
and should not be redefined:

       Keyword    Comments
       ------------------------------------------
       indices2() Internal version of indices()
       conti      Lists contravariant indices
       covi       Lists covariant indices of a indexed object
       deri       Lists derivative indices of an indexed object
       name       Returns the name of an indexed object
       concan
       irpmon
       lc0
       _lc2kdt0
       _lcprod
       _extlc


File: maxima.info,  Node: ctensor,  Next: atensor,  Prev: itensor,  Up: Top

29 ctensor
**********

* Menu:

* Introduction to ctensor::
* Definitions for ctensor::


File: maxima.info,  Node: Introduction to ctensor,  Next: Definitions for ctensor,  Prev: ctensor,  Up: ctensor

29.1 Introduction to ctensor
============================

`ctensor' is a component tensor manipulation package.  To use the
`ctensor' package, type `load(ctensor)'.  To begin an interactive
session with `ctensor', type `csetup()'.  You are first asked to
specify the dimension of the manifold. If the dimension is 2, 3 or 4
then the list of coordinates defaults to `[x,y]', `[x,y,z]' or
`[x,y,z,t]' respectively.  These names may be changed by assigning a
new list of coordinates to the variable `ct_coords' (described below)
and the user is queried about this. Care must be taken to avoid the
coordinate names conflicting with other object definitions.

   Next, the user enters the metric either directly or from a file by
specifying its ordinal position.  The metric is stored in the matrix
`lg'. Finally, the metric inverse is computed and stored in the matrix
`ug'. One has the option of carrying out all calculations in a power
series.

   A sample protocol is begun below for the static, spherically
symmetric metric (standard coordinates) which will be applied to the
problem of deriving Einstein's vacuum equations (which lead to the
Schwarzschild solution) as an example. Many of the functions in
`ctensor' will be displayed for the standard metric as examples.

     (%i1) load(ctensor);
     (%o1)      /share/tensor/ctensor.mac
     (%i2) csetup();
     Enter the dimension of the coordinate system:
     4;
     Do you wish to change the coordinate names?
     n;
     Do you want to
     1. Enter a new metric?

     2. Enter a metric from a file?

     3. Approximate a metric with a Taylor series?
     1;

     Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
     Answer 1, 2, 3 or 4
     1;
     Row 1 Column 1:
     a;
     Row 2 Column 2:
     x^2;
     Row 3 Column 3:
     x^2*sin(y)^2;
     Row 4 Column 4:
     -d;

     Matrix entered.
     Enter functional dependencies with the DEPENDS function or 'N' if none
     depends([a,d],x);
     Do you wish to see the metric?
     y;
                               [ a  0       0        0  ]
                               [                        ]
                               [     2                  ]
                               [ 0  x       0        0  ]
                               [                        ]
                               [         2    2         ]
                               [ 0  0   x  sin (y)   0  ]
                               [                        ]
                               [ 0  0       0       - d ]
     (%o2)                                done
     (%i3) christof(mcs);
                                                 a
                                                  x
     (%t3)                          mcs        = ---
                                       1, 1, 1   2 a

                                                  1
     (%t4)                           mcs        = -
                                        1, 2, 2   x

                                                  1
     (%t5)                           mcs        = -
                                        1, 3, 3   x

                                                 d
                                                  x
     (%t6)                          mcs        = ---
                                       1, 4, 4   2 d

                                                   x
     (%t7)                          mcs        = - -
                                       2, 2, 1     a

                                                cos(y)
     (%t8)                         mcs        = ------
                                      2, 3, 3   sin(y)

                                                    2
                                               x sin (y)
     (%t9)                      mcs        = - ---------
                                   3, 3, 1         a

     (%t10)                   mcs        = - cos(y) sin(y)
                                 3, 3, 2

                                                 d
                                                  x
     (%t11)                         mcs        = ---
                                       4, 4, 1   2 a
     (%o11)                               done


File: maxima.info,  Node: Definitions for ctensor,  Prev: Introduction to ctensor,  Up: ctensor

29.2 Definitions for ctensor
============================

29.2.1 Initialization and setup
-------------------------------

 -- Function: csetup ()
     A function in the `ctensor' (component tensor) package which
     initializes the package and allows the user to enter a metric
     interactively. See `ctensor' for more details.

 -- Function: cmetric (<dis>)
 -- Function: cmetric ()
     A function in the `ctensor' (component tensor) package that
     computes the metric inverse and sets up the package for further
     calculations.

     If `cframe_flag' is false, the function computes the inverse metric
     `ug' from the (user-defined) matrix `lg'. The metric determinant is
     also computed and stored in the variable `gdet'. Furthermore, the
     package determines if the metric is diagonal and sets the value of
     `diagmetric' accordingly. If the optional argument <dis> is
     present and not equal to `false', the user is prompted to see the
     metric inverse.

     If `cframe_flag' is `true', the function expects that the values of
     `fri' (the inverse frame matrix) and `lfg' (the frame metric) are
     defined. From these, the frame matrix `fr' and the inverse frame
     metric `ufg' are computed.


 -- Function: ct_coordsys (<coordinate_system>, <extra_arg>)
 -- Function: ct_coordsys (<coordinate_system>)
     Sets up a predefined coordinate system and metric. The argument
     <coordinate_system> can be one of the following symbols:


            SYMBOL               Dim Coordinates       Description/comments
            --------------------------------------------------------------------------
            cartesian2d           2  [x,y]             Cartesian 2D coordinate system
            polar                 2  [r,phi]           Polar coordinate system
            elliptic              2  [u,v]             Elliptic coordinate system
            confocalelliptic      2  [u,v]             Confocal elliptic coordinates
            bipolar               2  [u,v]             Bipolar coordinate system
            parabolic             2  [u,v]             Parabolic coordinate system
            cartesian3d           3  [x,y,z]           Cartesian 3D coordinate system
            polarcylindrical      3  [r,theta,z]       Polar 2D with cylindrical z
            ellipticcylindrical   3  [u,v,z]           Elliptic 2D with cylindrical z
            confocalellipsoidal   3  [u,v,w]           Confocal ellipsoidal
            bipolarcylindrical    3  [u,v,z]           Bipolar 2D with cylindrical z
            paraboliccylindrical  3  [u,v,z]           Parabolic 2D with cylindrical z
            paraboloidal          3  [u,v,phi]         Paraboloidal coordinates
            conical               3  [u,v,w]           Conical coordinates
            toroidal              3  [u,v,phi]         Toroidal coordinates
            spherical             3  [r,theta,phi]     Spherical coordinate system
            oblatespheroidal      3  [u,v,phi]         Oblate spheroidal coordinates
            oblatespheroidalsqrt  3  [u,v,phi]
            prolatespheroidal     3  [u,v,phi]         Prolate spheroidal coordinates
            prolatespheroidalsqrt 3  [u,v,phi]
            ellipsoidal           3  [r,theta,phi]     Ellipsoidal coordinates
            cartesian4d           4  [x,y,z,t]         Cartesian 4D coordinate system
            spherical4d           4  [r,theta,eta,phi] Spherical 4D coordinate system
            exteriorschwarzschild 4  [t,r,theta,phi]   Schwarzschild metric
            interiorschwarzschild 4  [t,z,u,v]         Interior Schwarzschild metric
            kerr_newman           4  [t,r,theta,phi]   Charged axially symmetric metric

     `coordinate_system' can also be a list of transformation functions,
     followed by a list containing the coordinate variables. For
     instance, you can specify a spherical metric as follows:


          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) ct_coordsys([r*cos(theta)*cos(phi),r*cos(theta)*sin(phi),
                r*sin(theta),[r,theta,phi]]);
          (%o2)                                done
          (%i3) lg:trigsimp(lg);
                                     [ 1  0         0        ]
                                     [                       ]
                                     [     2                 ]
          (%o3)                      [ 0  r         0        ]
                                     [                       ]
                                     [         2    2        ]
                                     [ 0  0   r  cos (theta) ]
          (%i4) ct_coords;
          (%o4)                           [r, theta, phi]
          (%i5) dim;
          (%o5)                                  3

     Transformation functions can also be used when `cframe_flag' is
     `true':


          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) cframe_flag:true;
          (%o2)                                true
          (%i3) ct_coordsys([r*cos(theta)*cos(phi),r*cos(theta)*sin(phi),
                r*sin(theta),[r,theta,phi]]);
          (%o3)                                done
          (%i4) fri;
                [ cos(phi) cos(theta)  - cos(phi) r sin(theta)  - sin(phi) r cos(theta) ]
                [                                                                       ]
          (%o4) [ sin(phi) cos(theta)  - sin(phi) r sin(theta)   cos(phi) r cos(theta)  ]
                [                                                                       ]
                [     sin(theta)            r cos(theta)                   0            ]
          (%i5) cmetric();
          (%o5)                                false
          (%i6) lg:trigsimp(lg);
                                     [ 1  0         0        ]
                                     [                       ]
                                     [     2                 ]
          (%o6)                      [ 0  r         0        ]
                                     [                       ]
                                     [         2    2        ]
                                     [ 0  0   r  cos (theta) ]

     The optional argument <extra_arg> can be any one of the following:

     `cylindrical' tells `ct_coordsys' to attach an additional
     cylindrical coordinate.

     `minkowski' tells `ct_coordsys' to attach an additional coordinate
     with negative metric signature.

     `all' tells `ct_coordsys' to call `cmetric' and `christof(false)'
     after setting up the metric.

     If the global variable `verbose' is set to `true', `ct_coordsys'
     displays the values of `dim', `ct_coords', and either `lg' or
     `lfg' and `fri', depending on the value of `cframe_flag'.


 -- Function: init_ctensor ()
     Initializes the `ctensor' package.

     The `init_ctensor' function reinitializes the `ctensor' package.
     It removes all arrays and matrices used by `ctensor', resets all
     flags, resets `dim' to 4, and resets the frame metric to the
     Lorentz-frame.


29.2.2 The tensors of curved space
----------------------------------

The main purpose of the `ctensor' package is to compute the tensors of
curved space(time), most notably the tensors used in general relativity.

   When a metric base is used, `ctensor' can compute the following
tensors:


      lg  -- ug
        \      \
         lcs -- mcs -- ric -- uric
                   \      \       \
                    \      tracer - ein -- lein
                     \
                      riem -- lriem -- weyl
                          \
                           uriem

   `ctensor' can also work using moving frames. When `cframe_flag' is
set to `true', the following tensors can be calculated:


      lfg -- ufg
          \
      fri -- fr -- lcs -- mcs -- lriem -- ric -- uric
           \                       |  \      \       \
            lg -- ug               |   weyl   tracer - ein -- lein
                                   |\
                                   | riem
                                   |
                                   \uriem

 -- Function: christof (<dis>)
     A function in the `ctensor' (component tensor) package.  It
     computes the Christoffel symbols of both kinds.  The argument
     <dis> determines which results are to be immediately displayed.
     The Christoffel symbols of the first and second kinds are stored
     in the arrays `lcs[i,j,k]' and `mcs[i,j,k]' respectively and
     defined to be symmetric in the first two indices. If the argument
     to `christof' is `lcs' or `mcs' then the unique non-zero values of
     `lcs[i,j,k]' or `mcs[i,j,k]', respectively, will be displayed. If
     the argument is `all' then the unique non-zero values of
     `lcs[i,j,k]' and `mcs[i,j,k]' will be displayed.  If the argument
     is `false' then the display of the elements will not occur. The
     array elements `mcs[i,j,k]' are defined in such a manner that the
     final index is contravariant.


 -- Function: ricci (<dis>)
     A function in the `ctensor' (component tensor) package. `ricci'
     computes the covariant (symmetric) components `ric[i,j]' of the
     Ricci tensor.  If the argument <dis> is `true', then the non-zero
     components are displayed.


 -- Function: uricci (<dis>)
     This function first computes the covariant components `ric[i,j]'
     of the Ricci tensor.  Then the mixed Ricci tensor is computed
     using the contravariant metric tensor.  If the value of the
     argument <dis> is `true', then these mixed components, `uric[i,j]'
     (the index `i' is covariant and the index `j' is contravariant),
     will be displayed directly.  Otherwise, `ricci(false)' will simply
     compute the entries of the array `uric[i,j]' without displaying
     the results.


 -- Function: scurvature ()
     Returns the scalar curvature (obtained by contracting the Ricci
     tensor) of the Riemannian manifold with the given metric.


 -- Function: einstein (<dis>)
     A function in the `ctensor' (component tensor) package.
     `einstein' computes the mixed Einstein tensor after the
     Christoffel symbols and Ricci tensor have been obtained (with the
     functions `christof' and `ricci').  If the argument <dis> is
     `true', then the non-zero values of the mixed Einstein tensor
     `ein[i,j]' will be displayed where `j' is the contravariant index.
     The variable `rateinstein' will cause the rational simplification
     on these components. If `ratfac' is `true' then the components will
     also be factored.


 -- Function: leinstein (<dis>)
     Covariant Einstein-tensor. `leinstein' stores the values of the
     covariant Einstein tensor in the array `lein'. The covariant
     Einstein-tensor is computed from the mixed Einstein tensor `ein'
     by multiplying it with the metric tensor. If the argument <dis> is
     `true', then the non-zero values of the covariant Einstein tensor
     are displayed.


 -- Function: riemann (<dis>)
     A function in the `ctensor' (component tensor) package.  `riemann'
     computes the Riemann curvature tensor from the given metric and
     the corresponding Christoffel symbols. The following index
     conventions are used:

                          l      _l       _l       _l   _m    _l   _m
           R[i,j,k,l] =  R    = |      - |      + |    |   - |    |
                          ijk     ij,k     ik,j     mk   ij    mj   ik

     This notation is consistent with the notation used by the `itensor'
     package and its `icurvature' function.  If the optional argument
     <dis> is `true', the non-zero components `riem[i,j,k,l]' will be
     displayed.  As with the Einstein tensor, various switches set by
     the user control the simplification of the components of the
     Riemann tensor.  If `ratriemann' is `true', then rational
     simplification will be done. If `ratfac' is `true' then each of
     the components will also be factored.

     If the variable `cframe_flag' is `false', the Riemann tensor is
     computed directly from the Christoffel-symbols. If `cframe_flag' is
     `true', the covariant Riemann-tensor is computed first from the
     frame field coefficients.


 -- Function: lriemann (<dis>)
     Covariant Riemann-tensor (`lriem[]').

     Computes the covariant Riemann-tensor as the array `lriem'. If the
     argument <dis> is `true', unique nonzero values are displayed.

     If the variable `cframe_flag' is `true', the covariant Riemann
     tensor is computed directly from the frame field coefficients.
     Otherwise, the (3,1) Riemann tensor is computed first.

     For information on index ordering, see `riemann'.


 -- Function: uriemann (<dis>)
     Computes the contravariant components of the Riemann curvature
     tensor as array elements `uriem[i,j,k,l]'.  These are displayed if
     <dis> is `true'.


 -- Function: rinvariant ()
     Forms the Kretchmann-invariant (`kinvariant') obtained by
     contracting the tensors

          lriem[i,j,k,l]*uriem[i,j,k,l].

     This object is not automatically simplified since it can be very
     large.


 -- Function: weyl (<dis>)
     Computes the Weyl conformal tensor.  If the argument <dis> is
     `true', the non-zero components `weyl[i,j,k,l]' will be displayed
     to the user.  Otherwise, these components will simply be computed
     and stored.  If the switch `ratweyl' is set to `true', then the
     components will be rationally simplified; if `ratfac' is `true'
     then the results will be factored as well.


29.2.3 Taylor series expansion
------------------------------

The `ctensor' package has the ability to truncate results by assuming
that they are Taylor-series approximations. This behavior is controlled
by the `ctayswitch' variable; when set to true, `ctensor' makes use
internally of the function `ctaylor' when simplifying results.

   The `ctaylor' function is invoked by the following `ctensor'
functions:


         Function     Comments
         ---------------------------------
         christof()   For mcs only
         ricci()
         uricci()
         einstein()
         riemann()
         weyl()
         checkdiv()

 -- Function: ctaylor ()
     The `ctaylor' function truncates its argument by converting it to
     a Taylor-series using `taylor', and then calling `ratdisrep'. This
     has the combined effect of dropping terms higher order in the
     expansion variable `ctayvar'. The order of terms that should be
     dropped is defined by `ctaypov'; the point around which the series
     expansion is carried out is specified in `ctaypt'.

     As an example, consider a simple metric that is a perturbation of
     the Minkowski metric. Without further restrictions, even a diagonal
     metric produces expressions for the Einstein tensor that are far
     too complex:


          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) ratfac:true;
          (%o2)                                true
          (%i3) derivabbrev:true;
          (%o3)                                true
          (%i4) ct_coords:[t,r,theta,phi];
          (%o4)                         [t, r, theta, phi]
          (%i5) lg:matrix([-1,0,0,0],[0,1,0,0],[0,0,r^2,0],[0,0,0,r^2*sin(theta)^2]);
                                  [ - 1  0  0         0        ]
                                  [                            ]
                                  [  0   1  0         0        ]
                                  [                            ]
          (%o5)                   [          2                 ]
                                  [  0   0  r         0        ]
                                  [                            ]
                                  [              2    2        ]
                                  [  0   0  0   r  sin (theta) ]
          (%i6) h:matrix([h11,0,0,0],[0,h22,0,0],[0,0,h33,0],[0,0,0,h44]);
                                      [ h11   0    0    0  ]
                                      [                    ]
                                      [  0   h22   0    0  ]
          (%o6)                       [                    ]
                                      [  0    0   h33   0  ]
                                      [                    ]
                                      [  0    0    0   h44 ]
          (%i7) depends(l,r);
          (%o7)                               [l(r)]
          (%i8) lg:lg+l*h;
                   [ h11 l - 1      0          0                 0            ]
                   [                                                          ]
                   [     0      h22 l + 1      0                 0            ]
                   [                                                          ]
          (%o8)    [                        2                                 ]
                   [     0          0      r  + h33 l            0            ]
                   [                                                          ]
                   [                                    2    2                ]
                   [     0          0          0       r  sin (theta) + h44 l ]
          (%i9) cmetric(false);
          (%o9)                                done
          (%i10) einstein(false);
          (%o10)                               done
          (%i11) ntermst(ein);
          [[1, 1], 62]
          [[1, 2], 0]
          [[1, 3], 0]
          [[1, 4], 0]
          [[2, 1], 0]
          [[2, 2], 24]
          [[2, 3], 0]
          [[2, 4], 0]
          [[3, 1], 0]
          [[3, 2], 0]
          [[3, 3], 46]
          [[3, 4], 0]
          [[4, 1], 0]
          [[4, 2], 0]
          [[4, 3], 0]
          [[4, 4], 46]
          (%o12)                               done

     However, if we recompute this example as an approximation that is
     linear in the variable `l', we get much simpler expressions:


          (%i14) ctayswitch:true;
          (%o14)                               true
          (%i15) ctayvar:l;
          (%o15)                                 l
          (%i16) ctaypov:1;
          (%o16)                                 1
          (%i17) ctaypt:0;
          (%o17)                                 0
          (%i18) christof(false);
          (%o18)                               done
          (%i19) ricci(false);
          (%o19)                               done
          (%i20) einstein(false);
          (%o20)                               done
          (%i21) ntermst(ein);
          [[1, 1], 6]
          [[1, 2], 0]
          [[1, 3], 0]
          [[1, 4], 0]
          [[2, 1], 0]
          [[2, 2], 13]
          [[2, 3], 2]
          [[2, 4], 0]
          [[3, 1], 0]
          [[3, 2], 2]
          [[3, 3], 9]
          [[3, 4], 0]
          [[4, 1], 0]
          [[4, 2], 0]
          [[4, 3], 0]
          [[4, 4], 9]
          (%o21)                               done
          (%i22) ratsimp(ein[1,1]);
                                   2      2  4               2     2
          (%o22) - (((h11 h22 - h11 ) (l )  r  - 2 h33 l    r ) sin (theta)
                                        r               r r

                                          2               2      4    2
                            - 2 h44 l    r  - h33 h44 (l ) )/(4 r  sin (theta))
                                     r r                r

     This capability can be useful, for instance, when working in the
     weak field limit far from a gravitational source.


29.2.4 Frame fields
-------------------

When the variable `cframe_flag' is set to true, the `ctensor' package
performs its calculations using a moving frame.

 -- Function: frame_bracket (<fr>, <fri>, <diagframe>)
     The frame bracket (`fb[]').

     Computes the frame bracket according to the following definition:

             c          c         c        d     e
          ifb   = ( ifri    - ifri    ) ifr   ifr
             ab         d,e       e,d      a     b


29.2.5 Algebraic classification
-------------------------------

A new feature (as of November, 2004) of `ctensor' is its ability to
compute the Petrov classification of a 4-dimensional spacetime metric.
For a demonstration of this capability, see the file
`share/tensor/petrov.dem'.

 -- Function: nptetrad ()
     Computes a Newman-Penrose null tetrad (`np') and its raised-index
     counterpart (`npi'). See `petrov' for an example.

     The null tetrad is constructed on the assumption that a
     four-diemensional orthonormal frame metric with metric signature
     (-,+,+,+) is being used.  The components of the null tetrad are
     related to the inverse frame matrix as follows:


          np  = (fri  + fri ) / sqrt(2)
            1       1      2

          np  = (fri  - fri ) / sqrt(2)
            2       1      2

          np  = (fri  + %i fri ) / sqrt(2)
            3       3         4

          np  = (fri  - %i fri ) / sqrt(2)
            4       3         4


 -- Function: psi (<dis>)
     Computes the five Newman-Penrose coefficients `psi[0]'...`psi[4]'.
     If `psi' is set to `true', the coefficients are displayed.  See
     `petrov' for an example.

     These coefficients are computed from the Weyl-tensor in a
     coordinate base.  If a frame base is used, the Weyl-tensor is
     first converted to a coordinate base, which can be a
     computationally expensive procedure. For this reason, in some
     cases it may be more advantageous to use a coordinate base in the
     first place before the Weyl tensor is computed. Note however, that
     constructing a Newman-Penrose null tetrad requires a frame base.
     Therefore, a meaningful computation sequence may begin with a
     frame base, which is then used to compute `lg' (computed
     automatically by `cmetric' and then `ug'. At this point, you can
     switch back to a coordinate base by setting `cframe_flag' to false
     before beginning to compute the Christoffel symbols. Changing to a
     frame base at a later stage could yield inconsistent results, as
     you may end up with a mixed bag of tensors, some computed in a
     frame base, some in a coordinate base, with no means to
     distinguish between the two.


 -- Function: petrov ()
     Computes the Petrov classification of the metric characterized by
     `psi[0]'...`psi[4]'.

     For example, the following demonstrates how to obtain the
     Petrov-classification of the Kerr metric:

          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) (cframe_flag:true,gcd:spmod,ctrgsimp:true,ratfac:true);
          (%o2)                                true
          (%i3) ct_coordsys(exteriorschwarzschild,all);
          (%o3)                                done
          (%i4) ug:invert(lg)$
          (%i5) weyl(false);
          (%o5)                                done
          (%i6) nptetrad(true);
          (%t6) np =

                 [  sqrt(r - 2 m)           sqrt(r)                                     ]
                 [ ---------------   ---------------------      0             0         ]
                 [ sqrt(2) sqrt(r)   sqrt(2) sqrt(r - 2 m)                              ]
                 [                                                                      ]
                 [  sqrt(r - 2 m)            sqrt(r)                                    ]
                 [ ---------------  - ---------------------     0             0         ]
                 [ sqrt(2) sqrt(r)    sqrt(2) sqrt(r - 2 m)                             ]
                 [                                                                      ]
                 [                                              r      %i r sin(theta)  ]
                 [        0                    0             -------   ---------------  ]
                 [                                           sqrt(2)       sqrt(2)      ]
                 [                                                                      ]
                 [                                              r       %i r sin(theta) ]
                 [        0                    0             -------  - --------------- ]
                 [                                           sqrt(2)        sqrt(2)     ]

                                       sqrt(r)          sqrt(r - 2 m)
          (%t7) npi = matrix([- ---------------------, ---------------, 0, 0],
                                sqrt(2) sqrt(r - 2 m)  sqrt(2) sqrt(r)

                    sqrt(r)            sqrt(r - 2 m)
          [- ---------------------, - ---------------, 0, 0],
             sqrt(2) sqrt(r - 2 m)    sqrt(2) sqrt(r)

                     1               %i
          [0, 0, ---------, --------------------],
                 sqrt(2) r  sqrt(2) r sin(theta)

                     1                 %i
          [0, 0, ---------, - --------------------])
                 sqrt(2) r    sqrt(2) r sin(theta)

          (%o7)                                done
          (%i7) psi(true);
          (%t8)                              psi  = 0
                                                0

          (%t9)                              psi  = 0
                                                1

                                                    m
          (%t10)                             psi  = --
                                                2    3
                                                    r

          (%t11)                             psi  = 0
                                                3

          (%t12)                             psi  = 0
                                                4
          (%o12)                               done
          (%i12) petrov();
          (%o12)                                 D

     The Petrov classification function is based on the algorithm
     published in "Classifying geometries in general relativity: III
     Classification in practice" by Pollney, Skea, and d'Inverno,
     Class. Quant. Grav. 17 2885-2902 (2000).  Except for some simple
     test cases, the implementation is untested as of December 19,
     2004, and is likely to contain errors.


29.2.6 Torsion and nonmetricity
-------------------------------

`ctensor' has the ability to compute and include torsion and
nonmetricity coefficients in the connection coefficients.

   The torsion coefficients are calculated from a user-supplied tensor
`tr', which should be a rank (2,1) tensor. From this, the torsion
coefficients `kt' are computed according to the following formulae:


                   m          m      m
            - g  tr   - g   tr   - tr   g
               im  kj    jm   ki     ij  km
     kt   = -------------------------------
       ijk                 2


       k     km
     kt   = g   kt
       ij         ijm

   Note that only the mixed-index tensor is calculated and stored in the
array `kt'.

   The nonmetricity coefficients are calculated from the user-supplied
nonmetricity vector `nm'. From this, the nonmetricity coefficients
`nmc' are computed as follows:


                  k    k        km
            -nm  D  - D  nm  + g   nm  g
        k      i  j    i   j         m  ij
     nmc  = ------------------------------
        ij                2

   where D stands for the Kronecker-delta.

   When `ctorsion_flag' is set to `true', the values of `kt' are
substracted from the mixed-indexed connection coefficients computed by
`christof' and stored in `mcs'. Similarly, if `cnonmet_flag' is set to
`true', the values of `nmc' are substracted from the mixed-indexed
connection coefficients.

   If necessary, `christof' calls the functions `contortion' and
`nonmetricity' in order to compute `kt' and `nm'.

 -- Function: contortion (<tr>)
     Computes the (2,1) contortion coefficients from the torsion tensor
     <tr>.


 -- Function: nonmetricity (<nm>)
     Computes the (2,1) nonmetricity coefficients from the nonmetricity
     vector <nm>.


29.2.7 Miscellaneous features
-----------------------------

 -- Function: ctransform (<M>)
     A function in the `ctensor' (component tensor) package which will
     perform a coordinate transformation upon an arbitrary square
     symmetric matrix <M>. The user must input the functions which
     define the transformation.  (Formerly called `transform'.)


 -- Function: findde (<A>, <n>)
     returns a list of the unique differential equations (expressions)
     corresponding to the elements of the <n> dimensional square array
     <A>. Presently, <n> may be 2 or 3. `deindex' is a global list
     containing the indices of <A> corresponding to these unique
     differential equations. For the Einstein tensor (`ein'), which is
     a two dimensional array, if computed for the metric in the example
     below, `findde' gives the following independent differential
     equations:

          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) derivabbrev:true;
          (%o2)                                true
          (%i3) dim:4;
          (%o3)                                  4
          (%i4) lg:matrix([a,0,0,0],[0,x^2,0,0],[0,0,x^2*sin(y)^2,0],[0,0,0,-d]);
                                    [ a  0       0        0  ]
                                    [                        ]
                                    [     2                  ]
                                    [ 0  x       0        0  ]
          (%o4)                     [                        ]
                                    [         2    2         ]
                                    [ 0  0   x  sin (y)   0  ]
                                    [                        ]
                                    [ 0  0       0       - d ]
          (%i5) depends([a,d],x);
          (%o5)                            [a(x), d(x)]
          (%i6) ct_coords:[x,y,z,t];
          (%o6)                            [x, y, z, t]
          (%i7) cmetric();
          (%o7)                                done
          (%i8) einstein(false);
          (%o8)                                done
          (%i9) findde(ein,2);
                                                      2
          (%o9) [d  x - a d + d, 2 a d d    x - a (d )  x - a  d d  x + 2 a d d
                  x                     x x         x        x    x            x

                                                                  2          2
                                                          - 2 a  d , a  x + a  - a]
                                                               x      x
          (%i10) deindex;
          (%o10)                     [[1, 1], [2, 2], [4, 4]]


 -- Function: cograd ()
     Computes the covariant gradient of a scalar function allowing the
     user to choose the corresponding vector name as the example under
     `contragrad' illustrates.

 -- Function: contragrad ()
     Computes the contravariant gradient of a scalar function allowing
     the user to choose the corresponding vector name as the example
     below for the Schwarzschild metric illustrates:


          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) derivabbrev:true;
          (%o2)                                true
          (%i3) ct_coordsys(exteriorschwarzschild,all);
          (%o3)                                done
          (%i4) depends(f,r);
          (%o4)                               [f(r)]
          (%i5) cograd(f,g1);
          (%o5)                                done
          (%i6) listarray(g1);
          (%o6)                            [0, f , 0, 0]
                                                r
          (%i7) contragrad(f,g2);
          (%o7)                                done
          (%i8) listarray(g2);
                                         f  r - 2 f  m
                                          r        r
          (%o8)                      [0, -------------, 0, 0]
                                               r


 -- Function: dscalar ()
     computes the tensor d'Alembertian of the scalar function once
     dependencies have been declared upon the function. For example:

          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) derivabbrev:true;
          (%o2)                                true
          (%i3) ct_coordsys(exteriorschwarzschild,all);
          (%o3)                                done
          (%i4) depends(p,r);
          (%o4)                               [p(r)]
          (%i5) factor(dscalar(p));
                                    2
                              p    r  - 2 m p    r + 2 p  r - 2 m p
                               r r           r r        r          r
          (%o5)               --------------------------------------
                                                 2
                                                r


 -- Function: checkdiv ()
     computes the covariant divergence of the mixed second rank tensor
     (whose first index must be covariant) by printing the
     corresponding n components of the vector field (the divergence)
     where n = `dim'. If the argument to the function is `g' then the
     divergence of the Einstein tensor will be formed and must be zero.
     In addition, the divergence (vector) is given the array name `div'.

 -- Function: cgeodesic (<dis>)
     A function in the `ctensor' (component tensor) package.
     `cgeodesic' computes the geodesic equations of motion for a given
     metric.  They are stored in the array `geod[i]'.  If the argument
     <dis> is `true' then these equations are displayed.


 -- Function: bdvac (<f>)
     generates the covariant components of the vacuum field equations of
     the Brans- Dicke gravitational theory. The scalar field is
     specified by the argument <f>, which should be a (quoted) function
     name with functional dependencies, e.g., `'p(x)'.

     The components of the second rank covariant field tensor are
     represented by the array `bd'.


 -- Function: invariant1 ()
     generates the mixed Euler- Lagrange tensor (field equations) for
     the invariant density of R^2. The field equations are the
     components of an array named `inv1'.


 -- Function: invariant2 ()
     *** NOT YET IMPLEMENTED ***

     generates the mixed Euler- Lagrange tensor (field equations) for
     the invariant density of `ric[i,j]*uriem[i,j]'. The field
     equations are the components of an array named `inv2'.


 -- Function: bimetric ()
     *** NOT YET IMPLEMENTED ***

     generates the field equations of Rosen's bimetric theory. The field
     equations are the components of an array named `rosen'.


29.2.8 Utility functions
------------------------

 -- Function: diagmatrixp (<M>)
     Returns `true' if <M> is a diagonal matrix or (2D) array.


 -- Function: symmetricp (<M>)
     Returns `true' if <M> is a symmetric matrix or (2D) array.


 -- Function: ntermst (<f>)
     gives the user a quick picture of the "size" of the doubly
     subscripted tensor (array) <f>.  It prints two element lists where
     the second element corresponds to NTERMS of the components
     specified by the first elements.  In this way, it is possible to
     quickly find the non-zero expressions and attempt simplification.


 -- Function: cdisplay (<ten>)
     displays all the elements of the tensor <ten>, as represented by a
     multidimensional array. Tensors of rank 0 and 1, as well as other
     types of variables, are displayed as with `ldisplay'. Tensors of
     rank 2 are displayed as 2-dimensional matrices, while tensors of
     higher rank are displayed as a list of 2-dimensional matrices. For
     instance, the Riemann-tensor of the Schwarzschild metric can be
     viewed as:

          (%i1) load(ctensor);
          (%o1)       /share/tensor/ctensor.mac
          (%i2) ratfac:true;
          (%o2)                                true
          (%i3) ct_coordsys(exteriorschwarzschild,all);
          (%o3)                                done
          (%i4) riemann(false);
          (%o4)                                done
          (%i5) cdisplay(riem);
                         [ 0               0                    0            0      ]
                         [                                                          ]
                         [                              2                           ]
                         [      3 m (r - 2 m)   m    2 m                            ]
                         [ 0  - ------------- + -- - ----       0            0      ]
                         [            4          3     4                            ]
                         [           r          r     r                             ]
                         [                                                          ]
              riem     = [                                 m (r - 2 m)              ]
                  1, 1   [ 0               0               -----------       0      ]
                         [                                      4                   ]
                         [                                     r                    ]
                         [                                                          ]
                         [                                              m (r - 2 m) ]
                         [ 0               0                    0       ----------- ]
                         [                                                   4      ]
                         [                                                  r       ]

                                          [    2 m (r - 2 m)       ]
                                          [ 0  -------------  0  0 ]
                                          [          4             ]
                                          [         r              ]
                               riem     = [                        ]
                                   1, 2   [ 0        0        0  0 ]
                                          [                        ]
                                          [ 0        0        0  0 ]
                                          [                        ]
                                          [ 0        0        0  0 ]

                                          [         m (r - 2 m)    ]
                                          [ 0  0  - -----------  0 ]
                                          [              4         ]
                                          [             r          ]
                               riem     = [                        ]
                                   1, 3   [ 0  0        0        0 ]
                                          [                        ]
                                          [ 0  0        0        0 ]
                                          [                        ]
                                          [ 0  0        0        0 ]

                                          [            m (r - 2 m) ]
                                          [ 0  0  0  - ----------- ]
                                          [                 4      ]
                                          [                r       ]
                               riem     = [                        ]
                                   1, 4   [ 0  0  0        0       ]
                                          [                        ]
                                          [ 0  0  0        0       ]
                                          [                        ]
                                          [ 0  0  0        0       ]

                                         [       0         0  0  0 ]
                                         [                         ]
                                         [       2 m               ]
                                         [ - ------------  0  0  0 ]
                              riem     = [    2                    ]
                                  2, 1   [   r  (r - 2 m)          ]
                                         [                         ]
                                         [       0         0  0  0 ]
                                         [                         ]
                                         [       0         0  0  0 ]

                             [     2 m                                         ]
                             [ ------------  0        0               0        ]
                             [  2                                              ]
                             [ r  (r - 2 m)                                    ]
                             [                                                 ]
                             [      0        0        0               0        ]
                             [                                                 ]
                  riem     = [                         m                       ]
                      2, 2   [      0        0  - ------------        0        ]
                             [                     2                           ]
                             [                    r  (r - 2 m)                 ]
                             [                                                 ]
                             [                                         m       ]
                             [      0        0        0         - ------------ ]
                             [                                     2           ]
                             [                                    r  (r - 2 m) ]

                                          [ 0  0       0        0 ]
                                          [                       ]
                                          [            m          ]
                                          [ 0  0  ------------  0 ]
                               riem     = [        2              ]
                                   2, 3   [       r  (r - 2 m)    ]
                                          [                       ]
                                          [ 0  0       0        0 ]
                                          [                       ]
                                          [ 0  0       0        0 ]

                                          [ 0  0  0       0       ]
                                          [                       ]
                                          [               m       ]
                                          [ 0  0  0  ------------ ]
                               riem     = [           2           ]
                                   2, 4   [          r  (r - 2 m) ]
                                          [                       ]
                                          [ 0  0  0       0       ]
                                          [                       ]
                                          [ 0  0  0       0       ]

                                                [ 0  0  0  0 ]
                                                [            ]
                                                [ 0  0  0  0 ]
                                                [            ]
                                     riem     = [ m          ]
                                         3, 1   [ -  0  0  0 ]
                                                [ r          ]
                                                [            ]
                                                [ 0  0  0  0 ]

                                                [ 0  0  0  0 ]
                                                [            ]
                                                [ 0  0  0  0 ]
                                                [            ]
                                     riem     = [    m       ]
                                         3, 2   [ 0  -  0  0 ]
                                                [    r       ]
                                                [            ]
                                                [ 0  0  0  0 ]

                                         [   m                      ]
                                         [ - -   0   0       0      ]
                                         [   r                      ]
                                         [                          ]
                                         [        m                 ]
                                         [  0   - -  0       0      ]
                              riem     = [        r                 ]
                                  3, 3   [                          ]
                                         [  0    0   0       0      ]
                                         [                          ]
                                         [              2 m - r     ]
                                         [  0    0   0  ------- + 1 ]
                                         [                 r        ]

                                              [ 0  0  0    0   ]
                                              [                ]
                                              [ 0  0  0    0   ]
                                              [                ]
                                   riem     = [            2 m ]
                                       3, 4   [ 0  0  0  - --- ]
                                              [             r  ]
                                              [                ]
                                              [ 0  0  0    0   ]

                                          [       0        0  0  0 ]
                                          [                        ]
                                          [       0        0  0  0 ]
                                          [                        ]
                               riem     = [       0        0  0  0 ]
                                   4, 1   [                        ]
                                          [      2                 ]
                                          [ m sin (theta)          ]
                                          [ -------------  0  0  0 ]
                                          [       r                ]

                                          [ 0        0        0  0 ]
                                          [                        ]
                                          [ 0        0        0  0 ]
                                          [                        ]
                               riem     = [ 0        0        0  0 ]
                                   4, 2   [                        ]
                                          [         2              ]
                                          [    m sin (theta)       ]
                                          [ 0  -------------  0  0 ]
                                          [          r             ]

                                        [ 0  0          0          0 ]
                                        [                            ]
                                        [ 0  0          0          0 ]
                                        [                            ]
                             riem     = [ 0  0          0          0 ]
                                 4, 3   [                            ]
                                        [                2           ]
                                        [         2 m sin (theta)    ]
                                        [ 0  0  - ---------------  0 ]
                                        [                r           ]

                           [        2                                             ]
                           [   m sin (theta)                                      ]
                           [ - -------------         0                0         0 ]
                           [         r                                            ]
                           [                                                      ]
                           [                         2                            ]
                           [                    m sin (theta)                     ]
                riem     = [        0         - -------------         0         0 ]
                    4, 4   [                          r                           ]
                           [                                                      ]
                           [                                          2           ]
                           [                                   2 m sin (theta)    ]
                           [        0                0         ---------------  0 ]
                           [                                          r           ]
                           [                                                      ]
                           [        0                0                0         0 ]

          (%o5)                                done

 -- Function: deleten (<L>, <n>)
     Returns a new list consisting of <L> with the <n>'th element
     deleted.

29.2.9 Variables used by `ctensor'
----------------------------------

 -- Option variable: dim
     Default value: 4

     An option in the `ctensor' (component tensor) package.  `dim' is
     the dimension of the manifold with the default 4. The command
     `dim: n' will reset the dimension to any other value `n'.


 -- Option variable: diagmetric
     Default value: `false'

     An option in the `ctensor' (component tensor) package.  If
     `diagmetric' is `true' special routines compute all geometrical
     objects (which contain the metric tensor explicitly) by taking
     into consideration the diagonality of the metric. Reduced run
     times will, of course, result. Note: this option is set
     automatically by `csetup' if a diagonal metric is specified.


 -- Option variable: ctrgsimp
     Causes trigonometric simplifications to be used when tensors are
     computed. Presently, `ctrgsimp' affects only computations
     involving a moving frame.


 -- Option variable: cframe_flag
     Causes computations to be performed relative to a moving frame as
     opposed to a holonomic metric. The frame is defined by the inverse
     frame array `fri' and the frame metric `lfg'. For computations
     using a Cartesian frame, `lfg' should be the unit matrix of the
     appropriate dimension; for computations in a Lorentz frame, `lfg'
     should have the appropriate signature.


 -- Option variable: ctorsion_flag
     Causes the contortion tensor to be included in the computation of
     the connection coefficients. The contortion tensor itself is
     computed by `contortion' from the user-supplied tensor `tr'.


 -- Option variable: cnonmet_flag
     Causes the nonmetricity coefficients to be included in the
     computation of the connection coefficients. The nonmetricity
     coefficients are computed from the user-supplied nonmetricity
     vector `nm' by the function `nonmetricity'.


 -- Option variable: ctayswitch
     If set to `true', causes some `ctensor' computations to be carried
     out using Taylor-series expansions. Presently, `christof', `ricci',
     `uricci', `einstein', and `weyl' take into account this setting.


 -- Option variable: ctayvar
     Variable used for Taylor-series expansion if `ctayswitch' is set to
     `true'.


 -- Option variable: ctaypov
     Maximum power used in Taylor-series expansion when `ctayswitch' is
     set to `true'.


 -- Option variable: ctaypt
     Point around which Taylor-series expansion is carried out when
     `ctayswitch' is set to `true'.


 -- System variable: gdet
     The determinant of the metric tensor `lg'. Computed by `cmetric'
     when `cframe_flag' is set to `false'.


 -- Option variable: ratchristof
     Causes rational simplification to be applied by `christof'.


 -- Option variable: rateinstein
     Default value: `true'

     If `true' rational simplification will be performed on the
     non-zero components of Einstein tensors; if `ratfac' is `true'
     then the components will also be factored.


 -- Option variable: ratriemann
     Default value: `true'

     One of the switches which controls simplification of Riemann
     tensors; if `true', then rational simplification will be done; if
     `ratfac' is `true' then each of the components will also be
     factored.


 -- Option variable: ratweyl
     Default value: `true'

     If `true', this switch causes the `weyl' function to apply
     rational simplification to the values of the Weyl tensor. If
     `ratfac' is `true', then the components will also be factored.

 -- Variable: lfg
     The covariant frame metric. By default, it is initialized to the
     4-dimensional Lorentz frame with signature (+,+,+,-). Used when
     `cframe_flag' is `true'.

 -- Variable: ufg
     The inverse frame metric. Computed from `lfg' when `cmetric' is
     called while `cframe_flag' is set to `true'.

 -- Variable: riem
     The (3,1) Riemann tensor. Computed when the function `riemann' is
     invoked. For information about index ordering, see the description
     of `riemann'.

     if `cframe_flag' is `true', `riem' is computed from the covariant
     Riemann-tensor `lriem'.


 -- Variable: lriem
     The covariant Riemann tensor. Computed by `lriemann'.


 -- Variable: uriem
     The contravariant Riemann tensor. Computed by `uriemann'.


 -- Variable: ric
     The mixed Ricci-tensor. Computed by `ricci'.


 -- Variable: uric
     The contravariant Ricci-tensor. Computed by `uricci'.


 -- Variable: lg
     The metric tensor. This tensor must be specified (as a `dim' by
     `dim' matrix) before other computations can be performed.


 -- Variable: ug
     The inverse of the metric tensor. Computed by `cmetric'.


 -- Variable: weyl
     The Weyl tensor. Computed by `weyl'.


 -- Variable: fb
     Frame bracket coefficients, as computed by `frame_bracket'.


 -- Variable: kinvariant
     The Kretchmann invariant. Computed by `rinvariant'.


 -- Variable: np
     A Newman-Penrose null tetrad. Computed by `nptetrad'.


 -- Variable: npi
     The raised-index Newman-Penrose null tetrad. Computed by
     `nptetrad'.  Defined as `ug.np'. The product `np.transpose(npi)'
     is constant:

          (%i39) trigsimp(np.transpose(npi));
                                        [  0   - 1  0  0 ]
                                        [                ]
                                        [ - 1   0   0  0 ]
          (%o39)                        [                ]
                                        [  0    0   0  1 ]
                                        [                ]
                                        [  0    0   1  0 ]


 -- Variable: tr
     User-supplied rank-3 tensor representing torsion. Used by
     `contortion'.

 -- Variable: kt
     The contortion tensor, computed from `tr' by `contortion'.

 -- Variable: nm
     User-supplied nonmetricity vector. Used by `nonmetricity'.

 -- Variable: nmc
     The nonmetricity coefficients, computed from `nm' by
     `nonmetricity'.


 -- System variable: tensorkill
     Variable indicating if the tensor package has been initialized.
     Set and used by `csetup', reset by `init_ctensor'.


 -- Option variable: ct_coords
     Default value: `[]'

     An option in the `ctensor' (component tensor) package.
     `ct_coords' contains a list of coordinates.  While normally
     defined when the function `csetup' is called, one may redefine the
     coordinates with the assignment `ct_coords: [j1, j2, ..., jn]'
     where the j's are the new coordinate names.  See also `csetup'.


29.2.10 Reserved names
----------------------

The following names are used internally by the `ctensor' package and
should not be redefined:

       Name         Description
       ---------------------------------------
       _lg()        Evaluates to lfg if frame metric used, lg otherwise
       _ug()        Evaluates to ufg if frame metric used, ug otherwise
       cleanup()    Removes items drom the deindex list
       contract4()  Used by psi()
       filemet()    Used by csetup() when reading the metric from a file
       findde1()    Used by findde()
       findde2()    Used by findde()
       findde3()    Used by findde()
       kdelt()      Kronecker-delta (not generalized)
       newmet()     Used by csetup() for setting up a metric interactively
       setflags()   Used by init_ctensor()
       readvalue()
       resimp()
       sermet()     Used by csetup() for entering a metric as Taylor-series
       txyzsum()
       tmetric()    Frame metric, used by cmetric() when cframe_flag:true
       triemann()   Riemann-tensor in frame base, used when cframe_flag:true
       tricci()     Ricci-tensor in frame base, used when cframe_flag:true
       trrc()       Ricci rotation coefficients, used by christof()
       yesp()

29.2.11 Changes
---------------

In November, 2004, the `ctensor' package was extensively rewritten.
Many functions and variables have been renamed in order to make the
package compatible with the commercial version of Macsyma.

       New Name     Old Name        Description
       --------------------------------------------------------------------------
       ctaylor()    DLGTAYLOR()     Taylor-series expansion of an expression
       lgeod[]      EM              Geodesic equations
       ein[]        G[]             Mixed Einstein-tensor
       ric[]        LR[]            Mixed Ricci-tensor
       ricci()      LRICCICOM()     Compute the mixed Ricci-tensor
       ctaypov      MINP            Maximum power in Taylor-series expansion
       cgeodesic()  MOTION          Compute geodesic equations
       ct_coords    OMEGA           Metric coordinates
       ctayvar      PARAM           Taylor-series expansion variable
       lriem[]      R[]             Covariant Riemann-tensor
       uriemann()   RAISERIEMANN()  Compute the contravariant Riemann-tensor
       ratriemann   RATRIEMAN       Rational simplification of the Riemann-tensor
       uric[]       RICCI[]         Contravariant Ricci-tensor
       uricci()     RICCICOM()      Compute the contravariant Ricci-tensor
       cmetric()    SETMETRIC()     Set up the metric
       ctaypt       TAYPT           Point for Taylor-series expansion
       ctayswitch   TAYSWITCH       Taylor-series setting switch
       csetup()     TSETUP()        Start interactive setup session
       ctransform() TTRANSFORM()    Interactive coordinate transformation
       uriem[]      UR[]            Contravariant Riemann-tensor
       weyl[]       W[]             (3,1) Weyl-tensor


File: maxima.info,  Node: atensor,  Next: Series,  Prev: ctensor,  Up: Top

30 atensor
**********

* Menu:

* Introduction to atensor::
* Definitions for atensor::


File: maxima.info,  Node: Introduction to atensor,  Next: Definitions for atensor,  Prev: atensor,  Up: atensor

30.1 Introduction to atensor
============================

`atensor' is an algebraic tensor manipulation package. To use `atensor',
type `load(atensor)', followed by a call to the `init_atensor' function.

   The essence of `atensor' is a set of simplification rules for the
noncommutative (dot) product operator ("`.'"). `atensor' recognizes
several algebra types; the corresponding simplification rules are put
into effect when the `init_atensor' function is called.

   The capabilities of `atensor' can be demonstrated by defining the
algebra of quaternions as a Clifford-algebra Cl(0,2) with two basis
vectors. The three quaternionic imaginary units are then the two basis
vectors and their product, i.e.:

         i = v     j = v     k = v  . v
              1         2         1    2

   Although the `atensor' package has a built-in definition for the
quaternion algebra, it is not used in this example, in which we
endeavour to build the quaternion multiplication table as a matrix:


     (%i1) load(atensor);
     (%o1)       /share/tensor/atensor.mac
     (%i2) init_atensor(clifford,0,0,2);
     (%o2)                                done
     (%i3) atensimp(v[1].v[1]);
     (%o3)                                 - 1
     (%i4) atensimp((v[1].v[2]).(v[1].v[2]));
     (%o4)                                 - 1
     (%i5) q:zeromatrix(4,4);
                                     [ 0  0  0  0 ]
                                     [            ]
                                     [ 0  0  0  0 ]
     (%o5)                           [            ]
                                     [ 0  0  0  0 ]
                                     [            ]
                                     [ 0  0  0  0 ]
     (%i6) q[1,1]:1;
     (%o6)                                  1
     (%i7) for i thru adim do q[1,i+1]:q[i+1,1]:v[i];
     (%o7)                                done
     (%i8) q[1,4]:q[4,1]:v[1].v[2];
     (%o8)                               v  . v
                                          1    2
     (%i9) for i from 2 thru 4 do for j from 2 thru 4 do
           q[i,j]:atensimp(q[i,1].q[1,j]);
     (%o9)                                done
     (%i10) q;
                        [    1        v         v      v  . v  ]
                        [              1         2      1    2 ]
                        [                                      ]
                        [   v         - 1     v  . v    - v    ]
                        [    1                 1    2      2   ]
     (%o10)             [                                      ]
                        [   v      - v  . v     - 1      v     ]
                        [    2        1    2              1    ]
                        [                                      ]
                        [ v  . v      v        - v       - 1   ]
                        [  1    2      2          1            ]

   `atensor' recognizes as base vectors indexed symbols, where the
symbol is that stored in `asymbol' and the index runs between 1 and
`adim'.  For indexed symbols, and indexed symbols only, the bilinear
forms `sf', `af', and `av' are evaluated. The evaluation substitutes
the value of `aform[i,j]' in place of `fun(v[i],v[j])' where `v'
represents the value of `asymbol' and `fun' is either `af' or `sf'; or,
it substitutes `v[aform[i,j]]' in place of `av(v[i],v[j])'.

   Needless to say, the functions `sf', `af' and `av' can be redefined.

   When the `atensor' package is loaded, the following flags are set:

     dotscrules:true;
     dotdistrib:true;
     dotexptsimp:false;

   If you wish to experiment with a nonassociative algebra, you may also
consider setting `dotassoc' to `false'. In this case, however,
`atensimp' will not always be able to obtain the desired
simplifications.


File: maxima.info,  Node: Definitions for atensor,  Prev: Introduction to atensor,  Up: atensor

30.2 Definitions for atensor
============================

 -- Function: init_atensor (<alg_type>, <opt_dims>)
 -- Function: init_atensor (<alg_type>)
     Initializes the `atensor' package with the specified algebra type.
     <alg_type> can be one of the following:

     `universal': The universal algebra has no commutation rules.

     `grassmann': The Grassman algebra is defined by the commutation
     relation `u.v+v.u=0'.

     `clifford': The Clifford algebra is defined by the commutation
     relation `u.v+v.u=-2*sf(u,v)' where `sf' is a symmetric
     scalar-valued function. For this algebra, <opt_dims> can be up to
     three nonnegative integers, representing the number of positive,
     degenerate, and negative dimensions of the algebra, respectively.
     If any <opt_dims> values are supplied, `atensor' will configure the
     values of `adim' and `aform' appropriately. Otherwise, `adim' will
     default to 0 and `aform' will not be defined.

     `symmetric': The symmetric algebra is defined by the commutation
     relation `u.v-v.u=0'.

     `symplectic': The symplectic algebra is defined by the commutation
     relation `u.v-v.u=2*af(u,v)' where `af' is an antisymmetric
     scalar-valued function. For the symplectic algebra, <opt_dims> can
     be up to two nonnegative integers, representing the nondegenerate
     and degenerate dimensions, respectively. If any <opt_dims> values
     are supplied, `atensor' will configure the values of `adim' and
     `aform' appropriately. Otherwise, `adim' will default to 0 and
     `aform' will not be defined.

     `lie_envelop': The algebra of the Lie envelope is defined by the
     commutation relation `u.v-v.u=2*av(u,v)' where `av' is an
     antisymmetric function.

     The `init_atensor' function also recognizes several predefined
     algebra types:

     `complex' implements the algebra of complex numbers as the
     Clifford algebra Cl(0,1). The call `init_atensor(complex)' is
     equivalent to `init_atensor(clifford,0,0,1)'.

     `quaternion' implements the algebra of quaternions. The call
     `init_atensor(quaternion)' is equivalent to
     `init_atensor(clifford,0,0,2)'.

     `pauli' implements the algebra of Pauli-spinors as the
     Clifford-algebra Cl(3,0). A call to `init_atensor(pauli)' is
     equivalent to `init_atensor(clifford,3)'.

     `dirac' implements the algebra of Dirac-spinors as the
     Clifford-algebra Cl(3,1). A call to `init_atensor(dirac)' is
     equivalent to `init_atensor(clifford,3,0,1)'.


 -- Function: atensimp (<expr>)
     Simplifies an algebraic tensor expression <expr> according to the
     rules configured by a call to `init_atensor'. Simplification
     includes recursive application of commutation relations and
     resolving calls to `sf', `af', and `av' where applicable. A
     safeguard is used to ensure that the function always terminates,
     even for complex expressions.


 -- Function: alg_type
     The algebra type. Valid values are `universal', `grassmann',
     `clifford', `symmetric', `symplectic' and `lie_envelop'.


 -- Variable: adim
     Default value: 0

     The dimensionality of the algebra. `atensor' uses the value of
     `adim' to determine if an indexed object is a valid base vector.
     See `abasep'.


 -- Variable: aform
     Default value: `ident(3)'

     Default values for the bilinear forms `sf', `af', and `av'. The
     default is the identity matrix `ident(3)'.


 -- Variable: asymbol
     Default value: `v'

     The symbol for base vectors..


 -- Function: sf (<u>, <v>)
     A symmetric scalar function that is used in commutation relations.
     The default implementation checks if both arguments are base
     vectors using `abasep' and if that is the case, substitutes the
     corresponding value from the matrix `aform'.


 -- Function: af (<u>, <v>)
     An antisymmetric scalar function that is used in commutation
     relations.  The default implementation checks if both arguments
     are base vectors using `abasep' and if that is the case,
     substitutes the corresponding value from the matrix `aform'.


 -- Function: av (<u>, <v>)
     An antisymmetric function that is used in commutation relations.
     The default implementation checks if both arguments are base
     vectors using `abasep' and if that is the case, substitutes the
     corresponding value from the matrix `aform'.

     For instance:

          (%i1) load(atensor);
          (%o1)       /share/tensor/atensor.mac
          (%i2) adim:3;
          (%o2)                                  3
          (%i3) aform:matrix([0,3,-2],[-3,0,1],[2,-1,0]);
                                         [  0    3   - 2 ]
                                         [               ]
          (%o3)                          [ - 3   0    1  ]
                                         [               ]
                                         [  2   - 1   0  ]
          (%i4) asymbol:x;
          (%o4)                                  x
          (%i5) av(x[1],x[2]);
          (%o5)                                 x
                                                 3


 -- Function: abasep (<v>)
     Checks if its argument is an `atensor' base vector. That is, if it
     is an indexed symbol, with the symbol being the same as the value
     of `asymbol', and the index having a numeric value between 1 and
     `adim'.



File: maxima.info,  Node: Series,  Next: Number Theory,  Prev: atensor,  Up: Top

31 Series
*********

* Menu:

* Introduction to Series::
* Definitions for Series::


File: maxima.info,  Node: Introduction to Series,  Next: Definitions for Series,  Prev: Series,  Up: Series

31.1 Introduction to Series
===========================

Maxima contains functions `taylor' and `powerseries' for finding the
series of differentiable functions.   It also has tools such as `nusum'
capable of finding the closed form of some series.   Operations such as
addition and multiplication work as usual on series. This section
presents the global variables which control the expansion.


File: maxima.info,  Node: Definitions for Series,  Prev: Introduction to Series,  Up: Series

31.2 Definitions for Series
===========================

 -- Option variable: cauchysum
     Default value: `false'

     When multiplying together sums with `inf' as their upper limit, if
     `sumexpand' is `true' and `cauchysum' is `true' then the Cauchy
     product will be used rather than the usual product.  In the Cauchy
     product the index of the inner summation is a function of the
     index of the outer one rather than varying independently.

     Example:

          (%i1) sumexpand: false$
          (%i2) cauchysum: false$
          (%i3) s: sum (f(i), i, 0, inf) * sum (g(j), j, 0, inf);
                                inf         inf
                                ====        ====
                                \           \
          (%o3)                ( >    f(i))  >    g(j)
                                /           /
                                ====        ====
                                i = 0       j = 0
          (%i4) sumexpand: true$
          (%i5) cauchysum: true$
          (%i6) ''s;
                           inf     i1
                           ====   ====
                           \      \
          (%o6)             >      >     g(i1 - i2) f(i2)
                           /      /
                           ====   ====
                           i1 = 0 i2 = 0


 -- Function: deftaylor (<f_1>(<x_1>), <expr_1>, ..., <f_n>(<x_n>),
          <expr_n>)
     For each function <f_i> of one variable <x_i>, `deftaylor' defines
     <expr_i> as the Taylor series about zero.  <expr_i> is typically a
     polynomial in <x_i> or a summation; more general expressions are
     accepted by `deftaylor' without complaint.

     `powerseries (<f_i>(<x_i>), <x_i>, 0)' returns the series defined
     by `deftaylor'.

     `deftaylor' returns a list of the functions <f_1>, ..., <f_n>.
     `deftaylor' evaluates its arguments.

     Example:

          (%i1) deftaylor (f(x), x^2 + sum(x^i/(2^i*i!^2), i, 4, inf));
          (%o1)                          [f]
          (%i2) powerseries (f(x), x, 0);
                                inf
                                ====      i1
                                \        x         2
          (%o2)                  >     -------- + x
                                /       i1    2
                                ====   2   i1!
                                i1 = 4
          (%i3) taylor (exp (sqrt (f(x))), x, 0, 4);
                                2         3          4
                               x    3073 x    12817 x
          (%o3)/T/     1 + x + -- + ------- + -------- + . . .
                               2     18432     307200


 -- Option variable: maxtayorder
     Default value: `true'

     When `maxtayorder' is `true', then during algebraic manipulation
     of (truncated) Taylor series, `taylor' tries to retain as many
     terms as are known to be correct.


 -- Function: niceindices (<expr>)
     Renames the indices of sums and products in <expr>.  `niceindices'
     attempts to rename each index to the value of `niceindicespref[1]',
     unless that name appears in the summand or multiplicand, in which
     case `niceindices' tries the succeeding elements of
     `niceindicespref' in turn, until an unused variable is found.  If
     the entire list is exhausted, additional indices are constructed
     by appending integers to the value of `niceindicespref[1]', e.g.,
     `i0', `i1', `i2', ....

     `niceindices' returns an expression.  `niceindices' evaluates its
     argument.

     Example:

          (%i1) niceindicespref;
          (%o1)                  [i, j, k, l, m, n]
          (%i2) product (sum (f (foo + i*j*bar), foo, 1, inf), bar, 1, inf);
                           inf    inf
                          /===\   ====
                           ! !    \
          (%o2)            ! !     >      f(bar i j + foo)
                           ! !    /
                          bar = 1 ====
                                  foo = 1
          (%i3) niceindices (%);
                               inf  inf
                              /===\ ====
                               ! !  \
          (%o3)                ! !   >    f(i j l + k)
                               ! !  /
                              l = 1 ====
                                    k = 1


 -- Option variable: niceindicespref
     Default value: `[i, j, k, l, m, n]'

     `niceindicespref' is the list from which `niceindices' takes the
     names of indices for sums and products.

     The elements of `niceindicespref' are typically names of variables,
     although that is not enforced by `niceindices'.

     Example:

          (%i1) niceindicespref: [p, q, r, s, t, u]$
          (%i2) product (sum (f (foo + i*j*bar), foo, 1, inf), bar, 1, inf);
                           inf    inf
                          /===\   ====
                           ! !    \
          (%o2)            ! !     >      f(bar i j + foo)
                           ! !    /
                          bar = 1 ====
                                  foo = 1
          (%i3) niceindices (%);
                               inf  inf
                              /===\ ====
                               ! !  \
          (%o3)                ! !   >    f(i j q + p)
                               ! !  /
                              q = 1 ====
                                    p = 1


 -- Function: nusum (<expr>, <x>, <i_0>, <i_1>)
     Carries out indefinite hypergeometric summation of <expr> with
     respect to <x> using a decision procedure due to R.W. Gosper.
     <expr> and the result must be expressible as products of integer
     powers, factorials, binomials, and rational functions.

     The terms "definite" and "indefinite summation" are used
     analogously to "definite" and "indefinite integration".  To sum
     indefinitely means to give a symbolic result for the sum over
     intervals of variable length, not just e.g. 0 to inf.  Thus, since
     there is no formula for the general partial sum of the binomial
     series, `nusum' can't do it.

     `nusum' and `unsum' know a little about sums and differences of
     finite products.  See also `unsum'.

     Examples:

          (%i1) nusum (n*n!, n, 0, n);

          Dependent equations eliminated:  (1)
          (%o1)                     (n + 1)! - 1
          (%i2) nusum (n^4*4^n/binomial(2*n,n), n, 0, n);
                               4        3       2              n
                2 (n + 1) (63 n  + 112 n  + 18 n  - 22 n + 3) 4      2
          (%o2) ------------------------------------------------ - ------
                              693 binomial(2 n, n)                 3 11 7
          (%i3) unsum (%, n);
                                        4  n
                                       n  4
          (%o3)                   ----------------
                                  binomial(2 n, n)
          (%i4) unsum (prod (i^2, i, 1, n), n);
                              n - 1
                              /===\
                               ! !   2
          (%o4)              ( ! !  i ) (n - 1) (n + 1)
                               ! !
                              i = 1
          (%i5) nusum (%, n, 1, n);

          Dependent equations eliminated:  (2 3)
                                      n
                                    /===\
                                     ! !   2
          (%o5)                      ! !  i  - 1
                                     ! !
                                    i = 1


 -- Function: pade (<taylor_series>, <numer_deg_bound>,
          <denom_deg_bound>)
     Returns a list of all rational functions which have the given
     Taylor series expansion where the sum of the degrees of the
     numerator and the denominator is less than or equal to the
     truncation level of the power series, i.e.  are "best"
     approximants, and which additionally satisfy the specified degree
     bounds.

     <taylor_series> is a univariate Taylor series.  <numer_deg_bound>
     and <denom_deg_bound> are positive integers specifying degree
     bounds on the numerator and denominator.

     <taylor_series> can also be a Laurent series, and the degree
     bounds can be `inf' which causes all rational functions whose total
     degree is less than or equal to the length of the power series to
     be returned.  Total degree is defined as `<numer_deg_bound> +
     <denom_deg_bound>'.  Length of a power series is defined as
     `"truncation level" + 1 - min(0, "order of series")'.

          (%i1) taylor (1 + x + x^2 + x^3, x, 0, 3);
                                        2    3
          (%o1)/T/             1 + x + x  + x  + . . .
          (%i2) pade (%, 1, 1);
                                           1
          (%o2)                       [- -----]
                                         x - 1
          (%i3) t: taylor(-(83787*x^10 - 45552*x^9 - 187296*x^8
                             + 387072*x^7 + 86016*x^6 - 1507328*x^5
                             + 1966080*x^4 + 4194304*x^3 - 25165824*x^2
                             + 67108864*x - 134217728)
                 /134217728, x, 0, 10);
                              2    3       4       5       6        7
                       x   3 x    x    15 x    23 x    21 x    189 x
          (%o3)/T/ 1 - - + ---- - -- - ----- + ----- - ----- - ------
                       2    16    32   1024    2048    32768   65536

                                            8         9          10
                                      5853 x    2847 x    83787 x
                                    + ------- + ------- - --------- + . . .
                                      4194304   8388608   134217728
          (%i4) pade (t, 4, 4);
          (%o4)                          []

     There is no rational function of degree 4 numerator/denominator,
     with this power series expansion.  You must in general have degree
     of the numerator and degree of the denominator adding up to at
     least the degree of the power series, in order to have enough
     unknown coefficients to solve.

          (%i5) pade (t, 5, 5);
                               5                4                 3
          (%o5) [- (520256329 x  - 96719020632 x  - 489651410240 x

                            2
           - 1619100813312 x  - 2176885157888 x - 2386516803584)

                         5                 4                  3
          /(47041365435 x  + 381702613848 x  + 1360678489152 x

                            2
           + 2856700692480 x  + 3370143559680 x + 2386516803584)]


 -- Option variable: powerdisp
     Default value: `false'

     When `powerdisp' is `true', a sum is displayed with its terms in
     order of increasing power.  Thus a polynomial is displayed as a
     truncated power series, with the constant term first and the
     highest power last.

     By default, terms of a sum are displayed in order of decreasing
     power.


 -- Function: powerseries (<expr>, <x>, <a>)
     Returns the general form of the power series expansion for <expr>
     in the variable <x> about the point <a> (which may be `inf' for
     infinity).

     If `powerseries' is unable to expand <expr>, `taylor' may give the
     first several terms of the series.

     When `verbose' is `true', `powerseries' prints progress messages.

          (%i1) verbose: true$
          (%i2) powerseries (log(sin(x)/x), x, 0);
          can't expand
                                           log(sin(x))
          so we'll try again after applying the rule:
                                                  d
                                                / -- (sin(x))
                                                [ dx
                                  log(sin(x)) = i ----------- dx
                                                ]   sin(x)
                                                /
          in the first simplification we have returned:
                                       /
                                       [
                                       i cot(x) dx - log(x)
                                       ]
                                       /
                              inf
                              ====        i1  2 i1             2 i1
                              \      (- 1)   2     bern(2 i1) x
                               >     ------------------------------
                              /                i1 (2 i1)!
                              ====
                              i1 = 1
          (%o2)                -------------------------------------
                                                2


 -- Option variable: psexpand
     Default value: `false'

     When `psexpand' is `true', an extended rational function
     expression is displayed fully expanded.  The switch `ratexpand'
     has the same effect.

     When `psexpand' is `false', a multivariate expression is displayed
     just as in the rational function package.

     When `psexpand' is  `multi', then terms with the same total degree
     in the variables are grouped together.


 -- Function: revert (<expr>, <x>)
 -- Function: revert2 (<expr>, <x>, <n>)
     These functions return the reversion of <expr>, a Taylor series
     about zero in the variable <x>.  `revert' returns a polynomial of
     degree equal to the highest power in <expr>.  `revert2' returns a
     polynomial of degree <n>, which may be greater than, equal to, or
     less than the degree of <expr>.

     `load ("revert")' loads these functions.

     Examples:

          (%i1) load ("revert")$
          (%i2) t: taylor (exp(x) - 1, x, 0, 6);
                             2    3    4    5     6
                            x    x    x    x     x
          (%o2)/T/      x + -- + -- + -- + --- + --- + . . .
                            2    6    24   120   720
          (%i3) revert (t, x);
                         6       5       4       3       2
                     10 x  - 12 x  + 15 x  - 20 x  + 30 x  - 60 x
          (%o3)/R/ - --------------------------------------------
                                          60
          (%i4) ratexpand (%);
                               6    5    4    3    2
                              x    x    x    x    x
          (%o4)             - -- + -- - -- + -- - -- + x
                              6    5    4    3    2
          (%i5) taylor (log(x+1), x, 0, 6);
                              2    3    4    5    6
                             x    x    x    x    x
          (%o5)/T/       x - -- + -- - -- + -- - -- + . . .
                             2    3    4    5    6
          (%i6) ratsimp (revert (t, x) - taylor (log(x+1), x, 0, 6));
          (%o6)                           0
          (%i7) revert2 (t, x, 4);
                                    4    3    2
                                   x    x    x
          (%o7)                  - -- + -- - -- + x
                                   4    3    2


 -- Function: taylor (<expr>, <x>, <a>, <n>)
 -- Function: taylor (<expr>, [<x_1>, <x_2>, ...], <a>, <n>)
 -- Function: taylor (<expr>, [<x>, <a>, <n>, 'asymp])
 -- Function: taylor (<expr>, [<x_1>, <x_2>, ...], [<a_1>, <a_2>, ...],
          [<n_1>, <n_2>, ...])
     `taylor (<expr>, <x>, <a>, <n>)' expands the expression <expr> in
     a truncated Taylor or Laurent series in the variable <x> around
     the point <a>, containing terms through `(<x> - <a>)^<n>'.

     If <expr> is of the form `<f>(<x>)/<g>(<x>)' and `<g>(<x>)' has no
     terms up to degree <n> then `taylor' attempts to expand `<g>(<x>)'
     up to degree `2 <n>'.  If there are still no nonzero terms,
     `taylor' doubles the degree of the expansion of `<g>(<x>)' so long
     as the degree of the expansion is less than or equal to `<n>
     2^taylordepth'.

     `taylor (<expr>, [<x_1>, <x_2>, ...], <a>, <n>)' returns a
     truncated power series of degree <n> in all variables <x_1>,
     <x_2>, ...  about the point `(<a>, <a>, ...)'.

     `taylor (<expr>, [<x_1>, <a_1>, <n_1>], [<x_2>, <a_2>, <n_2>],
     ...)'  returns a truncated power series in the variables <x_1>,
     <x_2>, ...  about the point `(<a_1>, <a_2>, ...)', truncated at
     <n_1>, <n_2>, ....

     `taylor (<expr>, [<x_1>, <x_2>, ...], [<a_1>, <a_2>, ...], [<n_1>,
     <n_2>, ...])' returns a truncated power series in the variables
     <x_1>, <x_2>, ...  about the point `(<a_1>, <a_2>, ...)',
     truncated at <n_1>, <n_2>, ....

     `taylor (<expr>, [<x>, <a>, <n>, 'asymp])' returns an expansion of
     <expr> in negative powers of `<x> - <a>'.  The highest order term
     is `(<x> - <a>)^<-n>'.

     When `maxtayorder' is `true', then during algebraic manipulation
     of (truncated) Taylor series, `taylor' tries to retain as many
     terms as are known to be correct.

     When `psexpand' is `true', an extended rational function
     expression is displayed fully expanded.  The switch `ratexpand'
     has the same effect.  When `psexpand' is `false', a multivariate
     expression is displayed just as in the rational function package.
     When `psexpand' is  `multi', then terms with the same total degree
     in the variables are grouped together.

     See also the `taylor_logexpand' switch for controlling expansion.

     Examples:

          (%i1) taylor (sqrt (sin(x) + a*x + 1), x, 0, 3);
                                     2             2
                       (a + 1) x   (a  + 2 a + 1) x
          (%o1)/T/ 1 + --------- - -----------------
                           2               8

                                             3      2             3
                                         (3 a  + 9 a  + 9 a - 1) x
                                       + -------------------------- + . . .
                                                     48
          (%i2) %^2;
                                              3
                                             x
          (%o2)/T/           1 + (a + 1) x - -- + . . .
                                             6
          (%i3) taylor (sqrt (x + 1), x, 0, 5);
                                 2    3      4      5
                            x   x    x    5 x    7 x
          (%o3)/T/      1 + - - -- + -- - ---- + ---- + . . .
                            2   8    16   128    256
          (%i4) %^2;
          (%o4)/T/                  1 + x + . . .
          (%i5) product ((1 + x^i)^2.5, i, 1, inf)/(1 + x^2);
                                   inf
                                  /===\
                                   ! !    i     2.5
                                   ! !  (x  + 1)
                                   ! !
                                  i = 1
          (%o5)                   -----------------
                                        2
                                       x  + 1
          (%i6) ev (taylor(%, x,  0, 3), keepfloat);
                                         2           3
          (%o6)/T/    1 + 2.5 x + 3.375 x  + 6.5625 x  + . . .
          (%i7) taylor (1/log (x + 1), x, 0, 3);
                                         2       3
                           1   1   x    x    19 x
          (%o7)/T/         - + - - -- + -- - ----- + . . .
                           x   2   12   24    720
          (%i8) taylor (cos(x) - sec(x), x, 0, 5);
                                          4
                                     2   x
          (%o8)/T/                - x  - -- + . . .
                                         6
          (%i9) taylor ((cos(x) - sec(x))^3, x, 0, 5);
          (%o9)/T/                    0 + . . .
          (%i10) taylor (1/(cos(x) - sec(x))^3, x, 0, 5);
                                                         2          4
                      1     1       11      347    6767 x    15377 x
          (%o10)/T/ - -- + ---- + ------ - ----- - ------- - --------
                       6      4        2   15120   604800    7983360
                      x    2 x    120 x

                                                                    + . . .
          (%i11) taylor (sqrt (1 - k^2*sin(x)^2), x, 0, 6);
                         2  2       4      2   4
                        k  x    (3 k  - 4 k ) x
          (%o11)/T/ 1 - ----- - ----------------
                          2            24

                                              6       4       2   6
                                         (45 k  - 60 k  + 16 k ) x
                                       - -------------------------- + . . .
                                                    720
          (%i12) taylor ((x + 1)^n, x, 0, 4);
                                2       2     3      2         3
                              (n  - n) x    (n  - 3 n  + 2 n) x
          (%o12)/T/ 1 + n x + ----------- + --------------------
                                   2                 6

                                         4      3       2         4
                                       (n  - 6 n  + 11 n  - 6 n) x
                                     + ---------------------------- + . . .
                                                    24
          (%i13) taylor (sin (y + x), x, 0, 3, y, 0, 3);
                         3                 2
                        y                 y
          (%o13)/T/ y - -- + . . . + (1 - -- + . . .) x
                        6                 2

                              3                       2
                         y   y            2      1   y            3
                    + (- - + -- + . . .) x  + (- - + -- + . . .) x  + . . .
                         2   12                  6   12
          (%i14) taylor (sin (y + x), [x, y], 0, 3);
                               3        2      2      3
                              x  + 3 y x  + 3 y  x + y
          (%o14)/T/   y + x - ------------------------- + . . .
                                          6
          (%i15) taylor (1/sin (y + x), x, 0, 3, y, 0, 3);
                    1   y              1    1               1            2
          (%o15)/T/ - + - + . . . + (- -- + - + . . .) x + (-- + . . .) x
                    y   6               2   6                3
                                       y                    y

                                                     1            3
                                                + (- -- + . . .) x  + . . .
                                                      4
                                                     y
          (%i16) taylor (1/sin (y + x), [x, y], 0, 3);
                                       3         2       2        3
                      1     x + y   7 x  + 21 y x  + 21 y  x + 7 y
          (%o16)/T/ ----- + ----- + ------------------------------- + . . .
                    x + y     6                   360


 -- Option variable: taylordepth
     Default value: 3

     If there are still no nonzero terms, `taylor' doubles the degree
     of the expansion of `<g>(<x>)' so long as the degree of the
     expansion is less than or equal to `<n> 2^taylordepth'.


 -- Function: taylorinfo (<expr>)
     Returns information about the Taylor series <expr>.  The return
     value is a list of lists.  Each list comprises the name of a
     variable, the point of expansion, and the degree of the expansion.

     `taylorinfo' returns `false' if <expr> is not a Taylor series.

     Example:

          (%i1) taylor ((1 - y^2)/(1 - x), x, 0, 3, [y, a, inf]);
                            2                       2
          (%o1)/T/ - (y - a)  - 2 a (y - a) + (1 - a )

                   2                        2
           + (1 - a  - 2 a (y - a) - (y - a) ) x

                   2                        2   2
           + (1 - a  - 2 a (y - a) - (y - a) ) x

                   2                        2   3
           + (1 - a  - 2 a (y - a) - (y - a) ) x  + . . .
          (%i2) taylorinfo(%);
          (%o2)               [[y, a, inf], [x, 0, 3]]


 -- Function: taylorp (<expr>)
     Returns `true' if <expr> is a Taylor series, and `false' otherwise.


 -- Option variable: taylor_logexpand
     Default value: `true'

     `taylor_logexpand' controls expansions of logarithms in `taylor'
     series.

     When `taylor_logexpand' is `true', all logarithms are expanded
     fully so that zero-recognition problems involving logarithmic
     identities do not disturb the expansion process.  However, this
     scheme is not always mathematically correct since it ignores
     branch information.

     When `taylor_logexpand' is set to `false', then the only expansion
     of logarithms that occur is that necessary to obtain a formal
     power series.


 -- Option variable: taylor_order_coefficients
     Default value: `true'

     `taylor_order_coefficients' controls the ordering of coefficients
     in a Taylor series.

     When `taylor_order_coefficients' is `true', coefficients of taylor
     series are ordered canonically.


 -- Function: taylor_simplifier (<expr>)
     Simplifies coefficients of the power series <expr>.  `taylor'
     calls this function.


 -- Option variable: taylor_truncate_polynomials
     Default value: `true'

     When `taylor_truncate_polynomials' is `true', polynomials are
     truncated based upon the input truncation levels.

     Otherwise, polynomials input to `taylor' are considered to have
     infinite precison.


 -- Function: taytorat (<expr>)
     Converts <expr> from `taylor' form to canonical rational
     expression (CRE) form.  The effect is the same as `rat (ratdisrep
     (<expr>))', but faster.


 -- Function: trunc (<expr>)
     Annotates the internal representation of the general expression
     <expr> so that it is displayed as if its sums were truncated
     Taylor series.  <expr> is not otherwise modified.

     Example:

          (%i1) expr: x^2 + x + 1;
                                      2
          (%o1)                      x  + x + 1
          (%i2) trunc (expr);
                                          2
          (%o2)                  1 + x + x  + . . .
          (%i3) is (expr = trunc (expr));
          (%o3)                         true


 -- Function: unsum (<f>, <n>)
     Returns the first backward difference `<f>(<n>) - <f>(<n> - 1)'.
     Thus `unsum' in a sense is the inverse of `sum'.

     See also `nusum'.

     Examples:

          (%i1) g(p) := p*4^n/binomial(2*n,n);
                                               n
                                            p 4
          (%o1)               g(p) := ----------------
                                      binomial(2 n, n)
          (%i2) g(n^4);
                                        4  n
                                       n  4
          (%o2)                   ----------------
                                  binomial(2 n, n)
          (%i3) nusum (%, n, 0, n);
                               4        3       2              n
                2 (n + 1) (63 n  + 112 n  + 18 n  - 22 n + 3) 4      2
          (%o3) ------------------------------------------------ - ------
                              693 binomial(2 n, n)                 3 11 7
          (%i4) unsum (%, n);
                                        4  n
                                       n  4
          (%o4)                   ----------------
                                  binomial(2 n, n)


 -- Option variable: verbose
     Default value: `false'

     When `verbose' is `true', `powerseries' prints progress messages.



File: maxima.info,  Node: Number Theory,  Next: Symmetries,  Prev: Series,  Up: Top

32 Number Theory
****************

* Menu:

* Definitions for Number Theory::


File: maxima.info,  Node: Definitions for Number Theory,  Prev: Number Theory,  Up: Number Theory

32.1 Definitions for Number Theory
==================================

 -- Function: bern (<n>)
     Returns the <n>'th Bernoulli number for integer <n>.  Bernoulli
     numbers equal to zero are suppressed if `zerobern' is `false'.

     See also `burn'.

          (%i1) zerobern: true$
          (%i2) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
                            1  1       1      1        1
          (%o2)       [1, - -, -, 0, - --, 0, --, 0, - --]
                            2  6       30     42       30
          (%i3) zerobern: false$
          (%i4) map (bern, [0, 1, 2, 3, 4, 5, 6, 7, 8]);
                      1  1    1   5     691   7    3617  43867
          (%o4) [1, - -, -, - --, --, - ----, -, - ----, -----]
                      2  6    30  66    2730  6    510    798


 -- Function: bernpoly (<x>, <n>)
     Returns the <n>'th Bernoulli polynomial in the variable <x>.


 -- Function: bfzeta (<s>, <n>)
     Returns the Riemann zeta function for the argument <s>.  The
     return value is a big float (bfloat); <n> is the number of digits
     in the return value.

     `load ("bffac")' loads this function.


 -- Function: bfhzeta (<s>, <h>, <n>)
     Returns the Hurwitz zeta function for the arguments <s> and <h>.
     The return value is a big float (bfloat); <n> is the number of
     digits in the return value.

     The Hurwitz zeta function is defined as

          sum ((k+h)^-s, k, 0, inf)

     `load ("bffac")' loads this function.


 -- Function: binomial (<x>, <y>)
     The binomial coefficient `<x>!/(<y>! (<x> - <y>)!)'.  If <x> and
     <y> are integers, then the numerical value of the binomial
     coefficient is computed.  If <y>, or <x - y>, is an integer, the
     binomial coefficient is expressed as a polynomial.

     Examples:

          (%i1) binomial (11, 7);
          (%o1)                          330
          (%i2) 11! / 7! / (11 - 7)!;
          (%o2)                          330
          (%i3) binomial (x, 7);
                  (x - 6) (x - 5) (x - 4) (x - 3) (x - 2) (x - 1) x
          (%o3)   -------------------------------------------------
                                        5040
          (%i4) binomial (x + 7, x);
                (x + 1) (x + 2) (x + 3) (x + 4) (x + 5) (x + 6) (x + 7)
          (%o4) -------------------------------------------------------
                                         5040
          (%i5) binomial (11, y);
          (%o5)                    binomial(11, y)


 -- Function: burn (<n>)
     Returns the <n>'th Bernoulli number for integer <n>.  `burn' may
     be more efficient than `bern' for large, isolated <n> (perhaps <n>
     greater than 105 or so), as `bern' computes all the Bernoulli
     numbers up to index <n> before returning.

     `burn' exploits the observation that (rational) Bernoulli numbers
     can be approximated by (transcendental) zetas with tolerable
     efficiency.

     `load ("bffac")' loads this function.


 -- Function: cf (<expr>)
     Converts <expr> into a continued fraction.  <expr> is an expression
     comprising continued fractions and square roots of integers.
     Operands in the expression may be combined with arithmetic
     operators.  Aside from continued fractions and square roots,
     factors in the expression must be integer or rational numbers.
     Maxima does not know about operations on continued fractions
     outside of `cf'.

     `cf' evaluates its arguments after binding `listarith' to `false'.
     `cf' returns a continued fraction, represented as a list.

     A continued fraction `a + 1/(b + 1/(c + ...))' is represented by
     the list `[a, b, c, ...]'.  The list elements `a', `b', `c', ...
     must evaluate to integers.  <expr> may also contain `sqrt (n)'
     where `n' is an integer.  In this case `cf' will give as many
     terms of the continued fraction as the value of the variable
     `cflength' times the period.

     A continued fraction can be evaluated to a number by evaluating
     the arithmetic representation returned by `cfdisrep'.  See also
     `cfexpand' for another way to evaluate a continued fraction.

     See also `cfdisrep', `cfexpand', and `cflength'.

     Examples:

        * <expr> is an expression comprising continued fractions and
          square roots of integers.

               (%i1) cf ([5, 3, 1]*[11, 9, 7] + [3, 7]/[4, 3, 2]);
               (%o1)               [59, 17, 2, 1, 1, 1, 27]
               (%i2) cf ((3/17)*[1, -2, 5]/sqrt(11) + (8/13));
               (%o2)        [0, 1, 1, 1, 3, 2, 1, 4, 1, 9, 1, 9, 2]

        * `cflength' controls how many periods of the continued fraction
          are computed for algebraic, irrational numbers.

               (%i1) cflength: 1$
               (%i2) cf ((1 + sqrt(5))/2);
               (%o2)                    [1, 1, 1, 1, 2]
               (%i3) cflength: 2$
               (%i4) cf ((1 + sqrt(5))/2);
               (%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
               (%i5) cflength: 3$
               (%i6) cf ((1 + sqrt(5))/2);
               (%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]

        * A continued fraction can be evaluated by evaluating the
          arithmetic representation returned by `cfdisrep'.

               (%i1) cflength: 3$
               (%i2) cfdisrep (cf (sqrt (3)))$
               (%i3) ev (%, numer);
               (%o3)                   1.731707317073171

        * Maxima does not know about operations on continued fractions
          outside of `cf'.

               (%i1) cf ([1,1,1,1,1,2] * 3);
               (%o1)                     [4, 1, 5, 2]
               (%i2) cf ([1,1,1,1,1,2]) * 3;
               (%o2)                  [3, 3, 3, 3, 3, 6]


 -- Function: cfdisrep (<list>)
     Constructs and returns an ordinary arithmetic expression of the
     form `a + 1/(b + 1/(c + ...))' from the list representation of a
     continued fraction `[a, b, c, ...]'.

          (%i1) cf ([1, 2, -3] + [1, -2, 1]);
          (%o1)                     [1, 1, 1, 2]
          (%i2) cfdisrep (%);
                                            1
          (%o2)                     1 + ---------
                                              1
                                        1 + -----
                                                1
                                            1 + -
                                                2


 -- Function: cfexpand (<x>)
     Returns a matrix of the numerators and denominators of the last
     (column 1) and next-to-last (column 2) convergents of the
     continued fraction <x>.

          (%i1) cf (rat (ev (%pi, numer)));

          `rat' replaced 3.141592653589793 by 103993//33102 = 3.141592653011902
          (%o1)                  [3, 7, 15, 1, 292]
          (%i2) cfexpand (%);
                                   [ 103993  355 ]
          (%o2)                    [             ]
                                   [ 33102   113 ]
          (%i3) %[1,1]/%[2,1], numer;
          (%o3)                   3.141592653011902


 -- Option variable: cflength
     Default value: 1

     `cflength' controls the number of terms of the continued fraction
     the function `cf' will give, as the value `cflength' times the
     period.  Thus the default is to give one period.

          (%i1) cflength: 1$
          (%i2) cf ((1 + sqrt(5))/2);
          (%o2)                    [1, 1, 1, 1, 2]
          (%i3) cflength: 2$
          (%i4) cf ((1 + sqrt(5))/2);
          (%o4)               [1, 1, 1, 1, 1, 1, 1, 2]
          (%i5) cflength: 3$
          (%i6) cf ((1 + sqrt(5))/2);
          (%o6)           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2]


 -- Function: divsum (<n>, <k>)
 -- Function: divsum (<n>)
     `divsum (<n>, <k>)' returns the sum of the divisors of <n> raised
     to the <k>'th power.

     `divsum (<n>)' returns the sum of the divisors of <n>.

          (%i1) divsum (12);
          (%o1)                          28
          (%i2) 1 + 2 + 3 + 4 + 6 + 12;
          (%o2)                          28
          (%i3) divsum (12, 2);
          (%o3)                          210
          (%i4) 1^2 + 2^2 + 3^2 + 4^2 + 6^2 + 12^2;
          (%o4)                          210


 -- Function: euler (<n>)
     Returns the <n>'th Euler number for nonnegative integer <n>.

     For the Euler-Mascheroni constant, see `%gamma'.

          (%i1) map (euler, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
          (%o1)    [1, 0, - 1, 0, 5, 0, - 61, 0, 1385, 0, - 50521]


 -- Constant: %gamma
     The Euler-Mascheroni constant, 0.5772156649015329 ....


 -- Function: factorial (<x>)
     Represents the factorial function. Maxima treats `factorial (<x>)'
     the same as `<x>!'.  See `!'.


 -- Function: fib (<n>)
     Returns the <n>'th Fibonacci number.  `fib(0)' equal to 0 and
     `fib(1)' equal to 1, and `fib (-<n>)' equal to `(-1)^(<n> + 1) *
     fib(<n>)'.

     After calling `fib', `prevfib' is equal to `fib (<x> - 1)', the
     Fibonacci number preceding the last one computed.

          (%i1) map (fib, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]);
          (%o1)         [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55]


 -- Function: fibtophi (<expr>)
     Expresses Fibonacci numbers in terms of the constant `%phi', which
     is `(1 + sqrt(5))/2', approximately 1.61803399.

     By default, Maxima does not know about `%phi'.  After executing
     `tellrat (%phi^2 - %phi - 1)' and `algebraic: true', `ratsimp' can
     simplify some expressions containing `%phi'.

          (%i1) fibtophi (fib (n));
                                     n             n
                                 %phi  - (1 - %phi)
          (%o1)                  -------------------
                                     2 %phi - 1
          (%i2) fib (n-1) + fib (n) - fib (n+1);
          (%o2)          - fib(n + 1) + fib(n) + fib(n - 1)
          (%i3) ratsimp (fibtophi (%));
          (%o3)                           0


 -- Function: ifactors (<n>)
     For a positive integer <n> returns the factorization of <n>. If
     `n=p1^e1..pk^nk' is the decomposition of <n> into prime factors,
     ifactors returns `[[p1, e1], ... , [pk, ek]]'.

     Factorization methods used are trial divisions by primes up to
     9973, Pollard's rho method and elliptic curve method.

          (%i1) ifactors(51575319651600);
          (%o1)     [[2, 4], [3, 2], [5, 2], [1583, 1], [9050207, 1]]
          (%i2) apply("*", map(lambda([u], u[1]^u[2]), %));
          (%o2)                        51575319651600


 -- Function: inrt (<x>, <n>)
     Returns the integer <n>'th root of the absolute value of <x>.

          (%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
          (%i2) map (lambda ([a], inrt (10^a, 3)), l);
          (%o2) [2, 4, 10, 21, 46, 100, 215, 464, 1000, 2154, 4641, 10000]


 -- Function: inv_mod (<n>, <m>)
     Computes the inverse of <n> modulo <m>.  `inv_mod (n,m)' returns
     `false', if <n> is a zero divisor modulo <m>.

          (%i1) inv_mod(3, 41);
          (%o1)                           14
          (%i2) ratsimp(3^-1), modulus=41;
          (%o2)                           14
          (%i3) inv_mod(3, 42);
          (%o3)                          false


 -- Function: jacobi (<p>, <q>)
     Returns the Jacobi symbol of <p> and <q>.

          (%i1) l: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
          (%i2) map (lambda ([a], jacobi (a, 9)), l);
          (%o2)         [1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]


 -- Function: lcm (<expr_1>, ..., <expr_n>)
     Returns the least common multiple of its arguments.  The arguments
     may be general expressions as well as integers.

     `load ("functs")' loads this function.


 -- Function: minfactorial (<expr>)
     Examines <expr> for occurrences of two factorials which differ by
     an integer.  `minfactorial' then turns one into a polynomial times
     the other.

          (%i1) n!/(n+2)!;
                                         n!
          (%o1)                       --------
                                      (n + 2)!
          (%i2) minfactorial (%);
                                          1
          (%o2)                    ---------------
                                   (n + 1) (n + 2)


 -- Function: power_mod (<a>, <n>, <m>)
     Uses a modular algorithm to compute `a^n mod m' where <a> and <n>
     are integers and <m> is a positive integer.  If <n> is negative,
     `inv_mod' is used to find the modular inverse.

          (%i1) power_mod(3, 15, 5);
          (%o1)                          2
          (%i2) mod(3^15,5);
          (%o2)                          2
          (%i3) power_mod(2, -1, 5);
          (%o3)                          3
          (%i4) inv_mod(2,5);
          (%o4)                          3


 -- Function: next_prime (<n>)
     Returns the smallest prime bigger than <n>.

          (%i1) next_prime(27);
          (%o1)                       29


 -- Function: partfrac (<expr>, <var>)
     Expands the expression <expr> in partial fractions with respect to
     the main variable <var>.  `partfrac' does a complete partial
     fraction decomposition.  The algorithm employed is based on the
     fact that the denominators of the partial fraction expansion (the
     factors of the original denominator) are relatively prime.  The
     numerators can be written as linear combinations of denominators,
     and the expansion falls out.

          (%i1) 1/(1+x)^2 - 2/(1+x) + 2/(2+x);
                                2       2        1
          (%o1)               ----- - ----- + --------
                              x + 2   x + 1          2
                                              (x + 1)
          (%i2) ratsimp (%);
                                           x
          (%o2)                 - -------------------
                                   3      2
                                  x  + 4 x  + 5 x + 2
          (%i3) partfrac (%, x);
                                2       2        1
          (%o3)               ----- - ----- + --------
                              x + 2   x + 1          2
                                              (x + 1)

 -- Function: primep (<n>)
     Primality test. If `primep (n)' returns `false', <n> is a
     composite number and if it returns `true', <n> is a prime number
     with very high probability.

     For <n> less than 341550071728321 a deterministic version of
     Miller-Rabin's test is used. If `primep (n)' returns `true', then
     <n> is a prime number.

     For <n> bigger than 34155071728321 `primep' uses
     `primep_number_of_tests' Miller-Rabin's pseudo-primality tests and
     one Lucas pseudo-primality test. The probability that <n> will
     pass one Miller-Rabin test is less than 1/4. Using the default
     value 25 for `primep_number_of_tests', the probability of <n>
     beeing composite is much smaller that 10^-15.


 -- Option variable: primep_number_of_tests
     Default value: 25

     Number of Miller-Rabin's tests used in `primep'.

 -- Function: prev_prime (<n>)
     Returns the greatest prime smaller than <n>.

          (%i1) prev_prime(27);
          (%o1)                       23

 -- Function: qunit (<n>)
     Returns the principal unit of the real quadratic number field
     `sqrt (<n>)' where <n> is an integer, i.e., the element whose norm
     is unity.  This amounts to solving Pell's equation `a^2 - <n> b^2
     = 1'.

          (%i1) qunit (17);
          (%o1)                     sqrt(17) + 4
          (%i2) expand (% * (sqrt(17) - 4));
          (%o2)                           1


 -- Function: totient (<n>)
     Returns the number of integers less than or equal to <n> which are
     relatively prime to <n>.


 -- Option variable: zerobern
     Default value: `true'

     When `zerobern' is `false', `bern' excludes the Bernoulli numbers
     which are equal to zero.  See `bern'.


 -- Function: zeta (<n>)
     Returns the Riemann zeta function if <x> is a negative integer, 0,
     1, or a positive even number, and returns a noun form `zeta (<n>)'
     for all other arguments, including rational noninteger, floating
     point, and complex arguments.

     See also `bfzeta' and `zeta%pi'.

          (%i1) map (zeta, [-4, -3, -2, -1, 0, 1, 2, 3, 4, 5]);
                                               2              4
                     1        1     1       %pi            %pi
          (%o1) [0, ---, 0, - --, - -, inf, ----, zeta(3), ----, zeta(5)]
                    120       12    2        6              90


 -- Option variable: zeta%pi
     Default value: `true'

     When `zeta%pi' is `true', `zeta' returns an expression
     proportional to `%pi^n' for even integer `n'.  Otherwise, `zeta'
     returns a noun form `zeta (n)' for even integer `n'.

          (%i1) zeta%pi: true$
          (%i2) zeta (4);
                                           4
                                        %pi
          (%o2)                         ----
                                         90
          (%i3) zeta%pi: false$
          (%i4) zeta (4);
          (%o4)                        zeta(4)



File: maxima.info,  Node: Symmetries,  Next: Groups,  Prev: Number Theory,  Up: Top

33 Symmetries
*************

* Menu:

* Definitions for Symmetries::