Simplify the implementation of ReciprocalFrame

This is fewer lines and more comments

galgebra/ga.py

@@ -4,7 +4,7 @@
 import warnings
 import operator
 import copy
-from itertools import combinations
+import itertools
 import functools
 from functools import reduce
 from typing import Tuple, TypeVar, Callable, Dict, Sequence
@@ -790,7 +790,7 @@ def indexes(self) -> GradedTuple[Tuple[int, ...]]:
  """ Index tuples of basis blades """
  basis_indexes = tuple(self.n_range)
  return GradedTuple(
- tuple(combinations(basis_indexes, i))
+ tuple(itertools.combinations(basis_indexes, i))
  for i in range(len(basis_indexes) + 1)
  )
 
@@ -1998,51 +1998,43 @@ def ReciprocalFrame(self, basis: Sequence[_mv.Mv], mode: str = 'norm') -> Tuple[
  Arbitrary strings are interpreted as ``"append"``, but in
  future will be an error
  """
- dim = len(basis)
-
- indexes = tuple(range(dim))
- index = [()]
-
- for i in indexes[-2:]:
- index.append(tuple(combinations(indexes, i + 1)))
-
- MFbasis = []
-
- for igrade in index[-2:]:
- grade = []
- for iblade in igrade:
- blade = self.mv(S(1), 'scalar')
- for ibasis in iblade:
- blade ^= basis[ibasis]
- blade = blade.trigsimp()
- grade.append(blade)
- MFbasis.append(grade)
- E = MFbasis[-1][0]
- E_sq = trigsimp((E * E).scalar())
 
- duals = copy.copy(MFbasis[-2])
+ def wedge_reduce(mvs):
+ """ wedge together a list of multivectors """
+ if not mvs:
+ return self.mv(S(1), 'scalar')
+ return functools.reduce(operator.xor, mvs).trigsimp()
+
+ E = wedge_reduce(basis)
+
+ # elements are such that `basis[i] ^ co_basis[i] == E`
+ co_basis = [
+ sign * wedge_reduce(basis_subset)
+ for sign, basis_subset in zip(
+ # alternating signs
+ itertools.cycle([S(1), S(-1)]),
+ # tuples with one basis missing
+ itertools.combinations(basis, len(basis) - 1),
+ )
+ ]
 
- duals.reverse()
- sgn = S(1)
- rbasis = []
- for dual in duals:
- recpv = (sgn * dual * E).trigsimp()
- rbasis.append(recpv)
- sgn = -sgn
+ # take the dual without normalization
+ r_basis = [(co_base * E).trigsimp() for co_base in co_basis]
 
+ # normalize
+ E_sq = trigsimp((E * E).scalar())
  if mode == 'norm':
- for i in range(dim):
- rbasis[i] = rbasis[i] / E_sq
+ r_basis = [r_base / E_sq for r_base in r_basis]
  else:
  if mode != 'append':
  # galgebra 0.5.0
  warnings.warn(
  "Mode {!r} not understood, falling back to {!r} but this "
  "is deprecated".format(mode, 'append'),
  DeprecationWarning, stacklevel=2)
- rbasis.append(E_sq)
+ r_basis.append(E_sq)
 
- return tuple(rbasis)
+ return tuple(r_basis)
 
  def Mlt(self, *args, **kwargs):
  return lt.Mlt(args[0], self, *args[1:], **kwargs)