How to simplify complicated expression further?

New to Julia
1

Hello everyone,

I am using Pluto and Symbolics.jl for symbolic computations and encountered a challenge in simplifying results. Here’s the code I am working with:

solutions = symbolic_linear_solve([L1 ~ 0, L2 ~ 0, L3 ~ 0, L4 ~ 0, L5 ~ 0, L6 ~ 0, L7 ~ 0], [B, D, E, G, F, H, U_D]);
U_D_sol = simplify(solutions[7])

The output for U_D_sol in Julia is as follows:

\text{U_D_sol} = \frac{ - 1.2800000000000002 \left( - 0.03145728000000003 \mathtt{U_{S}} \lambda - 0.031457280000000025 \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda + 0.09437184000000011 \mathtt{U_{S}} \alpha \lambda + 0.031457280000000025 \lambda^{2} \mathtt{U_{S}} + 0.09437184000000008 \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \alpha \lambda - 0.07077888000000003 \alpha^{2} \mathtt{U_{S}} \lambda - 0.09437184000000007 \lambda^{2} \mathtt{U_{S}} \alpha - 0.07077888000000006 \alpha^{2} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda - 0.0707788800000001 \alpha^{3} \mathtt{U_{S}} \lambda + 0.07077888000000004 \lambda^{2} \alpha^{2} \mathtt{U_{S}} - 0.031457280000000025 \alpha^{3} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda + 0.16515072000000014 \alpha^{4} \mathtt{U_{S}} \lambda + 0.03145728000000002 \lambda^{2} \alpha^{3} \mathtt{U_{S}} + 0.04718592000000002 \alpha^{4} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda - 0.017694720000000025 \alpha^{5} \mathtt{U_{S}} \lambda - 0.047185920000000034 \lambda^{2} \alpha^{4} \mathtt{U_{S}} + 0.03145728000000002 \alpha^{5} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda - 0.2595225600000001 \alpha^{6} \mathtt{U_{S}} \lambda + 0.04718592000000002 \lambda^{2} \alpha^{5} \mathtt{U_{S}} - 0.10223616000000002 \alpha^{6} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda + 0.21823488000000013 \alpha^{7} \mathtt{U_{S}} \lambda - 0.13369344000000005 \lambda^{2} \alpha^{6} \mathtt{U_{S}} + 0.07077888 \alpha^{7} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda + 0.07077888000000004 \alpha^{8} \mathtt{U_{S}} \lambda + 0.10616832000000004 \lambda^{2} \alpha^{7} \mathtt{U_{S}} + 0.031457280000000004 \alpha^{8} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda - 0.11599872000000008 \alpha^{9} \mathtt{U_{S}} \lambda + 0.04718592000000002 \lambda^{2} \alpha^{8} \mathtt{U_{S}} - 0.047185920000000006 \alpha^{9} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda - 0.07077888000000002 \lambda^{2} \alpha^{9} \mathtt{U_{S}} + 0.017694720000000018 \alpha^{11} \mathtt{U_{S}} \lambda + 0.007864320000000003 \alpha^{11} \mathtt{Bq_{\kappa}} \mathtt{U_{S}} \lambda + 0.011796480000000005 \lambda^{2} \alpha^{11} \mathtt{U_{S}} \right)}{ - 0.05120000000000002 \left( 0.7864320000000005 + 0.7864320000000005 \mathtt{Bq_{\kappa}} - 4.128768000000003 \alpha - 1.572864000000001 \lambda - 4.128768000000003 \mathtt{Bq_{\kappa}} \alpha - 0.7864320000000005 \mathtt{Bq_{\kappa}} \lambda + 7.077888000000004 \alpha^{2} + 6.488064000000003 \alpha \lambda + 0.7864320000000004 \lambda^{2} + 7.077888000000003 \alpha^{2} \mathtt{Bq_{\kappa}} + 2.3592960000000014 \mathtt{Bq_{\kappa}} \alpha \lambda - 1.2288000000000008 \alpha^{3} - 8.847360000000005 \alpha^{2} \lambda - 2.359296000000001 \lambda^{2} \alpha - 1.2288000000000006 \alpha^{3} \mathtt{Bq_{\kappa}} - 1.769472000000001 \alpha^{2} \mathtt{Bq_{\kappa}} \lambda - 8.847360000000004 \alpha^{4} + 0.44236800000000015 \alpha^{3} \lambda + 1.7694720000000008 \lambda^{2} \alpha^{2} - 8.847360000000004 \alpha^{4} \mathtt{Bq_{\kappa}} - 0.7864320000000002 \alpha^{3} \mathtt{Bq_{\kappa}} \lambda + 5.308416000000001 \alpha^{5} + 10.027008000000004 \alpha^{4} \lambda + 0.7864320000000002 \lambda^{2} \alpha^{3} + 5.308416000000001 \alpha^{5} \mathtt{Bq_{\kappa}} + 1.1796479999999998 \alpha^{4} \mathtt{Bq_{\kappa}} \lambda + 8.257536000000002 \alpha^{6} - 6.488064000000002 \alpha^{5} \lambda - 1.1796480000000003 \lambda^{2} \alpha^{4} + 8.257536000000002 \alpha^{6} \mathtt{Bq_{\kappa}} + 0.7864320000000002 \alpha^{5} \mathtt{Bq_{\kappa}} \lambda - 9.732096 \alpha^{7} - 2.949120000000003 \alpha^{6} \lambda + 1.1796480000000003 \lambda^{2} \alpha^{5} - 9.732095999999999 \alpha^{7} \mathtt{Bq_{\kappa}} - 2.5559039999999995 \alpha^{6} \mathtt{Bq_{\kappa}} \lambda + 1.179648000000003 \alpha^{7} \lambda - 3.3423360000000004 \lambda^{2} \alpha^{6} + 1.7694719999999997 \alpha^{7} \mathtt{Bq_{\kappa}} \lambda + 3.9321599999999974 \alpha^{9} + 3.244031999999999 \alpha^{8} \lambda + 2.6542079999999997 \lambda^{2} \alpha^{7} + 3.932159999999998 \alpha^{9} \mathtt{Bq_{\kappa}} + 0.786432 \alpha^{8} \mathtt{Bq_{\kappa}} \lambda - 1.179647999999996 \alpha^{10} - 0.1966080000000039 \alpha^{9} \lambda + 1.1796480000000003 \lambda^{2} \alpha^{8} - 1.1796479999999976 \alpha^{10} \mathtt{Bq_{\kappa}} - 1.179648 \alpha^{9} \mathtt{Bq_{\kappa}} \lambda - 0.44236799999999943 \alpha^{11} - 1.7694719999999973 \alpha^{10} \lambda - 1.769472 \lambda^{2} \alpha^{9} - 0.44236799999999965 \alpha^{11} \mathtt{Bq_{\kappa}} + 0.19660799999999845 \alpha^{12} + 0.1474560000000007 \alpha^{11} \lambda + 0.19660799999999917 \alpha^{12} \mathtt{Bq_{\kappa}} + 0.19660800000000003 \alpha^{11} \mathtt{Bq_{\kappa}} \lambda + 0.29491199999999895 \alpha^{12} \lambda + 0.294912 \lambda^{2} \alpha^{11} \right)}

