[ruby-talk:4595] Re: ANN: Clifford Algebras in Python

FYI -- in case someone is interested in doing a Ruby variant....

Even if not, the URLs cited below will be interesting for many of you
with an interest in math and physics.

In my non-expert opinion, Geometric calculus:multidimensional math =~=
Ruby:OOP

[e-mail] wrote:
> 
> Hello all three of you who would care about this module!
> 
> At long last, I felt confident enough in my understanding of
> Geometric (Clifford) Algebras to code up an implementation in
> Python (gracefully assisted by NumPy) which I call pyCA.
> (Note: this is numerical, not symbolic.)
> 
> The documentation strings are fairly complete (read: I was too
> lazy to write more documentation without a good reason).  pyCA
> is a fairly general framework for doing work in geometric algebras.
> It's not especially fast for large numbers of computations or
> efficient for many multivectors.  It also concedes to Python's
> operator precedence.  I plan to add modules over the framework
> to implement a convenient shell and implement specific algebras
> and subalgebras (3D Euclidean, quaternions, STA, PGA, etc.).
> 
> Well, have fun with the code below.  I'll place it on the Starship
> once I get around to getting my account re-activated.
> 
> Apologies for the line-breaks.
> 
> ----------%<-------------------------------------------------------
> """\
> pyCA 0.5 -- Geometric Algebra for Python
> Robert Kern
> [e-mail]
> 
> This module implements geometric algebras (a.k.a. Clifford algebras).
> For the uninitiated, a geometric algebra is an algebra of vectors of
> given dimensions and signature. The familiar inner (dot) product and the
> outer product, a generalized relative of the three-dimensional cross
> product are unified in an invertible geometric product.  Scalars, vectors,
> and higher-grade entities can be mixed freely and consistently in the
> form of mixed-grade multivectors.
> 
> For more information, see the Cambridge Clifford Algebras website:
> http://www.mrao.cam.ac.uk/~clifford
> and David Hestenes' website:
> http://modelingnts.la.asu.edu/GC_R&D.html
> 
> This implementation is inspired by Christian Perwass' XGAlib, a C++
> library for GA.  Although there is practically no code shared with
> that LGPLed library, I have still obtained permission from Mr. Perwass
> to release this code into the public domain.  In addition, I have
> obtained permission from James Lehmann to use a snippet of his code
> posted on comp.lang.python .  Thus:
> 
> I hereby release this code into the PUBLIC DOMAIN AS IS.  There is no
> support, warranty, or guarantee.  I will gladly accept comments, bug
> reports, and patches, however.
> 
> Two classes, Layout and MultiVector, and several helper functions are
> provided to implement the algebras.
> 
> Helper Functions
> ================
> 
>   Cl(p, q=0, names=None, firstIdx=0) -- generates the geometric algebra
>       Cl_p,q and returns the appropriate Layout and dictionary of basis
>       element names to their MultiVector instances.
> 
>   bases(layout) -- returns the dictionary of basis element names: instances
>       as above
> 
>   randomMV(layout, grades=None) -- random MultiVector using layout.  If
>       grades is a sequence of integrers, only those grades will be non-zero
> 
>   pretty() -- set repr() to pretty-print MultiVectors
> 
>   ugly() -- set repr() to return eval-able strings
> 
>   eps(newEps) -- set _eps, the epsilon for comparisons and zero-testing
> 
> Issues
> ======
> 
>  * Due to Python's order of operations, the bit operators & ^ << follow
>    the normal arithmetic operators + - * /, so
> 
>      1^e0 + 2^e1  !=  (1^e0) + (2^e1)
> 
>    as is probably intended.  Additionally (comments on this pun will be
>    summarily ignored),
> 
>      M = MultiVector(layout2D)  # null multivector
>      M << 1^e0 << 2^e1 == 10.0^e1 + 1.0^e01
>      M == 1.0
>      e0 == 2 + 1^e0
> 
>    as is definitely not intended.  However,
> 
>      M = MultiVector(layout2D)
>      M << (2^e0) << e1 << (3^e01) == M == 2^e0 + 1^e1 + 3^e01
>      e0 == 1^e0
>      e1 == 1^e1
>      e01 == 1^e01
> 
>    A shell that reorders the operators and adds a few conveniences is in
>    the works.
> 
>  * Since * is the inner product and the inner product with a scalar
>    vanishes by definition, an expression like
> 
>      1*e0 + 2*e1
> 
>    is null.  Use the outer product or full geometric product, to
>    multiply scalars with MultiVectors.  This can cause problems if
>    one has code that mixes Python numbers and MultiVectors.  If the
>    code multiplies two values that can each be either type without
>    checking, one can run into problems as "1 & 2" has a very different
>    result from the same multiplication with scalar MultiVectors.
> 
>  * If the LinearAlgebra module from the NumPy distribution is available,
>    taking the inverse of a MultiVector will use a method proposed by
>    Christian Perwass that involves the solution of a matrix equation.
>    A description of that method follows:
> 
>    Representing multivectors as 2**dims vectors (in the matrix sense),
>    we can carry out the geometric product with a multiplication table.
>    In pseudo-tensorish language (using summation notation):
> 
>      m_i * g_ijk * n_k = v_j
> 
>    Suppose m_i are known (M is the vector we are taking the inverse of),
>    the g_ijk have been computed for this algebra, and v_j = 1 if the
>    j'th element is the scalar element and 0 otherwise, we can compute the
>    dot product m_i * g_ijk.  This yields a rank-2 matrix.  We can
>    then use well-established computational linear algebra techniques
>    to solve this matrix equation for n_k.  The laInv method does precisely
>    that.
> 
>    The usual, analytic, method for computing inverses [M**-1 = ~M/(M&~M) iff
>    M&~M == |M|**2] fails for those multivectors where M&~M is not a scalar.
>    It is only used if the inv method is manually set to point to normalInv
>    or the module LinearAlgebra from NumPy is not available.
> 
>    My testing suggests that laInv works.  In the cases where normalInv works,
>    laInv returns the same result (within _eps).  In all cases,
>    M & M.laInv() == 1.0 (within _eps).  Use whichever you feel comfortable
>    with.
> 
>  * The basis vectors of any algebra will be orthonormal unless you supply
>    your own multiplication tables (which you are free to do after
>    the Layout constructor is called).  A derived class is in preparation
>    that will calculate these tables for you (and include methods for
>    generating reciprocal bases and the like).  I would greatly appreciate
>    help on the algorithm to generate these tables from an arbitrary
>    bilinear form.
> 
>  * No care is taken to preserve the typecode of the arrays.  The purpose
>    of this module is pedagogical.  If your application requires so many
>    multivectors that storage becomes important, the class structure here
>    is unsuitable for you anyways.  Instead, use the algorithms from this
>    module and implement application-specific data structures.
> 
>  * Conversely, explicit typecasting is rare.  MultiVectors will have
>    integer coefficients if you instantiate them that way.  Dividing them
>    by Python integers will have the same consequences as normal integer
>    division.  Public outcry will convince me to add the explicit casts
>    if this becomes a problem.
> 
>  * The way I make variables of the bases is by .update()'ing the global
>    (or local for functions) namespace dictionary.  I use a convenient
>    function for that in the global case.  I have not included it here
>    as I don't want to canonize something that is, well, let's face it,
>    evil.  But it's a venial kind of evil, so do it anyways.
> 
>  * Integration with PyOpenGL for some kick-ass demonstrations are left
>    as an exercise for the reader.
> 
> Happy hacking!
> Robert Kern
> [e-mail]
> """
> 
> import Numeric
> import types
> import operator
> import math
> 
> from Numeric import pi, e
> 
> try:
>     import LinearAlgebra
>     LA = LinearAlgebra
> except ImportError:
>     LA = None
> 
> _eps = 1e-15  # float epsilon for float comparisons
> _pretty = 0   # pretty-print global
> 
> class Layout:
>     """\
> Layout stores information regarding the geometric algebra itself and
> the internal representation of multivectors.  It is constructed like this:
> 
>   Layout(signature, bladeList, firstIdx=0, names=None)
> 
> The arguments to be passed are interpreted as follows:
> 
>   signature -- the signature of the vector space.  This should be
>       a list of positive and negative numbers where the sign determines
>       the sign of the inner product of the corresponding vector with itself.
>       The values are irrelevant except for sign.  This list also determines
>       the dimensionality of the vectors.  Signatures with zeroes are not
>       permitted at this time.
> 
>       Examples:
>         signature = [+1, -1, -1, -1]  # Hestenes', et al. Space-Time Algebra
>         signature = [+1, +1, +1]      # 3-D Euclidean signature
> 
>   bladeList -- list of tuples corresponding to the blades in the whole
>       algebra.  This list determines the order of coefficients in the
>       internal representation of multivectors.  The entry for the scalar
>       must be an empty tuple, and the entries for grade-1 vectors must be
>       singleton tuples.  Remember, the length of the list will be 2**dims.
> 
>       Example:
>         bladeList = [(), (0,), (1,), (0,1)]  # 2-D
> 
>   firstIdx -- the index of the first vector.  That is, some systems number
>       the base vectors starting with 0, some with 1.  Choose by passing
>       the correct number as firstIdx.  0 is the default.
> 
>   names -- list of names of each blade.  When pretty-printing multivectors,
>       use these symbols for the blades.  names should be in the same order
>       as bladeList.  You may use an empty string for scalars.  By default,
>       the name for each non-scalar blade is 'e' plus the indices of the blade
>       as given in bladeList.
> 
>       Example:
>         names = ['', 's0', 's1', 'i']  # 2-D
> 
> 
> Layout's Members:
> 
>   dims -- dimensionality of vectors (== len(signature))
> 
>   sig -- normalized signature (i.e. all values are +1 or -1)
> 
>   firstIdx -- starting point for vector indices
> 
>   bladeList -- list of blades
> 
>   gradeList -- corresponding list of the grades of each blade
> 
>   gaDims -- 2**dims
> 
>   names -- pretty-printing symbols for the blades
> 
>   even -- dictionary of even permutations of blades to the canonical blades
> 
>   odd -- dictionary of odd permutations of blades to the canonical blades
> 
>   gmt -- multiplication table for geometric product [1]
> 
>   imt -- multiplication table for inner product [1]
> 
>   omt -- multiplication table for outer product [1]
> 
> [1] The multiplication tables are NumPy arrays of rank 3 with indices like
>     the tensor g_ijk discussed above.
> """
> 
>     def __init__(self, sig, bladeList, firstIdx=0, names=None):
>         self.dims = len(sig)
>         self.sig = Numeric.divide(sig,
>                                   Numeric.absolute(sig)).astype(Numeric.Int)
>         self.firstIdx = firstIdx
> 
>         self.bladeList = map(tuple, bladeList)
>         self._checkList()
> 
>         self.gaDims = len(self.bladeList)
>         self.gradeList = map(len, self.bladeList)
> 
>         if names is None:
>             e = 'e'
>             self.names = []
> 
>             for i in range(self.gaDims):
>                 if self.gradeList[i] > 1:
>                     self.names.append(e + str(Numeric.add.reduce(map(str, self.bladeList[i]))))
>                 elif self.gradeList[i] == 1:
>                     self.names.append(e + str(self.bladeList[i][0]))
>                 else:
>                     self.names.append('')
> 
>         elif len(names) == self.gaDims:
>             self.names = names
>         else:
>             raise ValueError, "names list of length %i needs to be of length %i" % \
>                                (len(names), self.gaDims)
> 
>         self._genEvenOdd()
>         self._genTables()
> 
>     def __repr__(self):
>         s = "Layout(%s, %s, firstIdx=%s, names=%s)" % (list(self.sig),
>                                           self.bladeList,
>                                           self.firstIdx,
>                                           self.names)
>         return s
> 
>     def _sign(self, seq, orig):
>         """Determine {even,odd}-ness of permutation seq or orig.
> 
>         Returns 1 if even; -1 if odd.
>         """
> 
>         sign = 1
>         seq = list(seq)
> 
>         for i in range(len(seq)):
>             if seq[i] != orig[i]:
>                 j = seq.index(orig[i])
>                 sign = -sign
>                 seq[i], seq[j] = seq[j], seq[i]
>         return sign
> 
>     def _containsDups(self, list):
>         "Checks if list contains duplicates."
> 
>         for k in list:
>             if list.count(k) != 1:
>                 return 1
>         return 0
> 
>     def _genEvenOdd(self):
>         "Create mappings of even and odd permutations to their canonical blades."
> 
>         self.even = {}
>         self.odd = {}
> 
>         class NoMorePermutations(StandardError):
>             pass
> 
>         for i in range(self.gaDims):
>             blade = self.bladeList[i]
>             grade = self.gradeList[i]
> 
>             if grade in (0, 1):
>                 # handled easy cases
>                 self.even[blade] = blade
>                 continue
>             elif grade == 2:
>                 # another easy case
>                 self.even[blade] = blade
>                 self.odd[(blade[1], blade[0])] = blade
>                 continue
>             else:
>                 # general case, lifted from Chooser.py released on
>                 # comp.lang.python by James Lehmann.
>                 idx = range(grade)
> 
>                 try:
>                     for i in range(Numeric.multiply.reduce(range(1, grade+1))):
>                         # grade! permutations
> 
>                         done = 0
>                         j = grade - 1
> 
>                         while not done:
>                             idx[j] = idx[j] + 1
> 
>                             while idx[j] == grade:
>                                 idx[j] = 0
>                                 j = j - 1
>                                 idx[j] = idx[j] + 1
> 
>                                 if j == -1:
>                                     raise NoMorePermutations()
>                             j = grade - 1
> 
>                             if not self._containsDups(idx):
>                                 done = 1
> 
>                         perm = []
> 
>                         for k in idx:
>                             perm.append(blade[k])
> 
>                         perm = tuple(perm)
> 
>                         if self._sign(perm, blade) == 1:
>                             self.even[perm] = blade
>                         else:
>                             self.odd[perm] = blade
> 
>                 except NoMorePermutations:
>                     pass
> 
>                 self.even[blade] = blade
> 
>     def _checkList(self):
>         "Ensure validity of arguments."
> 
>         # check for uniqueness
>         for blade in self.bladeList:
>             if self.bladeList.count(blade) != 1:
>                 raise ValueError, "blades not unique"
> 
>         # check for right dimensionality
>         if len(self.bladeList) != 2**self.dims:
>             raise ValueError, "incorrect number of blades"
> 
>         # check for valid ranges of indices
>         valid =  range(self.firstIdx, self.firstIdx+self.dims)
>         try:
>             for blade in self.bladeList:
>                 for idx in blade:
>                     if (idx not in valid) or (list(blade).count(idx) != 1):
>                         raise ValueError
>         except (ValueError, TypeError):
>             raise ValueError, "invalid bladeList; must be a list of tuples"
> 
>     def _gmtElement(self, a, b):
>         "Returns the element of the geometric multiplication table given blades a, b."
> 
>         mul = 1         # multiplier
> 
>         newBlade = list(a) + list(b)
> 
>         unique = 0
>         while unique == 0:
>             for i in range(len(newBlade)):
>                 if newBlade.count(newBlade[i]) != 1:
>                     delta = newBlade[i+1:].index(newBlade[i])
>                     eo = 1 - 2*(delta % 2)
>                     # eo == 1 if the elements are an even number of flips away
>                     # eo == -1 otherwise
> 
>                     del newBlade[i+delta+1]
>                     del newBlade[i]
>                     mul = eo * mul * self.sig[a[i] - self.firstIdx]
>                     break
>             else:
>                 unique = 1
> 
>         newBlade = tuple(newBlade)
> 
>         if newBlade in self.bladeList:
>             idx = self.bladeList.index(newBlade)
>             # index of the product blade in the bladeList
>         elif newBlade in self.even.keys():
>             # even permutation
>             idx = self.bladeList.index(self.even[newBlade])
>         else:
>             # odd permutation
>             idx = self.bladeList.index(self.odd[newBlade])
>             mul = -mul
> 
>         element = Numeric.zeros((self.gaDims,), Numeric.Int)
>         element[idx] = mul
> 
>         return element, idx
> 
>     def _genTables(self):
>         "Generate the multiplication tables."
> 
>         # geometric multiplication table
>         gmt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
>                              Numeric.Int)
>         # inner product table
>         imt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
>                              Numeric.Int)
> 
>         # outer product table
>         omt = Numeric.zeros((self.gaDims, self.gaDims, self.gaDims),
>                              Numeric.Int)
> 
>         for i in range(self.gaDims):
>             for j in range(self.gaDims):
>                 gmt[i,:,j], idx = self._gmtElement(list(self.bladeList[i]),
>                                                   list(self.bladeList[j]))
> 
>                 if self.gradeList[idx] == abs(self.gradeList[i] - self.gradeList[j]) and \
>                    self.gradeList[i] != 0 and self.gradeList[j] != 0:
> 
>                     # A_r . B_s = <A_r B_s>_|r-s|
>                     # if r,s != 0
>                     imt[i,:,j] = gmt[i,:,j]
> 
>                 if self.gradeList[idx] == (self.gradeList[i] + self.gradeList[j]):
> 
>                     # A_r ^ B_s = <A_r B_s>_|r+s|
>                     omt[i,:,j] = gmt[i,:,j]
> 
>         self.gmt = gmt
>         self.imt = imt
>         self.omt = omt
> 
>     def pseudoScalar(self):
>         "Returns a MultiVector that is the pseudoscalar of this space."
> 
>         psIdx = self.gradeList.index(self.dims)
>         # pick out out blade with grade equal to dims
> 
>         pScalar = MultiVector(self)
>         pScalar.value[psIdx] = 1
> 
>         return pScalar
> 
>     def invPS(self):
>         "Returns the inverse of the pseudoscalar of the algebra."
> 
>         ps = self.pseudoScalar()
> 
>         return ps.inv()
> 
> class MultiVector:
>     """\
> The elements of the algebras, the multivectors, are implemented in the
> MultiVector class.  It has the following constructor:
> 
>   MultiVector(layout, value=None)
> 
> The meaning of the arguments is as follows:
> 
>   layout -- instance of Layout
> 
>   value -- a sequence, of length layout.gaDims, of coefficients of the base
>       blades
> 
> MultiVector's Members
> =====================
> 
>   layout -- instance of Layout
> 
>   value -- a NumPy array, of length layout.gaDims, of coefficients of the base
>       blades
> 
> MultiVector's Public Methods
> ============================
> 
>   __and__, __rand__ -- geometric product.  M & N
> 
>   __xor__, __rxor__ -- outer product.  M ^ N
> 
>   __mul__, __rmul__ -- inner product.  M * N
> 
>   __add__, __radd__ -- addition.  M + N
> 
>   __sub__, __rsub__ -- subtraction.  M - N
> 
>   __div__, __rdiv__ -- division.  M / N <==> M & (N.inv())
> 
>   __pow__ -- exponentiation of a multivector by an integer.  M ** n
> 
>   __rpow__ -- exponentiation of a scalar by a multivector.  x ** M
> 
>   __lshift__ -- in-place addition.  M << N
>                 Adds N to M, modifying M and returning M.
> 
>   mag2() -- magnitude squared.  (M & ~M)[()] == |M|**2
> 
>   __abs__ -- magnitude.  abs(M) == |M|
> 
>   adjoint(), __invert__ -- adjoint/reversion.  ~M
> 
>   __int__, __long__, __float__ -- cast to scalar Python number iff only the
>                                   multivector's scalar coefficient is non-zero
> 
> 
>   __getitem__, __setitem__ -- get/set coefficient to/from a Python number.
>                               Item can be an index into self.value or a
>                               blade (or a permutation of a blade).
> 
>   __delitem__ -- sets coefficient of item (same rules as above) to 0
> 
>   __getslice__ -- returns a new MultiVector with only the selected slice
>                   non-zero
> 
>   __setslice__ -- sets the slice of self.value to a value or values
> 
>   __delslice__ -- sets the slice of self.value to 0
> 
>   __call__ -- projects multivector onto the specified grade.  M(grade)
> 
>   __cmp__ -- returns 0 if the multivectors' coefficients are each equal to
>              within epsilon.  Otherwise, returns a well-defined, but
>              meaningless, ordering.
> 
>   __str__ -- pretty-printed representation
> 
>   __repr__ -- ugly, but eval-able by default, but pretty if global variable
>               _pretty is true
> 
>   isScalar() -- returns 1 if multivector is a scalar
> 
>   isBlade() -- returns 1 if multivector is a blade
> 
>   grades() -- returns the grades contained in the multivector
> 
>   normalInv() -- multivector inverse iff (M & ~M) == |M|**2.  ~M/|M|**2
> 
>   laInv() -- experimental, more general inverse using computational
>              linear algebra to solve for the inverse.
> 
>   inv() -- one of normalInv or laInv.  By default, inv is laInv if the
>            LinearAlgebra module is available, normalInv otherwise.
> 
>   dual() -- returns the dual of the multivector
> 
>   commutator(other) -- returns the commutator of two multivectors
> 
>   join(other) -- returns the join of two blades if it exists
> 
>   meet(other, subspace=None) -- returns the meet of two blades wrt some
>         subspace which defaults to the pseudoscalar
> 
>   gradeInvol() -- returns the grade involution of the multivector
> 
>   conjugate() -- returns the Clifford conjugate of the multivector
> 
> """
> 
>     def __init__(self, layout, value=None):
>         """Constructor.
> 
>         MultiVector(layout, value=None) --> MultiVector
>         """
> 
>         self.layout = layout
> 
>         if value is None:
>             self.value = Numeric.zeros((self.layout.gaDims,), Numeric.Float)
>         else:
>             self.value = Numeric.array(value)
>             if self.value.shape != (self.layout.gaDims,):
>                 raise ValueError, "value must be a sequence of length %s" % self.layout.gaDims
> 
>     def _checkOther(self, other, coerce=1):
>         """Ensure that the other argument has the same Layout or coerce value if
>         necessary/requested.
> 
>         _checkOther(other, coerce=1) --> newOther, isMultiVector
>         """
> 
>         if type(other) in (types.IntType, types.FloatType, types.LongType):
>             if coerce:
>                 # numeric scalar
>                 newOther = MultiVector(self.layout)
>                 newOther[()] = other
>                 return newOther, 1
>             else:
>                 return other, 0
> 
>         elif isinstance(other, self.__class__) and \
>              other.layout is not self.layout:
>             raise ValueError, "cannot operate on MultiVectors with different Layouts"
> 
>         return other, 1
> 
>     ## numeric special methods
>     # binary
> 
>     def __and__(self, other):
>         """Geometric product
> 
>         M & N --> MN
>         __and__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other, coerce=0)
> 
>         if mv:
>             newValue = Numeric.dot(self.value, Numeric.dot(self.layout.gmt,
>                                                            other.value))
>         else:
>             newValue = other * self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __rand__(self, other):
>         """Right-hand geometric product
> 
>         N & M --> NM
>         __rand__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other, coerce=0)
> 
>         if mv:
>             newValue = Numeric.dot(other.value, Numeric.dot(self.layout.gmt,
>                                                             self.value))
>         else:
>             newValue = other * self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __xor__(self, other):
>         """Outer product
> 
>         M ^ N
>         __xor__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other, coerce=0)
> 
>         if mv:
>             newValue = Numeric.dot(self.value, Numeric.dot(self.layout.omt,
>                                                            other.value))
>         else:
>             newValue = other * self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __rxor__(self, other):
>         """Right-hand outer product
> 
>         N ^ M
>         __rxor__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other, coerce=0)
> 
>         if mv:
>             newValue = Numeric.dot(other.value, Numeric.dot(self.layout.omt,
>                                                             self.value))
>         else:
>             newValue = other * self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __mul__(self, other):
>         """Inner product
> 
>         M * N
>         __mul__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         if mv:
>             newValue = Numeric.dot(self.value, Numeric.dot(self.layout.imt,
>                                                            other.value))
>         else:
>             return MultiVector(self.layout)  # l * M = M * l = 0 for scalar l
> 
>         return MultiVector(self.layout, newValue)
> 
>     __rmul__ = __mul__
> 
>     def __add__(self, other):
>         """Addition
> 
>         M + N
>         __add__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
>         newValue = self.value + other.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     __radd__ = __add__
> 
>     def __sub__(self, other):
>         """Subtraction
> 
>         M - N
>         __sub__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
>         newValue = self.value - other.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __rsub__(self, other):
>         """Right-hand subtraction
> 
>         N - M
>         __rsub__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
>         newValue = other.value - self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __div__(self, other):
>         """Division
>                        -1
>         M / N --> M & N
>         __div__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other, coerce=0)
> 
>         if mv:
>             return self & other.inv()
>         else:
>             newValue = self.value / other
>             return MultiVector(self.layout, newValue)
> 
>     def __rdiv__(self, other):
>         """Right-hand division
>                        -1
>         N / M --> N & M
>         __rdiv__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         return other & self.inv()
> 
>     def __pow__(self, other):
>         """Exponentiation of a multivector by an integer
>                     n
>         M ** n --> M
>         __pow__(other) --> MultiVector
>         """
> 
>         if type(other) not in (types.IntType, types.FloatType):
>             raise ValueError, "exponent must be a Python Int or Float"
> 
>         if abs(round(other) - other) > _eps:
>             raise ValueError, "exponent must have no fractional part"
> 
>         other = int(round(other))
> 
>         newMV = MultiVector(self.layout, list(self.value))
> 
>         for i in range(1, other):
>             newMV = newMV & self
> 
>         return newMV
> 
>     def __rpow__(self, other):
>         """Exponentiation of a real by a multivector
>                   M
>         r**M --> r
>         __rpow__(other) --> MultiVector
>         """
> 
>         # check that other is a Python number, not a MultiVector
> 
>         if type(other) not in (types.IntType, types.FloatType, types.LongType):
>             raise ValueError, "base must be a Python number"
> 
>         intMV = math.log(other) & self
>         # pow(x, y) == exp(y & log(x))
> 
>         newMV = MultiVector(self.layout)  # null
> 
>         nextTerm = MultiVector(self.layout)
>         nextTerm[()] = 1  # order 0 term of exp(x) Taylor expansion
> 
>         n = 1
> 
>         while nextTerm != 0:
>             # iterate until the added term is within _eps of 0
>             newMV << nextTerm
>             nextTerm = nextTerm & intMV / n
>             n = n + 1
>         else:
>             # squeeze out that extra little bit of accuracy
>             newMV << nextTerm
> 
>         return newMV
> 
>     def __lshift__(self, other):
>         """In-place addition
> 
>         M << N --> M + N
>         __lshift__(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         self.value = self.value + other.value
> 
>         return self
> 
> 
>     # unary
> 
>     def __neg__(self):
>         """Negation
> 
>         -M
>         __neg__() --> MultiVector
>         """
> 
>         newValue = -self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def __pos__(self):
>         """Positive (just a copy)
> 
>         +M
>         __pos__(self) --> MultiVector
>         """
> 
>         newValue = self.value + 0  # copy
> 
>         return MultiVector(self.layout, newValue)
> 
>     def mag2(self):
>         """Magnitude (modulus) squared
>            2
>         |M|
>         mag2() --> PyFloat | PyInt
>         """
> 
>         return (~self & self)[()]
> 
>     def __abs__(self):
>         """Magnitude (modulus)
> 
>         abs(M) --> |M|
>         __abs__() --> PyFloat
>         """
> 
>         return Numeric.sqrt(self.mag2())
> 
>     def adjoint(self):
>         """Adjoint / reversion
>                _
>         ~M --> M (any one of several conflicting notations)
>         ~(N & M) --> ~M & ~N
>         adjoint() --> MultiVector
>         """
>         # The multivector created by reversing all multiplications
> 
>         newMV = MultiVector(self.layout)
> 
>         grades = Numeric.array(self.layout.gradeList)
>         signs = Numeric.power(-1, grades*(grades-1)/2)
> 
>         newMV.value = signs * self.value
> 
>         return newMV
> 
>     __invert__ = adjoint
> 
>     # builtin
> 
>     def __int__(self):
>         """Coerce to an integer iff scalar.
> 
>         int(M)
>         __int__() --> PyInt
>         """
> 
>         return int(self.__float__())
> 
>     def __long__(self):
>         """Coerce to a long iff scalar.
> 
>         long(M)
>         __long__() --> PyLong
>         """
> 
>         return long(self.__float__())
> 
>     def __float__(self):
>         """"Coerce to a float iff scalar.
> 
>         float(M)
>         __float__() --> PyFloat
>         """
> 
>         if self.isScalar():
>             return float(self[()])
>         else:
>             raise ValueError, "non-scalar coefficients are non-zero"
> 
> 
>     ## sequence special methods
> 
>     def __len__(self):
>         """Returns length of value array.
> 
>         len(M)
>         __len__() --> PyInt
>         """
> 
>         return self.layout.gaDims
> 
>     def __getitem__(self, key):
>         """If key is a blade tuple (e.g. (0,1) or (1,3)), then return
>         the (real) value of that blade's coefficient.
>         Otherwise, treat key as an index into the list of coefficients.
> 
>         M[blade] --> PyFloat | PyInt
>         M[index] --> PyFloat | PyInt
>         __getitem__(key) --> PyFloat | PyInt
>         """
> 
>         if key in self.layout.bladeList:
>             return self.value[self.layout.bladeList.index(key)]
>         elif key in self.layout.even.keys():
>             return self.value[self.layout.bladeList.index(self.layout.even[key])]
>         elif key in self.layout.odd.keys():
>             return -self.value[self.layout.bladeList.index(self.layout.odd[key])]
>         else:
>             return self.value[key]
> 
>     def __setitem__(self, key, value):
>         """If key is a blade tuple (e.g. (0,1) or (1,3)), then set
>         the (real) value of that blade's coeficient.
>         Otherwise treat key as an index into the list of coefficients.
> 
>         M[blade] = PyFloat | PyInt
>         M[index] = PyFloat | PyInt
>         __setitem__(key, value)
>         """
> 
>         if key in self.layout.bladeList:
>             self.value[self.layout.bladeList.index(key)] = value
>         elif key in self.layout.even.keys():
>             self.value[self.layout.bladeList.index(self.layout.even[key])] = value
>         elif key in self.layout.odd.keys():
>             self.value[self.layout.bladeList.index(self.layout.odd[key])] = -value
>         else:
>             self.value[key] = value
> 
>     def __delitem__(self, key):
>         """Set the selected coefficient to 0.
> 
>         del M[blade]
>         del M[index]
>         __delitem__(key)
>         """
> 
>         if key in self.layout.bladeList:
>             self.value[self.layout.bladeList.index(key)] = 0
>         elif key in self.layout.even.keys():
>             self.value[self.layout.bladeList.index(self.layout.even[key])] = 0
>         elif key in self.layout.odd.keys():
>             self.value[self.layout.bladeList.index(self.layout.odd[key])] = 0
>         else:
>             self.value[key] = 0
> 
>     def __getslice__(self, i, j):
>         """Return a copy with only the slice non-zero.
> 
>         M[i:j]
>         __getslice__(i, j) --> MultiVector
>         """
> 
>         newMV = MultiVector(self.layout)
>         newMV.value[i:j] = self.value[i:j]
> 
>         return newMV
> 
>     def __setslice__(self, i, j, sequence):
>         """Paste a sequence into coefficients array.
> 
>         M[i:j] = sequence
>         __setslice__(i, j, sequence)
>         """
> 
>         self.value[i:j] = sequence
> 
>     def __delslice__(self, i, j):
>         """Set slice to zeros.
> 
>         del M[i:j]
>         __delslice__(i, j)
>         """
> 
>         self.value[i:j] = 0
> 
>     ## grade projection
> 
>     def __call__(self, grade):
>         """Return a new multi-vector projected onto the specified grade.
> 
>         M(grade) --> <M>
>                         grade
>         __call__(grade) --> MultiVector
>         """
> 
>         if grade not in self.layout.gradeList:
>             raise ValueError, "algebra does not have grade %s" % grade
> 
>         if type(grade) is not types.IntType:
>             raise ValueError, "grade must be an integer"
> 
>         mask = Numeric.equal(grade, self.layout.gradeList)
> 
>         newValue = Numeric.multiply(mask, self.value)
> 
>         return MultiVector(self.layout, newValue)
> 
>     ## fundamental special methods
> 
>     def __str__(self):
>         """Return pretty-printed representation.
> 
>         str(M)
>         __str__() --> PyString
>         """
> 
>         s = ''
> 
>         for i in range(self.layout.gaDims):
>             # if we have nothing yet, don't use + and - as operators but
>             # use - as an unary prefix if necessary
>             if s:
>                 seps = (' + ', ' - ')
>             else:
>                 seps = ('', '-')
> 
>             if self.layout.gradeList[i] == 0:
>                 # scalar
>                 if abs(self.value[i]) >= _eps:
>                     if self.value[i] > 0:
>                         s = '%s%s%s' % (s, seps[0], self.value[i])
>                     else:
>                         s = '%s%s%s' % (s, seps[1], -self.value[i])
> 
>             else:
>                 if abs(self.value[i]) >= _eps:
>                     # not a scalar
>                     if self.value[i] > 0:
>                         s = '%s%s%s^%s' % (s, seps[0], self.value[i],
>                                               self.layout.names[i])
>                     else:
>                         s = '%s%s%s^%s' % (s, seps[1], -self.value[i],
>                                               self.layout.names[i])
>         if s:
>             # non-zero
>             return s
>         else:
>             # return scalar 0
>             return '0'
> 
>     def __repr__(self):
>         """Return eval-able representation if global _pretty is false.
>         Otherwise, return str(self).
> 
>         repr(M)
>         __repr__() --> PyString
>         """
> 
>         if _pretty:
>             return self.__str__()
> 
>         s = "MultiVector(%s, value=%s)" % \
>              (repr(self.layout), list(self.value))
>         return s
> 
>     def __nonzero__(self):
>         """Instance is nonzero iff at least one of the coefficients
>         is nonzero.
> 
>         __nonzero() --> Boolean
>         """
> 
>         nonzeroes = Numeric.greater(Numeric.absolute(self.value), _eps)
> 
>         if nonzeroes:
>             return 1
>         else:
>             return 0
> 
>     def __cmp__(self, other):
>         """Compares two multivectors.
> 
>         This is mostly defined for equality testing (within epsilon).
>         In the case that the two multivectors have different Layouts,
>         we will raise an error.  In the case that they are not equal,
>         we will compare the tuple represenations of the coefficients
>         lists just so as to return something valid.  Therefore,
>         inequalities are well-nigh meaningless (since they are
>         meaningless for multivectors while equality is meaningful).
>         Oh, how I wish for rich comparisons.
> 
>         M == N
>         __cmp__(other) --> -1|0|1
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         if Numeric.alltrue(Numeric.less(Numeric.absolute( \
>                 self.value - other.value), _eps)):
>             # equal within epsilon
>             return 0
>         else:
>             return cmp(tuple(self.value), tuple(other.value))
> 
>     ## Geometric Algebraic functions
> 
>     def isScalar(self):
>         """Returns true iff self is a scalar.
> 
>         isScalar() --> Boolean
>         """
> 
>         indices = range(self.layout.gaDims)
>         indices.remove(self.layout.gradeList.index(0))
> 
>         for i in indices:
>             if abs(self.value[i]) < _eps:
>                 continue
>             else:
>                 return 0
> 
>         return 1
> 
>     def isBlade(self):
>         """Returns true if multivector is a blade.
> 
>         isBlade() --> Boolean
>         """
> 
>         grade = None
> 
>         for i in range(self.layout.gaDims):
>             if abs(self.value[i]) > _eps:
>                 if grade is None:
>                     grade = self.layout.gradeList[i]
>                 elif self.layout.gradeList[i] != grade:
>                     return 0
> 
>         return 1
> 
>     def grades(self):
>         """Return the grades contained in the multivector.
> 
>         grades() --> PyInt"""
> 
>         grades = []
> 
>         for i in range(self.layout.gaDims):
>             if abs(self.value[i]) > _eps and \
>                self.layout.gradeList[i] not in grades:
>                 grades.append(self.layout.gradeList[i])
> 
>         return grades
> 
>     def laInv(self):
>         """Return inverse using a computational linear algebra method
>         proposed by Christian Perwass.
>          -1
>         M
>         laInv() --> MultiVector
>         """
> 
>         identity = Numeric.zeros((self.layout.gaDims,))
>         identity[self.layout.gradeList.index(0)] = 1
> 
>         intermed = Numeric.dot(self.value, self.layout.gmt)
>         intermed = Numeric.transpose(intermed)
> 
>         if abs(LA.determinant(intermed)) < _eps:
>             raise ValueError, "multivector has no inverse"
> 
>         sol = LA.solve_linear_equations(intermed, identity)
> 
>         return MultiVector(self.layout, sol)
> 
> 
>     def normalInv(self):
>         """Returns the inverse of itself if M&~M == |M|**2.
>          -1
>         M   = ~M / (M & ~M)
>         normalInv() --> MultiVector
>         """
> 
>         Madjoint = ~self
>         MadjointM = (Madjoint & self)
> 
>         if MadjointM.isScalar() and abs(MadjointM[()]) > _eps:
>             # inverse exists
>             return Madjoint / MadjointM[()]
>         else:
>             raise ValueError, "no inverse exists for this multivector"
>     if LA:
>         inv = laInv
>     else:
>         inv = normalInv
> 
>     def dual(self):
>         """Returns the dual of the multivector.
> 
>         ~        -1
>         M = M * I     where I is the pseudoscalar.
>         dual() --> MultiVector
>         """
> 
>         return self * self.layout.invPS()
> 
>     def commutator(self, other):
>         """Returns the commutator product of two multivectors.
> 
>         [M, N] = M X N = (M&N - N&M)/2
>         commutator(other) --> MultiVector
>         """
> 
>         return ((self & other) - (other & self)) / 2
> 
>     def join(self, other):
>         """Returns the join of two blades.
> 
>         J = M ^ N   iff M, N are blades and M^N != 0
>         join(other) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         if self.isBlade() and other.isBlade():
>             J = self ^ other
> 
>             if J:
>                 return J
>             else:
>                 raise ValueError, "join does not exist, blades share common factors"
>         else:
>             raise ValueError, "not blades"
> 
>     def meet(self, other, subspace=None):
>         """Returns the meet of an r-blade and an s-blade iff r+s >= dims.
> 
>         Computation is done with respect to a subspace that defaults to
>         the pseudo scalar if none is given.
>                      -1
>         M \/i N = (Mi  ) * N
>         meet(other, subspace=None) --> MultiVector
>         """
> 
>         other, mv = self._checkOther(other)
> 
>         r = self.grades()
>         s = other.grades()
> 
>         if len(r) > 1 or len(s) > 1:
>             raise ValueError, "not blades"
> 
>         r = r[0]
>         s = s[0]
> 
>         if subspace is None:
>             subspace = self.layout.pseudoScalar()
> 
>         dims = subspace.grades()[0]
> 
>         if r + s < dims:
>             raise ValueError, "r+s < dims; meet not defined"
>         elif r + s == dims:
>             # special case; M v N == 0
>             return MultiVector(self.layout)
>         else:
>             return (self & subspace.inv()) * other
> 
>     def gradeInvol(self):
>         """Returns the grade involution of the multivector.
>          *                    i
>         M  = Sum[i, dims, (-1)  <M>  ]
>                                    i
>         gradeInvol() --> MultiVector
>         """
> 
>         signs = Numeric.power(-1, self.layout.gradeList)
> 
>         newValue = signs * self.value
> 
>         return MultiVector(self.layout, newValue)
> 
>     def conjugate(self):
>         """Returns the Clifford conjugate (reversion and grade involution).
>          *
>         M  --> (~M).gradeInvol()
>         conjugate() --> MultiVector
>         """
> 
>         return (~self).gradeInvol()
> 
> def comb(n, k):
>     """\
>     Returns /n\\
>             \\k/
> 
>     comb(n, k) --> PyInt"""
> 
>     def fact(n):
>         return Numeric.multiply.reduce(range(1, n+1))
> 
>     return fact(n) / (fact(k) * fact(n-k))
> 
> def elements(dims, firstIdx=0):
>     """Return a list of tuples representing all 2**dims of blades
>     in a dims-dimensional GA.
> 
>     elements(dims, firstIdx=0) --> bladeList"""
> 
>     indcs = range(firstIdx, firstIdx + dims)
> 
>     blades = [()]
> 
>     for k in range(1, dims+1):
>         # k = grade
> 
>         if k == 1:
>             for i in indcs:
>                 blades.append((i,))
>             continue
> 
>         curBladeX = indcs[:k]
> 
>         marker = -2
>         # curBladeX[marker:] will be replaced with a range each time
>         # the final element overflows; you'll see
> 
>         for i in range(comb(dims, k)):
>             if curBladeX[-1] < firstIdx+dims-1:
>                 # increment last index
>                 blades.append(tuple(curBladeX))
>                 curBladeX[-1] = curBladeX[-1] + 1
> 
>             else:
>                 if curBladeX[marker:] == range(firstIdx+dims+marker,
>                                                firstIdx+dims):
>                     # decrement marker
>                     marker = marker - 1
>                     if marker < -k:
>                         blades.append(tuple(curBladeX))
>                         continue
> 
>                 # replace
>                 blades.append(tuple(curBladeX))
>                 curBladeX[marker:] = range(curBladeX[marker]+1,
>                                            curBladeX[marker]+1 - marker)
> 
>     return blades
> 
> 
> def Cl(p, q=0, names=None, firstIdx=0):
>     """Returns a Layout and basis blades for the geometric algebra Cl_p,q.
> 
>     The notation Cl_p,q means that the algebra is p+q dimensional, with
>     the first p vectors with positive signature and the final q vectors
>     negative.
> 
>     Cl(p, q=0, names=None, firstIdx=0) --> Layout, {'name': basisElement, ...}"""
> 
>     sig = [+1]*p + [-1]*q
>     bladeList = elements(p+q, firstIdx)
> 
>     layout = Layout(sig, bladeList, firstIdx=firstIdx, names=names)
>     blades = bases(layout)
> 
>     return layout, blades
> 
> def bases(layout):
>     """Returns a dictionary mapping basis element names to their MultiVector
>     instances.
> 
>     bases(layout) --> {'name': baseElement, ...}"""
> 
>     dict = {}
>     for i in range(layout.gaDims):
>         if layout.gradeList[i] != 0:
>             v = Numeric.zeros((layout.gaDims,))
>             v[i] = 1
>             dict[layout.names[i]] = MultiVector(layout, v)
>     return dict
> 
> def randomMV(layout, min=-2.0, max=2.0, grades=None):
>     """Random MultiVector with given layout.
> 
>     Coefficients are between min and max, and if grades is a list of integers,
>     only those grades will be non-zero.  Uses the RandomArray module.
> 
>     randomMV(layout, min=-2.0, max=2.0, grades=None)"""
> 
>     from RandomArray import uniform
> 
>     if grades is None:
>         return MultiVector(layout, uniform(min, max, (layout.gaDims,)))
>     else:
>         newValue = Numeric.zeros((layout.gaDims,)).astype(Numeric.Float)
>         for i in range(layout.gaDims):
>             if layout.gradeList[i] in grades:
>                 newValue[i] = uniform(min, max)
>         return MultiVector(layout, newValue)
> 
> def pretty():
>     """Makes repr(M) default to pretty-print.
> 
>     pretty()"""
> 
>     global _pretty
>     _pretty = 1
> 
> def ugly():
>     """Makes repr(M) default to eval-able representation.
> 
>     ugly()"""
> 
>     global _pretty
>     _pretty = 0
> 
> def eps(newEps):
>     """Set the epsilon for float comparisons.
> 
>     eps(newEps)"""
> 
>     global _eps
>     _eps = newEps
> 
> ------------%<------------------------------------------
> 
> --
> Robert Kern
> [e-mail]
> 
> "In the fields of hell where the grass grows high
>  Are the graves of dreams allowed to die."
>   -- Richard Harter


-- 
Conrad Schneiker
(This note is unofficial and subject to improvement without notice.)