A product for fermionic variables

Question

I'm trying to write a dot product that can handle fermionic variables, i.e., let $a,b$ be fermionic variables, $a\, b=-b\, a$.

There is already a package that can handle fermionic+bosonic variables, grassmann.m, which you can find here
http://people.brandeis.edu/~headrick/Mathematica/grassmann.m
But this package is very slow and for long expressions it stops working (add 12 fermionic variables and it will stop working properly).

What I have written is the following:

ClearAll[dot, fermion, boson, grading]
dot[a___, b_Plus, c___] := dot[a, #, c] & /@ b
dot[a___, b_ d___, c___] := dot[a, b, d, c]
grading[_] := 0;
fermion[a__] := (grading[#] = 1) & /@ {a};
boson[a__] := (grading[#] = 2) & /@ {a};
dot[a___, dot[b___], c___] := dot[a, b, c]
dot[b___, a_, c___, a_, d___] /; OddQ[grading[a]] := 0;
dot[a___, b_ /; grading[b] == 0, c___] := b dot[a, c]
dot[a___, b_, c_, d___] /; (! OrderedQ[{b, c}]) := 
dot[a, c, b, d] (-1)^(grading[b]*grading[c])

Explaining a little bit: I have different gradings: for fermion (1), boson (2) and functions (0). I also want to impose dot to have the attribute Flat, but if I try to add

SetAttributes[dot, Flat]

I end up with an infinite loop when sorting (Question 1: why?, I don't get it. It starts setting the head of any element inside dot a head dot!)

The first two lines are obvious: I want dot to be distributive with respect to Plus, and I found out this to be the simplest way; the second one might be unnecessary, and/or replaced by

dot[a___,b_*c___,d___]/;grading[b]==0:=b dot[a,c,d]

But I haven't tried this option and so far it's working, so I don't dare to change the code! Line 8 plus line 3 do the trick.

Line 7 is a known identity.

Finally, in line 8 I try to sort the list inside dot. This is the hard part. If I impose Flat as an attribute of dot, I end up with an infinite loop (Question 1.) I haven't tried further (I made it finally work two hours ago), but I'm wondering (Question 2) if there is a more efficient way to sort the elements. With grassmann.m sometimes my laptop takes several minutes (15 - 20) to compute something. I don't want that to happen here.

If you look up the package grassmann.m, you will see that most of my code comes from there, except from this last line (and the distributive property.) Also, since grassmann.m already uses NonCummutativeMultiple, which has the Flat attribute, it sometimes conflicts with certain pattern matching replacements, that's why I tried to define my own variables.

As I mentioned, if you have too many fermionic variables, like 12 (6 Weyls in four dimensions), grassmann.m stops assigning grading to some of them, thus, not doing what I want it to do. (Question 3: how can I avoid such problem?, it is related to the memory leak pointed out in one of "Jeremy's additions"?)

Finally, sometimes I will need to generalize the grading in order to add a "form" grading, so this little code should be a good start.

I've looked at other questions, like
Boson commutation relations
Sorting function for non commuting bosons
But I don't understand the answers, and even less how to apply them to my problem. But maybe they are relevant.

Answer 1

After a lot of time working, I've come up with the conclusion that my original code works like a charm.

Defining a dot product as,

ClearAll[dot, fermion, boson, grading]
dot[a___, b_Plus, c___] := dot[a, #, c] & /@ b
dot[a___, b_ d___, c___] := dot[a, b, d, c]
grading[_] := 0;
fermion[a__] := (grading[#] = 1) & /@ {a};
boson[a__] := (grading[#] = 2) & /@ {a};
dot[a___, dot[b___], c___] := dot[a, b, c]
dot[b___, a_, c___, a_, d___] /; OddQ[grading[a]] := 0;
dot[a___, b_ /; EvenQ[grading[b]], c___] := b dot[a, c]
dot[a___, b_, c_, d___] /; (! OrderedQ[{b, c}]) := 
dot[a, c, b, d] (-1)^(grading[b]*grading[c])

is quite efficient. Note that I define both fermions and bosons, but, at the end of the day, only the fermions remain in the noncommutative (nc) product. The reason is simple: I can have $x^n$ for an x that is a boson, and the fastest way to work with this kind of guys is just using a normal product.

Flatten sometimes leads to a loop in certain pattern matching, which I haven't encounter using this dot function. The memory leak in grassmann.m happens because grassmann.m is keeping track of the grading of the expressions. After several computations, grassmann.m forgets the grading and treat expressions as a 0-grading term.

Finally, you can generalize this dot for different grading objects (like some $\mathbb{Z}_n$ grading, as the one found in certain supercosets) and different gradings, such as form grading (such as de Rham differentials) by replacing the 9th line:

dot[a___, b_ , c___] /; (grading1[b]==0&&grading2[b]==0&&...) :=  b dot[a, c] 

To use this code, define your fermions and bosons as,

fermion[t[a__], tb[a__]];
boson[x[a__]];

t[a__] can have as many indices as you want. It can be some t[i,r,$\alpha$], for position i, r-charge r and Lorentz index $\alpha$. From there you can construct any tensor, scalar, or whatever you want.

I hope this is helpful for people out there struggling whit some $\mathcal N=8$ SUGRA.