Set the default algebraic rules for scalars

def _algebraic_rules_scalar():
    """Set the default algebraic rules for scalars"""
    a = wc("a", head=SCALAR_VAL_TYPES)
    b = wc("b", head=SCALAR_VAL_TYPES)
    x = wc("x", head=SCALAR_TYPES)
    y = wc("y", head=SCALAR_TYPES)
    z = wc("z", head=SCALAR_TYPES)

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

    ScalarTimes._binary_rules.update(check_rules_dict([
        ('R001', (
            pattern_head(a, b),
            lambda a, b: a * b)),
        ('R002', (
            pattern_head(x, x),
            lambda x: x**2)),
        ('R003', (
            pattern_head(Zero, x),
            lambda x: Zero)),
        ('R004', (
            pattern_head(x, Zero),
            lambda x: Zero)),
        ('R005', (
            pattern_head(
                pattern(ScalarPower, x, y),
                pattern(ScalarPower, x, z)),
            lambda x, y, z: x**(y+z))),
        ('R006', (
            pattern_head(x, pattern(ScalarPower, x, -1)),
            lambda x: One)),
    ]))

    ScalarPower._rules.update(check_rules_dict([
        ('R001', (
            pattern_head(a, b),
            lambda a, b: a**b)),
        ('R002', (
            pattern_head(x, 0),
            lambda x: One)),
        ('R003', (
            pattern_head(x, 1),
            lambda x: x)),
        ('R004', (
            pattern_head(pattern(ScalarPower, x, y), z),
            lambda x, y, z: x**(y*z))),
    ]))

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

    ScalarIndexedSum._rules.update(check_rules_dict([
        ('R001', (  # sum over zero -> zero
            pattern_head(Zero, indranges__),
            lambda indranges: Zero)),
        ('R002', (  # pull constant prefactor out of sum
            pattern_head(pattern(ScalarTimes, x, y), indranges__),
            lambda x, y, indranges:
                pull_constfactor_from_sum(x, y, indranges))),
    ]))