'Use anticommutation to get canonical form of operators.
Employ associativity of normal ordered product: {ab{cd}} = {abcd}
but note that {ab}{cd} /= {abcd}.
We also employ distributivity: {ab + cd} = {ab} + {cd}.
Canonical form also implies expand() {ab(c+d)} = {abc} + {abd}.'

def __new__(cls, arg):

    arg = sympify(arg)
    arg = arg.expand()
    if arg.is_Add:
        return Add(*[cls(term) for term in arg.args])
    if arg.is_Mul:
        (c_part, seq) = arg.args_cnc()
        if c_part:
            coeff = Mul(*c_part)
            if (not seq):
                return coeff
        else:
            coeff = S.One
        newseq = []
        foundit = False
        for fac in seq:
            if isinstance(fac, NO):
                newseq.extend(fac.args)
                foundit = True
            else:
                newseq.append(fac)
        if foundit:
            return (coeff * cls(Mul(*newseq)))
        if isinstance(seq[0], BosonicOperator):
            raise NotImplementedError
        try:
            (newseq, sign) = _sort_anticommuting_fermions(seq)
        except ViolationOfPauliPrinciple:
            return S.Zero
        if (sign % 2):
            return ((S.NegativeOne * coeff) * cls(Mul(*newseq)))
        elif sign:
            return (coeff * cls(Mul(*newseq)))
        else:
            pass
        if (coeff != S.One):
            return (coeff * cls(Mul(*newseq)))
        return Expr.__new__(cls, Mul(*newseq))
    if isinstance(arg, NO):
        return arg
    return arg
