def _algebraic_rules_operator():
    """Set the default algebraic rules for the operations defined in this
    module"""
    u = wc("u", head=SCALAR_TYPES)
    v = wc("v", head=SCALAR_TYPES)

    n = wc("n", head=(int, str, SymbolicLabelBase))
    m = wc("m", head=(int, str, SymbolicLabelBase))

    A = wc("A", head=Operator)
    B = wc("B", head=Operator)

    A_plus = wc("A", head=OperatorPlus)
    A_times = wc("A", head=OperatorTimes)

    ls = wc("ls", head=LocalSpace)
    h1 = wc("h1", head=HilbertSpace)
    H_ProductSpace = wc("H", head=ProductSpace)

    localsigma = wc(
        'localsigma', head=LocalSigma, kwargs={'hs': ls})

    ra = wc("ra", head=(int, str, SymbolicLabelBase))
    rb = wc("rb", head=(int, str, SymbolicLabelBase))
    rc = wc("rc", head=(int, str, SymbolicLabelBase))
    rd = wc("rd", head=(int, str, SymbolicLabelBase))

    indranges__ = wc("indranges__", head=IndexRangeBase)

    ScalarTimesOperator._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(1, A),
            lambda A: A)),
        ('R002', (
            pattern_head(0, A),
            lambda A: ZeroOperator)),
        ('R003', (
            pattern_head(u, ZeroOperator),
            lambda u: ZeroOperator)),
        ('R004', (
            pattern_head(u, pattern(ScalarTimesOperator, v, A)),
            lambda u, v, A: (u * v) * A)),
        ('R005', (
            pattern_head(-1, A_plus),
            lambda A: OperatorPlus.create(*[-1 * op for op in A.args]))),
    ]))

    OperatorTimes._binary_rules.update(check_rules_dict([
        ('R001', (
            pattern_head(pattern(ScalarTimesOperator, u, A), B),
            lambda u, A, B: u * (A * B))),

        ('R002', (
            pattern_head(ZeroOperator, B),
            lambda B: ZeroOperator)),
        ('R003', (
            pattern_head(A, ZeroOperator),
            lambda A: ZeroOperator)),

        ('R004', (
            pattern_head(A, pattern(ScalarTimesOperator, u, B)),
            lambda A, u, B: u * (A * B))),

        ('R005', (
            pattern_head(
                pattern(LocalSigma, ra, rb, hs=ls),
                pattern(LocalSigma, rc, rd, hs=ls)),
            lambda ls, ra, rb, rc, rd: (
                KroneckerDelta(
                    BasisKet(rb, hs=ls).index,
                    BasisKet(rc, hs=ls).index) *
                LocalSigma.create(ra, rd, hs=ls)))),

        # Harmonic oscillator rules
        ('R009', (
            pattern_head(pattern(Create, hs=ls), localsigma),
            lambda ls, localsigma:
                sqrt(localsigma.index_j + 1) * localsigma.raise_jk(j_incr=1))),

        ('R010', (
            pattern_head(pattern(Destroy, hs=ls), localsigma),
            lambda ls, localsigma:
                sqrt(localsigma.index_j) * localsigma.raise_jk(j_incr=-1))),

        ('R011', (
            pattern_head(localsigma, pattern(Destroy, hs=ls)),
            lambda ls, localsigma:
                sqrt(localsigma.index_k + 1) * localsigma.raise_jk(k_incr=1))),

        ('R012', (
            pattern_head(localsigma, pattern(Create, hs=ls)),
            lambda ls, localsigma:
                sqrt(localsigma.index_k) * localsigma.raise_jk(k_incr=-1))),

        # Normal ordering for harmonic oscillator <=> all a^* to the left, a to
        # the right.
        ('R013', (
            pattern_head(pattern(Destroy, hs=ls), pattern(Create, hs=ls)),
            lambda ls: IdentityOperator + Create(hs=ls) * Destroy(hs=ls))),

        # Oscillator unitary group rules
        ('R014', (
            pattern_head(pattern(Phase, u, hs=ls), pattern(Phase, v, hs=ls)),
            lambda ls, u, v: Phase.create(u + v, hs=ls))),
        ('R015', (
            pattern_head(
                pattern(Displace, u, hs=ls),
                pattern(Displace, v, hs=ls)),
            lambda ls, u, v: (
                exp((u * v.conjugate() - u.conjugate() * v) / 2) *
                Displace.create(u + v, hs=ls)))),

        ('R016', (
            pattern_head(pattern(Destroy, hs=ls), pattern(Phase, u, hs=ls)),
            lambda ls, u:
                exp(I * u) * Phase.create(u, hs=ls) * Destroy(hs=ls))),
        ('R017', (
            pattern_head(pattern(Destroy, hs=ls), pattern(Displace, u, hs=ls)),
            lambda ls, u: Displace.create(u, hs=ls) * (Destroy(hs=ls) + u))),

        ('R018', (
            pattern_head(pattern(Phase, u, hs=ls), pattern(Create, hs=ls)),
            lambda ls, u:
                exp(I * u) * Create(hs=ls) * Phase.create(u, hs=ls))),
        ('R019', (
            pattern_head(pattern(Displace, u, hs=ls), pattern(Create, hs=ls)),
            lambda ls, u: (((Create(hs=ls) - u.conjugate()) *
                            Displace.create(u, hs=ls))))),

        ('R020', (
            pattern_head(pattern(Phase, u, hs=ls), localsigma),
            lambda ls, u, localsigma:
            exp(I * u * localsigma.index_j) * localsigma)),
        ('R021', (
            pattern_head(localsigma, pattern(Phase, u, hs=ls)),
            lambda ls, u, localsigma:
            exp(I * u * localsigma.index_k) * localsigma)),

        # Spin rules
        ('R022', (
            pattern_head(pattern(Jplus, hs=ls), localsigma),
            lambda ls, localsigma:
                Jpjmcoeff(ls, localsigma.index_j, shift=True) *
                localsigma.raise_jk(j_incr=1))),

        ('R023', (
            pattern_head(pattern(Jminus, hs=ls), localsigma),
            lambda ls, localsigma:
                Jmjmcoeff(ls, localsigma.index_j, shift=True) *
                localsigma.raise_jk(j_incr=-1))),

        ('R024', (
            pattern_head(pattern(Jz, hs=ls), localsigma),
            lambda ls, localsigma:
                Jzjmcoeff(ls, localsigma.index_j, shift=True) * localsigma)),

        ('R025', (
            pattern_head(localsigma, pattern(Jplus, hs=ls)),
            lambda ls, localsigma:
                Jmjmcoeff(ls, localsigma.index_k, shift=True) *
                localsigma.raise_jk(k_incr=-1))),

        ('R026', (
            pattern_head(localsigma, pattern(Jminus, hs=ls)),
            lambda ls, localsigma:
                Jpjmcoeff(ls, localsigma.index_k, shift=True) *
                localsigma.raise_jk(k_incr=+1))),

        ('R027', (
            pattern_head(localsigma, pattern(Jz, hs=ls)),
            lambda ls, localsigma:
                Jzjmcoeff(ls, localsigma.index_k, shift=True) * localsigma)),

        # Normal ordering for angular momentum <=> all J_+ to the left, J_z to
        # center and J_- to the right
        ('R028', (
            pattern_head(pattern(Jminus, hs=ls), pattern(Jplus, hs=ls)),
            lambda ls: -2 * Jz(hs=ls) + Jplus(hs=ls) * Jminus(hs=ls))),

        ('R029', (
            pattern_head(pattern(Jminus, hs=ls), pattern(Jz, hs=ls)),
            lambda ls: Jz(hs=ls) * Jminus(hs=ls) + Jminus(hs=ls))),

        ('R030', (
            pattern_head(pattern(Jz, hs=ls), pattern(Jplus, hs=ls)),
            lambda ls: Jplus(hs=ls) * Jz(hs=ls) + Jplus(hs=ls))),
    ]))

    Displace._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(0, hs=ls), lambda ls: IdentityOperator))
    ]))
    Phase._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(0, hs=ls), lambda ls: IdentityOperator))
    ]))
    Squeeze._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(0, hs=ls), lambda ls: IdentityOperator))
    ]))

    OperatorTrace._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(A, over_space=TrivialSpace),
            lambda A: A)),
        ('R002', (
            pattern_head(ZeroOperator, over_space=h1),
            lambda h1: ZeroOperator)),
        ('R003', (
            pattern_head(IdentityOperator, over_space=h1),
            lambda h1: h1.dimension * IdentityOperator)),
        ('R004', (
            pattern_head(A_plus, over_space=h1),
            lambda h1, A: OperatorPlus.create(
                *[OperatorTrace.create(o, over_space=h1)
                  for o in A.operands]))),
        ('R005', (
            pattern_head(pattern(Adjoint, A), over_space=h1),
            lambda h1, A: Adjoint.create(
                OperatorTrace.create(A, over_space=h1)))),
        ('R006', (
            pattern_head(pattern(ScalarTimesOperator, u, A), over_space=h1),
            lambda h1, u, A: u * OperatorTrace.create(A, over_space=h1))),
        ('R007', (
            pattern_head(A, over_space=H_ProductSpace),
            lambda H, A: decompose_space(H, A))),
        ('R008', (
            pattern_head(pattern(Create, hs=ls), over_space=ls),
            lambda ls: ZeroOperator)),
        ('R009', (
            pattern_head(pattern(Destroy, hs=ls), over_space=ls),
            lambda ls: ZeroOperator)),
        ('R010', (
            pattern_head(pattern(LocalSigma, n, m, hs=ls), over_space=ls),
            lambda ls, n, m:
                KroneckerDelta(
                    BasisKet(n, hs=ls).index,
                    BasisKet(m, hs=ls).index) *
                IdentityOperator)),
        ('R011', (
            pattern_head(A, over_space=ls),
            lambda ls, A: factor_for_trace(ls, A))),
    ]))

    Commutator._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(A, A), lambda A: ZeroOperator)),
        ('R002', (
            pattern_head(
                pattern(ScalarTimesOperator, u, A),
                pattern(ScalarTimesOperator, v, B)),
            lambda u, v, A, B: u * v * Commutator.create(A, B))),
        ('R003', (
            pattern_head(pattern(ScalarTimesOperator, v, A), B),
            lambda v, A, B: v * Commutator.create(A, B))),
        ('R004', (
            pattern_head(A, pattern(ScalarTimesOperator, v, B)),
            lambda v, A, B: v * Commutator.create(A, B))),

        # special known commutators
        ('R005', (
            pattern_head(pattern(Create, hs=ls), pattern(Destroy, hs=ls)),
            lambda ls: ScalarTimesOperator(-1, IdentityOperator))),
        # the remaining  rules basically defer to OperatorTimes; just writing
        # out the commutator will generate something simple
        ('R006', (
            pattern_head(
                wc('A', head=(
                    Create, Destroy, LocalSigma, Phase, Displace)),
                wc('B', head=(
                    Create, Destroy, LocalSigma, Phase, Displace))),
            lambda A, B: A * B - B * A)),
        ('R007', (
            pattern_head(
                wc('A', head=(LocalSigma, Jplus, Jminus, Jz)),
                wc('B', head=(LocalSigma, Jplus, Jminus, Jz))),
            lambda A, B: A * B - B * A)),
    ]))

    def pull_constfactor_from_sum(u, A, indranges):
        bound_symbols = set([r.index_symbol for r in indranges])
        if len(u.free_symbols.intersection(bound_symbols)) == 0:
            return u * OperatorIndexedSum.create(A, *indranges)
        else:
            raise CannotSimplify()

    OperatorIndexedSum._rules.update(check_rules_dict([
        ('R001', (  # sum over zero -> zero
            pattern_head(ZeroOperator, indranges__),
            lambda indranges: ZeroOperator)),
        ('R002', (  # pull constant prefactor out of sum
            pattern_head(pattern(ScalarTimesOperator, u, A), indranges__),
            lambda u, A, indranges:
                pull_constfactor_from_sum(u, A, indranges))),
    ]))
Set the default algebraic rules for the operations defined in this
    module