docs: 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))),
      ]))