In Mathematica, solving the same system yields a much simpler form:

\text{U_D_Mathematica} = \frac{\left( 4 \left( -1 + \lambda \right) - 5 \alpha^{3} + 4 \mathtt{Bq\_\kappa} \left( -1 + \alpha^{5} \right) + \alpha^{5} \left( 9 + 6 \lambda \right) \right) \mathtt{U\_S} \lambda}{9 \alpha \left( -1 + \lambda \right) + 4 \left( -1 + \lambda \right)^{2} - 10 \alpha^{3} \left( -1 + \lambda \right) + \mathtt{Bq\_\kappa} \left( 4 \left( -1 + \alpha^{5} \right) \lambda + \left( -1 + \alpha \right)^{4} \left( 4 + \alpha \left( 7 + 4 \alpha \right) \right) \right) + \alpha^{6} \left( 4 + 6 \lambda \right) + 3 \alpha^{5} \left( -3 + \lambda + 2 \lambda^{2} \right)}

I checked that both are equivalent using:

simplify(U_D_sol-U_D_Mathematica; expand=true) ~ 0

However, how can I simplify Julia’s solution, U_D_sol, to the form shown by Mathematica’s solution?

I’ve tried expand before simplify and it made it worse.
I also used Rational types for the coefficients of the system of equations. It helped simplify the solution since the coefficients were Int type, but then when I started manipulating the solution further, I encountered an OverflowError error.
Are there additional techniques in Symbolics.jl, such as specific simplification rules or settings, that could help simplify U_D_sol further?

