Provide common subexpression elimination like SymPy's `cse`

[gh-dorp92]

Explanation copied from:

• ​https://www.sympy.org/scipy-2017-codegen-tutorial/notebooks/07-the-hard-way.html

Common Subexpression Elimination

If you look carefully at the expressions in the two matrices you'll see repeated expressions. These are not ideal in the sense that the computer has to repeat the exact same calculation multiple times. For large expressions this can be a major issue. Compilers, such as gcc, can often eliminate common subexpressions on their own when different optimization flags are invoked but for complex expressions the algorithms in some compilers do not do a thorough job or compilation can take an extremely long time. SymPy has tools to perform common subexpression elimination which is both thorough and reasonably efficient. In particular if gcc is run with the lowest optimization setting -O0, cse can give large speedups.

For example if you have two expressions:

a = x*y + 5
b = x*y + 6

you can convert this to these three expressions:

z = x*y
a = z + 5
b = z + 6

and x*y only has to be computed once.

The cse() function in SymPy returns

• the subexpression, z = x*y, and
• the simplified expressions: a = z + 5, b = z + 6.

Here is how it works:

>>> sym.cse?

sub_exprs, simplified_rhs = sym.cse(rhs_of_odes)

for var, expr in sub_exprs:
    sym.pprint(sym.Eq(var, expr))

x₀ = 0.0158⋅y₀⋅y₁
x₁ = -x₀
x₂ = 27600000.0⋅y₀⋅y₄
x₃ = -x₂
x₄ = 24400.0⋅y₂⋅y₄
       2
x₅ = y₀ 
x₆ = 14520000.0⋅x₅
x₇ = -x₆
x₈ = 5.83⋅y₄
x₉ = 20900000.0⋅y₀⋅y₁₁
x₁₀ = -x₉
x₁₁ = 35500000.0⋅y₀⋅y₅
x₁₂ = -x₁₁
x₁₃ = 13600000.0⋅y₀⋅y₆
x₁₄ = -x₁₃
x₁₅ = 22900000.0⋅y₀⋅y₇
x₁₆ = -x₁₅
x₁₇ = y₀⋅y₈
x₁₈ = 13000000.0⋅x₁₇
x₁₉ = -x₁₈
x₂₀ = 13000000.0⋅y₀⋅y₉

Is there something like this in SageMath?