# Commits

committed b23d5df

First checkin.

• Participants
• Branches default

# File clifford.py

`+""" Geometric Algebra for Python`
`+`
`+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`
`+`
`+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, mvClass=MultiVector) -- generates `
`+      the geometric algebra Cl_p,q and returns the appropriate Layout `
`+      and dictionary of basis element names to their MultiVector instances.`
`+      Uses mvClass if provided.`
`+`
`+  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. Contributions, of course, are appreciated! I'm thinking`
`+   about using parsetools from`
`+   ftp://archaeopteryx.com/pub/parsetools/parsetools-0.9.2.tgz`
`+   `
`+ * 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.`
`+`
`+   Of course, a new issue arises with this method.  The inverses found`
`+   are sometimes dependant on the order of multiplication.  That is:`
`+`
`+     M.laInv() & M == 1.0`
`+     M & M.laInv() != 1.0`
`+   `
`+   XXX Thus, there are two other methods defined, leftInv and rightInv`
`+   which point to leftLaInv and rightLaInv if the LinearAlgebra module`
`+   exists.  The method inv points to rightInv.  Should LinearAlgebra`
`+   not be available, or the user choose, leftInv and rightInv will both`
`+   point to normalInv, which yields a left- and right-inverse that are`
`+   the same should either exist (the proof is fairly simple).`
`+   `
`+ * 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. Pretty please?`
`+`
`+ * 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.`
`+`
`+`
`+Changes 0.6-0.7`
`+===============`
`+`
`+ * Added a real license.`
`+ * Convert to numpy instead of Numeric.`
`+   `
`+Changes 0.5-0.6`
`+===============`
`+`
`+ * join() and meet() actually work now, but have numerical accuracy problems`
`+ * added clean() to MultiVector`
`+ * added leftInv() and rightInv() to MultiVector`
`+ * moved pseudoScalar() and invPS() to MultiVector (so we can derive `
`+   new classes from MultiVector)`
`+ * changed all of the instances of creating a new MultiVector to create`
`+   an instance of self.__class__ for proper inheritance`
`+ * fixed bug in laInv()`
`+ * fixed the massive confusion about how dot() works`
`+ * added left-contraction`
`+ * fixed embarassing bug in gmt generation`
`+ * added normal() and anticommutator() methods`
`+ * fixed dumb bug in elements() that limited it to 4 dimensions`
`+`
`+Happy hacking!`
`+Robert Kern`
`+robert.kern@gmail.com`
`+"""`
`+`
`+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`
`+`
`+def _myDot(a, b):`
`+    """Returns the inner product as *I* learned it.`
`+`
`+    a_i...k * b_k...m = c_i...m     in summation notation with the ...'s `
`+                                    representing arbitrary, omitted indices`
`+`
`+    The sum is over the last axis of the first argument and the first axis `
`+    of the last axis.`
`+`
`+    _myDot(a, b) --> NumPy array`
`+    """`
`+`
`+    a = Numeric.asarray(a)`
`+    b = Numeric.asarray(b)`
`+`
`+    tempAxes = tuple(range(1, len(b.shape)) + [0])`
`+    newB = Numeric.transpose(b, tempAxes)`
`+`
`+    # innerproduct sums over the *last* axes of *both* arguments`
`+    return Numeric.innerproduct(a, newB)`
`+    `
`+`
`+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]`
`+`
`+  lcmt -- multiplication table for the left-contraction [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 with permission.`
`+                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)):`
`+                index = newBlade[i]`
`+                if newBlade.count(index) != 1:`
`+                    delta = newBlade[i+1:].index(index)`
`+                    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[index - 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)`
`+`
`+        # left-contraction table`
`+        lcmt = 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]`
`+`
`+                if self.gradeList[idx] == (self.gradeList[j] - self.gradeList[i]):`
`+                    # A_r _| B_s = <A_r B_s>_(s-r) if s-r >= 0`
`+                    lcmt[i,:,j] = gmt[i,:,j]`
`+                    `
`+        self.gmt = gmt`
`+        self.imt = imt`
`+        self.omt = omt`
`+        self.lcmt = lcmt`
`+`
`+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`
`+    """`
`+`
`+    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 = self._newMV()`
`+                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`
`+`
`+    def _newMV(self, newValue=None):`
`+        """Returns a new MultiVector (or derived class instance).`
`+`
`+        _newMV(self, newValue=None)`
`+        """`
`+`
`+        return self.__class__(self.layout, newValue)`
`+`
`+`
`+    ## 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 self._newMV(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 self._newMV(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 self._newMV(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 self._newMV(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 self._newMV()  # l * M = M * l = 0 for scalar l`
`+`
`+        return self._newMV(newValue)`
`+`
`+    __rmul__ = __mul__`
`+    `
`+    def __add__(self, other):`
`+        """Addition`
`+        `
`+        M + N`
`+        __add__(other) --> MultiVector`
`+        """`
`+        `
`+        other, mv = self._checkOther(other)`
`+        newValue = self.value + other.value`
`+`
`+        return self._newMV(newValue)`
`+        `
`+    __radd__ = __add__`
`+`
`+    def __sub__(self, other):`
`+        """Subtraction`
`+        `
`+        M - N`
`+        __sub__(other) --> MultiVector`
`+        """`
`+        `
`+        other, mv = self._checkOther(other)`
`+        newValue = self.value - other.value`
`+`
`+        return self._newMV(newValue)`
`+`
`+    def __rsub__(self, other):`
`+        """Right-hand subtraction`
`+        `
`+        N - M`
`+        __rsub__(other) --> MultiVector`
`+        """`
`+        `
`+        other, mv = self._checkOther(other)`
`+        newValue = other.value - self.value`
`+`
`+        return self._newMV(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 self._newMV(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 = self._newMV(Numeric.array(self.value))  # copy`
`+`
`+        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 = self._newMV()  # null`
`+`
`+        nextTerm = self._newMV()`
`+        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 self._newMV(newValue)`
`+`
`+    def __pos__(self):`
`+        """Positive (just a copy)`
`+`
`+        +M`
`+        __pos__(self) --> MultiVector`
`+        """`
`+        `
`+        newValue = self.value + 0  # copy`
`+`
`+        return self._newMV(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`
`+`
`+        grades = Numeric.array(self.layout.gradeList)`
`+        signs = Numeric.power(-1, grades*(grades-1)/2)`
`+`
`+        newValue = signs * self.value  # elementwise multiplication`
`+        `
`+        return self._newMV(newValue)`
`+`
`+    __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 = self._newMV()`
`+        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 self._newMV(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))`
`+`
`+    def clean(self, eps=None):`
`+        """Sets coefficients whose absolute value is < eps to exactly 0.`
`+`
`+        eps defaults to the current value of the global _eps.`
`+`
`+        clean(eps=None)`
`+        """`
`+`
`+        if eps is None:`
`+            eps = _eps`
`+`
`+        mask = Numeric.greater(Numeric.absolute(self.value), eps)`
`+`
`+        # note element-wise multiplication`
`+        self.value = mask * self.value`
`+`
`+        return self`
`+`
`+    def round(self, eps=None):`
`+        """Rounds all coefficients according to Python's rounding rules.`
`+`
`+        eps defaults to the current value of the global _eps.`
`+`
`+        round(eps=None)`
`+        """`
`+`
`+        if eps is None:`
`+            eps = _eps`
`+`
`+        `
`+        self.value = Numeric.around(self.value, eps)`
`+`
`+        return self`
`+`
`+    ## Geometric Algebraic functions`
`+`
`+    def lc(self, other):`
`+        """Returns the left-contraction of two multivectors.`
`+`
`+        M _| N`
`+        lc(other) --> MultiVector`
`+        """`
`+`
`+        other, mv = self._checkOther(other, coerce=1)`
`+`
`+        newValue = Numeric.dot(self.value, Numeric.dot(self.layout.lcmt,`
`+                                                       other.value))`
`+`
`+        return self._newMV(newValue)`
`+`
`+    `
`+    `
`+    def pseudoScalar(self):`
`+        "Returns a MultiVector that is the pseudoscalar of this space."`
`+`
`+        psIdx = self.layout.gradeList.index(self.layout.dims)  `
`+        # pick out out blade with grade equal to dims`
`+`
`+        pScalar = self._newMV()`
`+        pScalar.value[psIdx] = 1`
`+`
`+        return pScalar`
`+`
`+    def invPS(self):`
`+        "Returns the inverse of the pseudoscalar of the algebra."`
`+`
`+        ps = self.pseudoScalar()`
`+`
`+        return ps.inv()`
`+`
`+    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, 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 normal(self):`
`+        """Return the (mostly) normalized multivector.`
`+`
`+        The _mostly_ comes from the fact that some multivectors have a `
`+        negative squared-magnitude.  So, without introducing formally`
`+        imaginary numbers, we can only fix the normalized multivector's`
`+        magnitude to +-1.`
`+        `
`+        M / |M|  up to a sign`
`+        normal() --> MultiVector`
`+        """`
`+`
`+        return self / Numeric.sqrt(abs(self.mag2()))`
`+        `
`+    def leftLaInv(self):`
`+        """Return left-inverse using a computational linear algebra method `
`+        proposed by Christian Perwass.`
`+         -1         -1`
`+        M    where M  & M  == 1`
`+        leftLaInv() --> MultiVector`
`+        """`
`+        `
`+        identity = Numeric.zeros((self.layout.gaDims,))`
`+        identity[self.layout.gradeList.index(0)] = 1`
`+`
`+        intermed = Numeric.dot(self.layout.gmt, self.value)`
`+        intermed = Numeric.transpose(intermed)`
`+`
`+        if abs(LA.determinant(intermed)) < _eps:`
`+            raise ValueError, "multivector has no left-inverse"`
`+`
`+        sol = LA.solve_linear_equations(intermed, identity)`
`+`
`+        return self._newMV(sol)`
`+        `
`+    def rightLaInv(self):`
`+        """Return right-inverse using a computational linear algebra method `
`+        proposed by Christian Perwass.`
`+         -1              -1`
`+        M    where M & M  == 1`
`+        rightLaInv() --> MultiVector`
`+        """`
`+`
`+        identity = Numeric.zeros((self.layout.gaDims,))`
`+        identity[self.layout.gradeList.index(0)] = 1`
`+`
`+        intermed = _myDot(self.value, self.layout.gmt)`
`+`
`+        if abs(LA.determinant(intermed)) < _eps:`
`+            raise ValueError, "multivector has no right-inverse"`
`+`
`+        sol = LA.solve_linear_equations(intermed, identity)`
`+`
`+        return self._newMV(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:`
`+        leftInv = leftLaInv`
`+        inv = rightInv = rightLaInv`
`+    else:`
`+        inv = leftInv = rightInv = normalInv`
`+`
`+    def dual(self, I=None):`
`+        """Returns the dual of the multivector against the given subspace I.`
`+        I defaults to the pseudoscalar.`
`+`
`+        ~        -1`
`+        M = M * I`
`+        dual(I=None) --> MultiVector`
`+        """`
`+        `
`+        if I is None:`
`+            Iinv = self.layout.invPS()`
`+        else:`
`+            Iinv = I.inv()`
`+        `
`+        return self * Iinv`
`+`
`+    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 anticommutator(self, other):`
`+        """Returns the anti-commutator product of two multivectors.`
`+`
`+        (M&N + N&M)/2`
`+        anticommutator(other) --> MultiVector`
`+        """`
`+`
`+        return ((self & other) + (other & self)) / 2`
`+`
`+    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 self._newMV(newValue)`
`+`
`+    def conjugate(self):`
`+        """Returns the Clifford conjugate (reversion and grade involution).`
`+         *`
`+        M  --> (~M).gradeInvol()`
`+        conjugate() --> MultiVector`
`+        """`
`+`
`+        return (~self).gradeInvol()`
`+`
`+    ## Subspace operations`
`+`
`+    def project(self, other):`
`+        """Projects the multivector onto the subspace represented by this blade.`
`+                            -1`
`+        P (M) = (M _| A) & A`
`+         A`
`+        project(M) --> MultiVector`
`+        """`
`+`
`+        other, mv = self._checkOther(other, coerce=1)`
`+`
`+        if not self.isBlade():`
`+            raise ValueError, "self is not a blade"`
`+`
`+        return other.lc(self) & self.inv()`
`+`
`+    def basis(self):`
`+        """Finds a vector basis of this subspace.`
`+`
`+        basis() --> [ MultiVector, MultiVector, ... ]`
`+        """`
`+`
`+        gr = self.grades()`
`+`
`+        if len(gr) != 1:`
`+            raise ValueError, "self is not a blade"`
`+`
`+        selfInv = self.inv()`
`+`
`+        selfInv.clean()`
`+`
`+        wholeBasis = []  # vector basis of the whole space`
`+`
`+        for i in range(self.layout.gaDims):`
`+            if self.layout.gradeList[i] == 1:`
`+                v = Numeric.zeros((self.layout.gaDims,), Numeric.Float)`
`+                v[i] = 1.`
`+                wholeBasis.append(self._newMV(v))`
`+`
`+        thisBasis = []  # vector basis of this subspace`
`+`
`+        J, mv = self._checkOther(1.)  # outer product of all of the vectors up `
`+                             # to the point of iteration`
`+`
`+        for ei in wholeBasis:`
`+            Pei = ei.lc(self) & selfInv`
`+`
`+            J.clean()`
`+            `
`+            J2 = J ^ Pei`
`+`
`+            if J2 != 0:`
`+                J = J2`
`+                thisBasis.append(Pei)`
`+                if len(thisBasis) == gr[0]:  # we have a complete set`
`+                    break`
`+        `
`+        return thisBasis`
`+`
`+    def join(self, other):`
`+        """Returns the join of two blades.`
`+              .`
`+        J = A ^ B`
`+        join(other) --> MultiVector`
`+        """`
`+`
`+        other, mv = self._checkOther(other)`
`+`
`+        grSelf = self.grades()`
`+        grOther = other.grades()`
`+`
`+        if len(grSelf) == len(grOther) == 1:`
`+            # both blades`
`+            `
`+            # try the outer product first`
`+            J = self ^ other`
`+`
`+            if J != 0:`
`+                return J.normal()`
`+`
`+            # try something else`
`+            M = (other & self.invPS()).lc(self)`
`+`
`+            if M != 0:`
`+                C = M.normal()`
`+                J = (self & C.rightInv()) ^ other`
`+                return J.normal()`
`+            `
`+            if grSelf[0] >= grOther[0]:`
`+                A = self`
`+                B = other`
`+            else:`
`+                A = other`
`+                B = self`
`+`
`+            if (A & B) == (A * B):`
`+                # B is a subspace of A or the same if grades are equal`
`+                return A.normal()`
`+`
`+            # ugly, but general way`
`+            # watch out for residues`
`+`
`+            # A is still the larger-dimensional subspace`
`+`
`+            Bbasis = B.basis()`
`+`
`+            # add the basis vectors of B one by one to the larger `
`+            # subspace except for the ones that make the outer`
`+            # product vanish`
`+`
`+            J = A`
`+`
`+            for ei in Bbasis:`
`+                J.clean()`
`+                J2 = J ^ ei`
`+`
`+                if J2 != 0:`
`+                    J = J2`
`+            `
`+            # for consistency's sake, we'll normalize the join`
`+            J = J.normal()`
`+                    `
`+            return J`
`+`
`+        else:`
`+            raise ValueError, "not blades"`
`+`
`+    def meet(self, other, subspace=None):`
`+        """Returns the meet of two blades.`
`+`
`+        Computation is done with respect to a subspace that defaults to `
`+        the join 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"`
`+`
`+        if subspace is None:`
`+            subspace = self.join(other)`
`+`
`+        return (self & subspace.inv()) * other`
`+`
`+`
`+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]`
`+`
`+        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:`
`+                marker = -2`
`+                tmp = curBladeX[:]  # copy`
`+                tmp.reverse()`
`+                `
`+                # locate where the steady increase begins`
`+                for j in range(k-1):`
`+                    if tmp[j] - tmp[j+1] == 1:`
`+                        marker = marker - 1`
`+                    else:`
`+                        break`
`+                        `
`+                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, mvClass=MultiVector):`
`+    """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, mvClass)`
`+`
`+    return layout, blades`
`+`
`+`
`+def bases(layout, mvClass=MultiVector):`
`+    """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]] = mvClass(layout, v)`
`+    return dict`
`+`
`+def randomMV(layout, min=-2.0, max=2.0, grades=None, mvClass=MultiVector):`
`+    """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 mvClass(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 mvClass(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`
`+`
`+`