Any suggestions are greatly appreciated! Thank you!

2

Hi, and welcome to the Julia community!

Could you provide the full code, in particular definitions of L1, B , etc.?

3

Thanks a lot for helping me out.

Here is a working example to reproduce the previously mention results:

using Symbolics, Nemo
@variables U_S U_D a λ α Bq_κ
@variables B D E F G H
L1 = 2*D - α^2*a^2*B - α^3*a^3*U_D
L2 = 2*H - α^2*a^2*F + 2*α^3*a^3*G - α^3*a^3*U_D + 1//5*α^5*a^5*E
L3 = 2*H - a^2*F + 2*a^3*G - a^3*U_S + 1//5*a^5*E
L4 = H + 1//2*a^2*F - 2*a^3*G + a^3*U_S - 2//5*a^5*E
L5 =  - D + H - 1//2*α^2*a^2*B + 1//2*α^2*a^2*F - 2*α^3*a^3*G - 2//5*α^5*a^5*E
L6 = - 6*D - 6*Bq_κ*D + 6*H*λ + α^2*a^2*B*Bq_κ + 3//5*α^5*a^5*E*λ
L7 = - B + F*λ
solutions = symbolic_linear_solve([L1 ~ 0, L2 ~ 0, L3 ~ 0, L4 ~ 0, L5 ~ 0, L6 ~ 0, L7 ~ 0], [B, D, E, G, F, H, U_D]);
U_D_sol = simplify(solutions[7])

u_D_numerator_Mathematica = λ*(-5*α^3 + 4*(-1 + λ) + α^5*(9 + 6*λ) + 4*(-1 + α^5)*Bq_κ)*U_S
u_D_denominator_Mathematica = (9*α*(-1 + λ) - 10*α^3*(-1 + λ) + 4*(-1 + λ)^2 + α^6*(4 + 6*λ) + 3*α^5*(-3 + λ + 2*λ^2) + ( (-1 + α)^4*(4 + α*(7 + 4*α)) + 4*(-1 + α^5)*λ)*Bq_κ )
u_D_Mathematica = u_D_numerator_Mathematica // u_D_denominator_Mathematica
simplify(U_D_sol-u_D_Mathematica; expand=true)

I hope this helps. Please let me know if you need any other information.

4

I did not manage to find an elegant solution, but you can simplify the expression using some custom rules.

There seem to be two problems.

  • Symbolics.jl does not cancel out common Rational factors.

    julia> @variables x y z;

julia> simplify((2//1 * x + 4//1 * y) / (6//1 * x + 8//1 * z))
((2//1)*x + (4//1)*y) / ((6//1)*x + (8//1)*z)

julia> simplify((2 * x + 4 * y) / (6 * x + 8 * z))
(x + 2y) / (3x + 4z)

    This seems like a bug or just something that still needs to be implemented. In any case, in our situation all Rationals in U_D_sol have denominator 1, so we can ‘convert’ them to Ints, as a workaround.

    rdiv1 = @rule ~x => (~x).num where (~x isa Rational && (~x).den == 1)

  • Mathematica gives a nice factorised result, but this is not the case for Symbolics.jl. In my opinion this is not a fault of the package, but simply a matter of taste with respect to what is considered ‘simplified’. If we want more factorisation, we will have to add factorisation rules. This turns out to be quite difficult. For example, I struggled to come up with a rule which simplifies 2a - 2b to 2(a - b). The difficulty here is that 2a - 2b seems to be represented as 2 \cdot a + (-2) \cdot b, so there is no factor 2 to match. As another example, (~t)^(~n) will not match for n = 1, as in that case you do not write the exponent. In any case, after a lot of effort, I came up with the factorisation rules below. They almost certainly contain redundancies (especially given the use of @acrule and ~~s), but every time I leave something out, an example breaks.

Code

Note that this will take a while to run as we have 512 (+1) rules.

using Symbolics
using SymbolicUtils.Rewriters
 
@variables U_S U_D a λ α Bq_κ
@variables B D E F G H
L1 = 2*D - α^2*a^2*B - α^3*a^3*U_D
L2 = 2*H - α^2*a^2*F + 2*α^3*a^3*G - α^3*a^3*U_D + 1//5*α^5*a^5*E
L3 = 2*H - a^2*F + 2*a^3*G - a^3*U_S + 1//5*a^5*E
L4 = H + 1//2*a^2*F - 2*a^3*G + a^3*U_S - 2//5*a^5*E
L5 =  - D + H - 1//2*α^2*a^2*B + 1//2*α^2*a^2*F - 2*α^3*a^3*G - 2//5*α^5*a^5*E
L6 = - 6*D - 6*Bq_κ*D + 6*H*λ + α^2*a^2*B*Bq_κ + 3//5*α^5*a^5*E*λ
L7 = - B + F*λ
solutions = symbolic_linear_solve([L1 ~ 0, L2 ~ 0, L3 ~ 0, L4 ~ 0, L5 ~ 0, L6 ~ 0, L7 ~ 0], [B, D, E, G, F, H, U_D]);
U_D_sol = simplify(solutions[7])

println("Symbolics, original:\n", U_D_sol, "\n")

rdiv1 = @rule ~x => (~x).num where (~x isa Rational && (~x).den == 1)  # n//1 => n
s = Postwalk(PassThrough(rdiv1))(U_D_sol) |> simplify  # Not the same as simplify(U_D_sol, rdiv1 |> PassThrough |> Postwalk), which somehow just returns U_D_sol
println("Symbolics, after simplifying coefficients:\n", s, "\n")

function generate_factorisation_rules()
	rules = []
	for left_sign in (true, false),  # + or - ('nicest' case first)
		right_sign in (true, false), 
		left_power in (false, true), # 1 (which is not written) or >= 2 (/ <= -1)
		right_power in (false, true),
		include_x in (false, true),
		include_y in (false, true),
		include_a in (false, true),
		include_b in (false, true),
		s_t_equal in (true, false)
		
		# Create all variations of
		#  ~~x * (~t)^(~n) * ~~y + ~~a * (~t)^(~m) * ~~b => *(~~x..., (~t)^(~n - min(~n, ~m)), ~~y...) + *(~~a..., (~t)^(~m - min(~n, ~m)), ~~b...)) * (~t)^min(~n, ~m)
		# But use ~s and ~t and require them to be equal or opposite: necessary for 2a - 2b (i.e. 2a + (-2)b) => 2(a - b)
		
		# (It seems only ... * -... (i.e. ... * (-1) * ...) and ... * ... need to be considered, 
		#  i.e. left_sign == false is not needed.)
		
		left_expr_factors = []
		left_sign || push!(left_expr_factors, :(-1))
		include_x && push!(left_expr_factors, :(~~x))
		
		if left_power
			push!(left_expr_factors, :((~s)^(~n)))
			n = :(~n)
		else
			push!(left_expr_factors, :(~s))
			n = 1
		end
		include_y && push!(left_expr_factors, :(~~y))
		if length(left_expr_factors) >= 2
			left_expr = :(*($(left_expr_factors...)))
		else
			left_expr = first(left_expr_factors)  # avoid unitary *
		end
		
		# Completely similar for the original second term
		# (You could probably unify the code for left and right.)
		right_expr_factors = []
		right_sign || push!(right_expr_factors, :(-1))
		include_a && push!(right_expr_factors, :(~~a))
		if right_power
			push!(right_expr_factors, :((~t)^(~m)))
			m = :(~m)
		else
			push!(right_expr_factors, :(~t))
			m = 1
		end
		include_b && push!(right_expr_factors, :(~~b))
		if length(right_expr_factors) >= 2
			right_expr = :(*($(right_expr_factors...)))
		else
			right_expr = first(right_expr_factors)
		end
		
		new_left_expr_factors = []
		left_sign || push!(new_left_expr_factors, :(-1))
		include_x && push!(new_left_expr_factors, :(~~x...))
		push!(new_left_expr_factors, :((~s)^($n - min($n, $m))))
		include_y && push!(new_left_expr_factors, :(~~y...))
		if length(new_left_expr_factors) >= 2
			new_left_expr = :(*($(new_left_expr_factors...)))
		else
			new_left_expr = first(new_left_expr_factors)
		end
		
		new_right_expr_factors = []
		right_sign || push!(new_right_expr_factors, :(-1))
		include_a && push!(new_right_expr_factors, :(~~a...))
		push!(new_right_expr_factors, :((~t)^($m - min($n, $m))))
		s_t_equal || push!(new_right_expr_factors, :((-1)^min($n, $m)))
		include_b && push!(new_right_expr_factors, :(~~b...))
		if length(new_right_expr_factors) >= 2
			new_right_expr = :(*($(new_right_expr_factors...)))
		else
			new_right_expr = first(new_right_expr_factors)
		end
		
		new_common_factor_expr = :((~s)^min($n, $m))
		predicate = s_t_equal ? :(~t === ~s) : :(~t === -~s)
		
        rule = :(($left_expr + $right_expr) => (($new_left_expr + $new_right_expr) * $new_common_factor_expr) where $predicate)
        push!(rules, rule)
    end
    return rules
end

factorisation_rules = [eval(:(@acrule $frule)) for frule in generate_factorisation_rules()]	 

s = Fixpoint(Postwalk(Chain(factorisation_rules)))(s)  # Don't simplify, as this just undoes the factorisation
s = Postwalk(PassThrough(rdiv1))(s)  # For some reason both Fixpoint(... rdiv1) and ...(Chain([factorisation_rules; rdiv1])) keep an n//1 
println("Symbolics, after factorisation:\n", s, "\n")

u_D_numerator_Mathematica = λ*(-5*α^3 + 4*(-1 + λ) + α^5*(9 + 6*λ) + 4*(-1 + α^5)*Bq_κ)*U_S
u_D_denominator_Mathematica = (9*α*(-1 + λ) - 10*α^3*(-1 + λ) + 4*(-1 + λ)^2 + α^6*(4 + 6*λ) + 3*α^5*(-3 + λ + 2*λ^2) + ( (-1 + α)^4*(4 + α*(7 + 4*α)) + 4*(-1 + α^5)*λ)*Bq_κ )
u_D_Mathematica = u_D_numerator_Mathematica // u_D_denominator_Mathematica
println("Mathematica:\n", u_D_Mathematica, "\n")

println("Same outcome if output below is 0:")
println(simplify(s - u_D_Mathematica; expand=true))
Output Symbolics, original:

(-(384//1)U_Sλ - (384//1)Bq_κU_Sλ + (384//1)U_S(λ^2) - (480//1)U_S(α^3)λ + (864//1)U_S(α^5)λ + (384//1)Bq_κU_S(α^5)λ + (576//1)U_S(α^5)(λ^2)) / ((384//1) + (384//1)Bq_κ - (864//1)α - (768//1)λ - (864//1)Bq_κα - (384//1)Bq_κλ + (864//1)αλ + (384//1)(λ^2) + (960//1)(α^3) + (960//1)Bq_κ(α^3) - (960//1)(α^3)λ - (864//1)(α^5) - (864//1)Bq_κ(α^5) + (384//1)(α^6) + (288//1)(α^5)λ + (384//1)Bq_κ(α^6) + (384//1)Bq_κ(α^5)λ + (576//1)(α^6)λ + (576//1)(α^5)(λ^2))


Symbolics, after simplifying coefficients:
(-4U_Sλ - 4Bq_κU_Sλ + 4U_S(λ^2) - 5U_S*(α^3)λ + 9U_S(α^5)λ + 4Bq_κU_S*(α^5)λ + 6U_S(α^5)(λ^2)) / (4 + 4Bq_κ - 9α - 8λ - 9Bq_κα - 4Bq_κλ + 9αλ + 4(λ^2) + 10(α^3) + 10Bq_κ*(α^3) - 10(α^3)λ - 9(α^5) - 9Bq_κ(α^5) + 4(α^6) + 3(α^5)λ + 4Bq_κ(α^6) + 4Bq_κ*(α^5)*λ + 6(α^6)λ + 6(α^5)(λ^2))


Symbolics, after factorisation:

((4(-1 + λ + Bq_κ*(-1 + α^5)) + (-5 + (α^2)(9 + 6λ))(α^3))U_Sλ) / (α*(10(1 + Bq_κ - λ)(α^2) - 9(1 - λ + α^4)) + 4(1 + (1 + Bq_κ)(α^6)) + Bq_κ*(4(1 - λ) + α*(-9 + (α^4)(-9 + 4λ))) + (-8 + 4λ + (3 + 6(α + λ))(α^5))*λ)


Mathematica:

((4(-1 + λ) - 5(α^3) + 4Bq_κ*(-1 + α^5) + (α^5)(9 + 6λ))U_Sλ) / (9α(-1 + λ) + 4((-1 + λ)^2) - 10(α^3)(-1 + λ) + Bq_κ(4(-1 + α^5)λ + (4 + α(7 + 4α))((-1 + α)^4)) + (α^6)(4 + 6λ) + 3(α^5)*(-3 + λ + 2(λ^2)))


Same outcome if output below is 0:

0

The final Symbolics.jl output is indeed nicely factorised and more compact than initially. I would still slightly prefer the Mathematica result, but again this is just a matter of taste.

A note on some weird simplify behaviour

For some reason using simplify seems to sometimes simply return the input, ignoring the rewriting that happens internally. Using a rewriter directly might work as you would expect (as in the code above), or not (as in the code below). For example,

julia> using Symbolics, SymbolicUtils.Rewriters

julia> @variables a b c

julia> rdiv1 = @rule ~x => (~x).num where (~x isa Rational && (~x).den == 1)
~x => ((~x).num where ~x isa Rational && (~x).den == 1)

julia> rdiv1(2//1)
2

julia> simplify(2//1 * a + 3//1 * b, Postwalk(PassThrough(rdiv1)))
(2//1)*a + (3//1)*b

julia> Postwalk(PassThrough(rdiv1))(2//1 * a + 3//1 * b)
(2//1)*a + (3//1)*b

If we ‘debug’ we can see that in the simplify case we are actually rewriting to 2a + 3b but for some reason we do not return it. In the direct rewriter case, it seems that no decomposition is taking place.

julia> rprint = @rule ~x => println(~x)
~x => println(~x)

julia> simplify(2//1 * a + 3//1 * b, Postwalk(RestartedChain([rprint, rdiv1])))  # Last line is return value, printed automatically in the REPL (unless ;)
2//1
2
a
2a
3//1
3
b
3b
2a + 3b
(2//1)*a + (3//1)*b

julia> Postwalk(Chain([rdiv1, rprint]))(2//1 * a + 3//1 * b)
(2//1)*a + (3//1)*b
(2//1)*a + (3//1)*b

Honestly, none of this makes sense to me…

2 Likes
5

Thanks a lot for this answer.

The factorization rules you provided are above my current Julia level to understand the details of it, but it is very informative to see how one would go ahead and implement one’s rules. Thanks for that.

6

I’m certainly not an expert, so don’t consider this the definitive way to define rules. rdiv1 is probably fine, but normally you wouldn’t have to bring in metaprogramming like I do in generate_factorisation_rules :slight_smile: .

The idea in this function is that I want to all variants of

@rule ~~x * (~t)^(~n) * ~~y + ~~a * (~t)^(~m) * ~~b => *(~~x..., (~t)^(~n - min(~n, ~m)), ~~y...) + *(~~a..., (~t)^(~m - min(~n, ~m)), ~~b...)) * (~t)^min(~n, ~m)

(as this does not match everything I want it to), so e.g.

@rule ~t * ~~y + (~t)^(~m) => (*(~~y...) + (~t)^(~m - 1)) * (~t)

for t \cdot y + t^m = (y + t^{m - 1}) \cdot t, and

@rule ~~x * ~s + ~~a * ~t => ((*(~~x...) + *(-1, ~~a...)) * (~t)) where ~t === -~s

for x \cdot t - a \cdot t = (x - a) \cdot t. (For reasons of uniformity, I’ll write the former as

@rule ~~x * ~s + (~t)^(~m) => ((*(~~x...) + (~t)^(~m - 1)) * (~s)) where ~t === ~s

.)
The * on the right hand side here is just the n-ary multiplication operator in prefix notation, and ... is for splatting a (single) Vector of arguments into n arguments (e.g. +([1, 2, 3]...) is equivalent with +(1, 2, 3) (though a bit less performant). The Vector in question is the one produced by ~~, which matches zero or more expressions and returns them as a Vector.

To get all variants (I think) we need, we should sometimes drop the ~~x, sometimes the ~m, etc. generate_factorisation_rules just does this programatically by looping over all options, and building up Expressions (i.e. pieces of code, within :( )). At the end we then evaluate them to obtain proper rules. See

for more information.

7

Thanks a lot for the explanation and your help.

1 Like